(*$B- Aufgabe 6.26 *)
PROGRAM kurve;
USES Graph;   (*******)
TYPE vektor = ARRAY[1..100] OF real;
VAR n,i,j,fall,xh,yh, xg,yg,xt,yt,oldx,oldy, farbe,keine: integer;
    h, si, xmin,xmax,ymin,ymax: real;  s,x,y,xs,ys: vektor;
    f : text; stri: string[20];

(*$I Alg6_4 Algorithmus 6.4: FUNCTION g *)
(*$I plotproz *)
(*$I alg6_5 Algorithmus 6.5: tridia *)
(*$I tridiazwei von Aufgabe 6.21 *)

PROCEDURE ableitungen(fall,n:integer; VAR x,y,ys:vektor);
VAR i:integer; fak: real;  a,b,c,h,d,u,v:vektor;
BEGIN
 FOR i := 1 TO n-1 DO h[i]:=x[i+1]-x[i];
 FOR i := 1 TO n-1 DO d[i]:=(y[i+1]-y[i])/sqr(h[i]);
 CASE fall OF
   1,2,3: (* defekter Spline*)
     BEGIN
       FOR i := 2 TO n-1 DO
       ys[i] := (y[i+1]-y[i-1])/(x[i+1]-x[i-1]);
       CASE fall OF
       1:   BEGIN ys[1]:=d[1]*h[1]; ys[n]:=d[n-1]*h[n-1] END;
       2:   BEGIN
               ys[1]:=1.5*d[1]*h[1]-0.5*ys[2];
               ys[n]:=1.5*d[n-1]*h[n-1]-0.5*ys[n-1];
            END;
       3:   BEGIN
               ys[1]:=(y[2]-y[n-1])/(h[1]+h[n-1]); ys[n]:=ys[1]
            END;
       END;
     END;

   21,22,23: (* echter Spline *)
    BEGIN
      FOR i := 2 TO n-1 DO
      BEGIN
         a[i] := 2*(1/h[i-1]+1/h[i]);
         b[i]:=1/h[i];
         c[i]:=b[i];
         ys[i] := 3*(d[i]+d[i-1]);
      END;
      CASE fall OF
      21:(* natuerlich *)
        BEGIN
           b[1]:=1/h[1]; c[1]:=b[1];a[1]:=2*b[1];
           a[n]:=2/h[n-1];
           ys[1]:=3*d[1]; ys[n]:=3*d[n-1];
           tridia(n,c,a,b,ys);
        END;
      22: (* de Boor *)
        BEGIN
          a[1]:=1/h[1]; c[1]:=a[1]; b[1]:=a[1]+1/h[2];
          a[n]:=1/h[n-1]; c[n-1]:=a[n]+1/h[n-2];
          ys[1]:=2*d[1]+(d[1]+d[2])/(h[1]+h[2])*h[1];
          ys[n]:=2*d[n-1]+(d[n-1]+d[n-2])/
                    (h[n-1]+h[n-2])*h[n-1];
          tridia(n,c,a,b,ys);
        END;
      23:(* periodisch *)
         BEGIN
            b[1]:=1/h[1]; c[1]:=b[1];a[1]:=(2*b[1]+1/h[n-1]);
            a[n-1]:= 2/h[n-2]+1/h[n-1];
            FOR i := 2 TO n-1 DO
            BEGIN u[i]:=0; v[i]:=3*(d[i]+d[i-1]) END;
            u[1]:=1; u[n-1]:=1; v[1]:=3*(d[1]+d[n-1]);
            tridiazwei(n-1,c,a,b,u,v);
            fak := (v[1]+v[n-1])/(u[1]+u[n-1]+h[n-1]);
            FOR i:= 1 TO n-1 DO ys[i]:=v[i]-fak*u[i];
            ys[n]:=ys[1];
         END;

      END;
    END;
 END;
END;

PROCEDURE ini;
BEGIN
   writeln('Inputfilenamen ? Fuer diese Aufgabe: DAT6.26');
   readln(stri); assign(f,stri); reset(f);
   writeln('Punkte x,y: mit CTRL-Z aufhoeren');
   n:=0;
   WHILE  NOT eof(f)DO
   BEGIN
      n := n+1;
      writeln('x,y'); read(f,x[n],y[n])
   END;
   s[1]:=0;
   FOR i:=2 TO n DO
    s[i]:=s[i-1]+sqrt(sqr(x[i]-x[i-1])+sqr(y[i]-y[i-1]));
   FOR i := 1 TO n DO writeln(s[i],x[i],y[i]);
END;

PROCEDURE menue;
BEGIN
 writeln('Fall fuer Ableitungen eingeben :');
 writeln('1:lineare Interpolation ');
 writeln('2:lineare Interpolation, natuerliche RB');
 writeln('3:lineare Interpolation, periodische RB');
 writeln('21:echter Spline, natuerliche RB');
 writeln('22:echter Spline, de Boor RB');
 writeln('23:echter Spline, periodische RB');
 readln(fall);
END;

BEGIN
    ini;  menue;
    ableitungen(fall,n,s,x,xs);
    ableitungen(fall,n,s,y,ys);
    writeln('xmin,xmax,ymin,ymax eingeben');
    readln(xmin,xmax,ymin,ymax);
    plotbegin( xg,yg,xt,yt,farbe,keine);
    achsen(xmin,xmax,ymin,ymax,xg,yg,xt,yt,farbe, keine);
    for i := 1 to n do kreuz(x[i],y[i],5);
    { gotoxy(1,24);  } (*****)
    pl(g(n,0,s,x,xs),g(n,0,s,y,ys),keine);
    FOR i := 1 TO n-1 DO
    BEGIN
       h := (s[i+1]-s[i])/20;
       FOR j := 0 TO 20 DO
       BEGIN
          si := s[i]+j*h;
          pl(g(n,si,s,x,xs), g(n,si,s,y,ys),farbe)
       END;
    END;
    readln
END.