(*     Interpolation:                                                                                                     *)
(*                 Parametrische Darstellung  ( x= φ(x) und y = ψ(x) )                 *)                                         
(*                       ( für geschlossene Kurven )                                              *)
(*                              1.  Klassische  Interpolation   ( Newton )                     *)
(*                              2.  Periodische kubische Spline-Funktion                   *)

Off[General::spell]  ;      Off[
      General::spell1]

lauf = 1;                         (*  Lauf mit der Nummer lauf  ,  
          alle vorhergehenden Läufe bleiben erhalten ! *)

sch =  "V";

If[sch == "N", 
      If[index ≥ 0 And index ≤ 
          m  ,    
TableForm[
          Table[{ind = PaddedForm[i, 2], "     xp[", ind, "] = ", 
              PaddedForm[xp[i], {3, 2}], 
"    yp[", 
              ind, "] = ", PaddedForm[yp[i], {3, 2}]}, {i, 0, m}], 
          TableSpacing -> {1, 0}]]]

(*  Stützpunkte 0,  .. , m  eingeben        *)

\(\(If[\ 
        sch == "V", {\ \ m = 
            12, 
xp[0] = \(-4.5\); \ \ \ yp[0] = 
            0. ; \ \ \ xp[1] = \(-3.0\); \ \ \ \ \ yp[1] = 
            3. ; \ \ \ \ \ \ \ \ \ \ xp[2] = \(-1.0\); \ \ \ \ \ \ yp[2] = 
            1.5; \ \ \ \ \ xp[3] = 0.0; \ \ \ \ \ yp[3] = 0.0; \ \ xp[4] = \ 
            1. ; \ \ \ \ \ \ \ \ \ \ yp[4] = \(-1.5\); \ \ \ \ \ xp[5] = \ 
            3; \ \ \ \ \ \ \ \ \ yp[5] = \(-3. \); \ \ \ \ \ \ \ \ xp[6] = \ 
            4.5; \ \ \ \ yp[6] = 
            0.0; \ \ xp[7] = \ \ 3; \ \ \ \ \ \ \ \ \ \ \ yp[7] = \ 
            3.0; \ \ \ \ \ \ \ xp[8] = \(+1. \); \ \ \ \ \ \ \ \ yp[8] = \ 
            1.5; 
\ \ \ \ \ \ xp[9] = \ 0.0; \ \ \ \ \ yp[9] = 
            0.0; \ \ xp[10] = \ \(-1.0\); \ \ \ yp[10] = \(-1.5\); \ \ \ xp[
              11] = \(-3.0\); \ \ \ \ \ yp[
              11] = \ \(-3.0\); \ \ \ \ \ \ 
\ \ \ \ \ \ xp[
              12] = \(-4.5\); \ \ \ \ yp[12] = 0.0;\ \ }];\)\)

(*  Vorgegebene Stützpunkte 0,  .. , m       *)

If[sch == "V", 
      TableForm[
        Table[{ind = PaddedForm[i, 2], "     xp[", ind, "] = ", 
            PaddedForm[xp[i], {3, 2}], "    yp[", ind, "] = ", 
            PaddedForm[yp[i], {3, 2}]}, {i, 0, m}], 
        TableSpacing -> {1, 0}]]

index  = 15;                  
  \(If[sch == "L", 
          If[index ≥ 0\ And\ index ≤ 
              m\ , {m = 
                m - 1; 
\ \ \ \ \ \ \ \ \ Do[{\ xp[j] = xp[j + 1], 
                  yp[j] = yp[j + 1]}, {j, index, m}]}]]\);

(*   Stützpunkte   0,..,m ,  nach dem Entfernen von Punkten    *)

If[sch == "L", 
      If[index ≥ 0 And index ≤ 
          m  ,    
TableForm[
          Table[{ind = PaddedForm[i, 2], "     xp[", ind, "] = ", 
              PaddedForm[xp[i], {3, 2}], 
"    yp[", 
              ind, "] = ", PaddedForm[yp[i], {3, 2}]}, {i, 0, m}], 
          TableSpacing -> {1, 0}]]]

index  =  17;                  
If[sch == "E",If[index ≥ 0 And index ≤ m,{m = m + 1;Do[{xp[j], yp[j]}, {j, 0, index - 1}];
   Do[{xp[j] = xp[j - 1], yp[j] = yp[j - 1]}, {j, 
                        m, index + 1, -1}];
       xp[index] =   -4.45                                         ;
        yp[index] =    - 0.5                                       }]];

(*   Stützpunkte   0,..,m ,  nach dem Einfügen von Punkten    *)

If[sch == "E", 
      If[index ≥ 0 And index ≤ 
          m,    
TableForm[
          Table[{ind = PaddedForm[i, 2], "     xp[", ind, "] = ", 
              PaddedForm[xp[i], {3, 2}], 
"    yp[", 
              ind, "] = ", PaddedForm[yp[i], {3, 2}]}, {i, 0, m}], 
          TableSpacing -> {1, 0}]]]

index  =   1 ;                  
    If[sch == "K",If[index ≥ 0 And index ≤ m,{       xp[index] =    -4.45                                        ;
        yp[index] =  0.5                                     }]];

(*   Stützpunkte   0,..,m ,  nach dem Einfügen von Punkten    *)

If[sch == "K", 
      If[index ≥ 0 And index ≤ 
          m,    
TableForm[
          Table[{ind = PaddedForm[i, 2], "     xp[", ind, "] = ", 
              PaddedForm[xp[i], {3, 2}], 
"    yp[", 
              ind, "] = ", PaddedForm[yp[i], {3, 2}]}, {i, 0, m}], 
          TableSpacing -> {1, 0}]]]

nummer[lauf] = m;

tp[0] = 0;

    Do[{Delt[k - 1] = 
          Sqrt[(xp[k] - xp[k - 1])^2 + (yp[k] - yp[k - 1])^2], 
        tp[k] = tp[k - 1] + Delt[k - 1]}, {k, 1, m}]

nd = 400;

Do[{ 
        tach[j] = 
          8. *ArcTan[1]/nd*
            j, 
                                          xach[j] = -3 Cos[tach[j]] - 1.5, 
        yach[j] = 3*Sin[1.5  tach[j]]}, {j, -100, 
        100}]
Do[{ 
        tach[j] = 
          8. *ArcTan[1]/nd*
            j, 
                                          xach[j] = -6 +  3*tach[j], 
        yach[j] = -1.5 xach[j]}, {j, 101, 
        150}]
Do[{ 
        tach[j] = 
          8. *ArcTan[1]/nd*
            j, 
                                          xach[j] = 6  -  3*tach[251 + j], 
        yach[j] = 
          1.5  xach[j]}, {j, -150, -101}]
 Do[{ 
        tach[j] = 
          8. *ArcTan[1]/nd*
            j, 
                                          xach[j] = 
          1.5 +  3*Cos[tach[-250 + j]], 
        yach[j] = 3*Sin[1.5  tach[-250 + j]]}, {j, 151, 360}]

Do[
        tj[j] = tp[0] + j*(tp[m] - tp[0])/nd, {j, 0, nd}];

<< Graphics`Colors`

liste1 = {Red, HotPink, Green, Apricot, Brown, DarkGreen, Cobalt, 
          Brick, Orange, Magenta, IndianRed, ForestGreen};

\(achtplot = 
      ListPlot[Table[{xach[j], yach[j]}, {j, \(-146\), 360}], 
        PlotJoined\  -> True, 
	
        PlotRange -> {{\(-6\), 5}, {\(-4\), 4}}, 	PlotStyle -> Brown, 
        AspectRatio -> 0.5, AxesLabel -> {"> x", "  ^ y"}]\)

\(punktplot = 
      ListPlot[Table[{xp[k], yp[k]}, {k, 0, m}], PlotJoined\  -> \ False, 	
        PlotRange -> {{\(-6\), 5}, {\(-4\), 4}}, 
	PlotStyle -> Blue, 
        Prolog\  -> \ AbsolutePointSize[5], AspectRatio -> 0.5, 
        AxesLabel -> {"> x", "  ^ y"}]\)

\(linienplot = 
      ListPlot[Table[{xp[k], yp[k]}, {k, 0, m}], PlotJoined\  -> \ True, 	
        PlotRange -> {{\(-6\), 5}, {\(-4\), 4}}, 
	PlotStyle -> Blue, 
        Prolog\  -> \ AbsolutePointSize[5], AspectRatio -> 0.5, 
        AxesLabel -> {"> x", "  ^ y"}]\)

Show[achtplot, linienplot, punktplot, 
      Prolog  ->  AbsolutePointSize[5]]

FunkDivDiff[xpform_, 
        ypform_] :=  (
Do[
          DivDiffret[k, 
              1] = (ypform[k + 1] - ypform[k])/(xpform[k + 1] - 
                  xpform[k]), {k, 0, m - 1}]; 
Do[
          Do[DivDiffret[k, 
                j] = (DivDiffret[k + 1, j - 1] - 
                    DivDiffret[k, j - 1])/(xpform[k + j] - 
                    xpform[k]), {k, 0, m - j}], {j, 2, m}]; 
        DivDiffmatret = 
          Table[DivDiffret[0, j], {j, 1, m}]; {DivDiffmatret})

DivDiffxmat = FunkDivDiff[tp, xp] ;
DivDiffx = DivDiffxmat[[1]];                    
DivDiffymat = FunkDivDiff[tp, yp] ;     
DivDiffy = DivDiffymat[[1]];

\(Do[{xnew[j] = xp[0], 
        prox = tj[j] - 
            tp[0], 
Do[{xnew[j] = 
              xnew[j] + prox*DivDiffx[\([i]\)], 
            prox = prox*\((tj[j] - tp[i])\)}, {i, 1, 
            m}], 
ynew[j] = yp[0], 
        proy = tj[j] - tp[0], 
	
        Do[{ynew[j] = ynew[j] + proy*DivDiffy[\([i]\)], 
            proy = proy*\((tj[j] - tp[i])\)}, {i, 1, m}]}, {j, 0, 
        nd}]\)

\(newtonplot = 
      ListPlot[Table[{xnew[j], ynew[j]}, {j, 0, nd}], 
        PlotJoined\  -> \ True, 
	
        PlotRange -> {{\(-6.0\), 5}, {\(-4\), 4}}, 	PlotStyle -> Red, 
        AspectRatio -> 0.5, AxesLabel -> {"> x", "  ^ y"}]\)

Show[achtplot, newtonplot, linienplot, punktplot, 
      Prolog  ->  AbsolutePointSize[5]]

(*           Periodische Spline - 
            Interpolation                         *)

FunkPerSpl[xpform_, ypform_]:= (
                                        
                                               
Du[0] = xpform[1] - xpform[0];   Du[m - 1] = xpform[m] - xpform[m - 1];
Dh[0] = 2*( Du[m - 1] + Du[0] );
    Do[{Du[k] = xpform[k + 1] - xpform[k], 
              Dh[k] = 2 (Du[k - 1] + Du[k]), Dp[k] = 0}, {k, 1, m - 2}];
Dp[m - 2] = Du[m - 2];  Dh[m - 1] = 2*( Du[m - 2] + Du[m - 1]);
                                        
Ch[0] = Sqrt[Dh[0]];   Cp[0] = Du[m - 1]/Ch[0];
      Do[{Cn[k - 1] = Du[k - 1]/Ch[k - 1], 
              Ch[k] = Sqrt[Dh[k] - Cn[k - 1]^2],    Cp[
                  k] = ( Dp[k] - Cp[k - 1]*Cn[k - 1] )/Ch[k]}, {k, 1, 
              m - 2}];
    Csum = 0;  Do[ Csum = Csum + Cp[i]^2, {i, 1, m - 2}];
    Cn[m - 2] = Cp[m - 2];  Ch[m - 1] = Sqrt[ Dh[m - 1] - Csum];
                      
                      
Dv[0] = ypform[1] - ypform[0];       Dv[m - 1] = ypform[m] - ypform[m - 1];
    Dr[0] = 3 (Dv[0]/Du[0] - Dv[m - 1]/Du[m - 1]);
    Z[0] = Dr[0]/Ch[0];
    \(Do[{Dv[k] = ypform[k + 1] - ypform[k], 
		Dr[k] = 
                3 \((Dv[k]/Du[k] - Dv[k - 1]/Du[k - 1])\), 
		Z[
                  k] = \((Dr[k] - Z[k - 1]*Cn[k - 1])\)/
                  Ch[k]}, 
\ \ \ \ \ \ {k, 1, m - 1}]\);
  Zsum = 0;   Do[ Zsum = Zsum + Cp[i]*Z[i], {i, 0, m - 2}];
  Z[m - 1] = ( Dr[m - 1] - Zsum )/Ch[m - 1];
                           
                           
  Cp[m - 2] = 0;
  B[m - 1] = Z[m - 1]/Ch[m - 1];
    Do[
            B[k] = ( Z[k] - B[k + 1]*Cn[k] - Cp[k]*B[m - 1])/Ch[k], {k, 
              m - 2, 0, -1}];
  B[m] = B[0];
                          
\(Do[{A[
                  k] = \((B[k + 1] - B[k])\)/\((3 
                       Du[k])\), 
		\ \ \ \ \ \ \ \ \ \ \ \ \ Cc[k] = 
                Dv[k]/Du[k] - Du[k]*\((B[k + 1] + 2  B[k])\)/3, \ 
              Dc[k] = ypform[k]}, {k, 0, m - 1}]\);
  Aret = Table[A[k],  {k, 0, m - 1}];      Bret = Table[B[k], {k, 0, m}];   
  Cret = Table[Cc[k], {k, 0, m - 1}];    Dret = Table[Dc[k], {k, 0, m - 1}];
{Aret, Bret, Cret, Dret})

ABCDxmat = FunkPerSpl[tp, xp]  ;
Aperx = ABCDxmat[[1]] ;                        
Bperx = ABCDxmat[[2]];                       (* Koeffizienen Bk  f  ür die x-Werte *)
Cperx = ABCDxmat[[3]];                        (* Koeffizienen Ck  f  ür die x-Werte *)
Dperx = ABCDxmat[[4]];                        (* Koeffizienen Dk  f  ür die x-Werte *)

ABCDymat = FunkPerSpl[tp, yp] ;
Apery = ABCDymat[[1]];                          
Bpery = ABCDymat[[2]];                       (* Koeffizienen Bk  f  ür die x-Werte *)
Cpery = ABCDymat[[3]];                        (* Koeffizienen Ck  f  ür die x-Werte *)
Dpery = ABCDymat[[4]] ;                       (* Koeffizienen Dk  f  ür die x-Werte *)

\(Do[{tint = tp[0], 
	
        Do[{Dnt[k] = tp[k + 1] - tp[k], 
            If[\((tj[j] ≥ 
                    tint\ )\)\  ∧ \ \ \((tj[j] ≤ 
                    tint + Dnt[k]\ )\), {knt = k, Break[]}\ , 
              tint = tint + Dnt[k]]}, {k, 0, m - 1}], 
	
        Dntmj = tj[j] - tp[knt], 
\ \ \ \ \ xper[j] = 
          Aperx[\([knt + 1]\)]*Dntmj^3 + Bperx[\([knt + 1]\)]*Dntmj^2 + 
            Cperx[\([knt + 1]\)]*Dntmj + Dperx[\([knt + 1]\)]\ , 
	
        yper[j] = 
          Apery[\([knt + 1]\)]*Dntmj^3 + Bpery[\([knt + 1]\)]*Dntmj^2 + 
            Cpery[\([knt + 1]\)]*Dntmj + \ \ Dpery[\([knt + 1]\)]\ }, 
	{j, 
        0, nd}]\)

\(persplplot = 
      ListPlot[Table[{xper[j], yper[j]}, {j, 0, nd}], 
        PlotJoined\  -> \ True, 
	
        PlotRange -> {{\(-6\), 5}, {\(-4\), 4}}, 	PlotStyle -> Green, 
        AspectRatio -> 0.5, AxesLabel -> {"> x", "  ^ y"}]\)

Show[achtplot, persplplot, linienplot, punktplot, 
      Prolog  ->  AbsolutePointSize[5]]

graf[lauf]  =  
      Show[achtplot, persplplot, newtonplot, linienplot, punktplot, 
        Prolog  ->  AbsolutePointSize[5]]

Do[{Print["       Lauf Nummer ", l , " mit ", 
          nummer[l], " Stützpunkten " ], Show[graf[l]]}, {l, 
        1, lauf}]