PROGRAM Zkpol1;
{Vlastnosti polynomialni interpolace (Procedure Polint z Numerical Recipes).
 Zada se rovnomerne n1, n2, n3 bodu jmezi xmin, xmax
 y dano funkci sqrt(x/alf) a spoctou se hodnoty extra a interpolace
 pro x od 2xmin-xmax do 2xmax-xmin a porovnaji se
 hodnoty interpolacniho polynomu se zadanou funkci
 vysledky v souboru ZKPOL1.DAT}
CONST
   np = 30;
   maxbod = 200;
TYPE
   RealArrayNP = ARRAY [1..np] OF real;

PROCEDURE polint(VAR xa,ya: RealArrayNP;
                         n: integer;
                         x: real;
                  VAR y,dy: real);
{
 Pocita hodnotu y interpolacniho polynomu  n-1  stupne v bodu x,
 dy je odhad chyby y
 zadane hodnoty funkce jsou xa[n],ya[n]
}
VAR
   ns,m,i: integer;
   w,hp,ho,dift,dif,den: real;
   c,d: ^RealArrayNP;
BEGIN
   new(c);
   new(d);
   ns := 1;
   dif := abs(x-xa[1]);
   FOR i := 1 TO n DO BEGIN
      dift := abs(x-xa[i]);
      IF dift < dif THEN BEGIN
         ns := i;
         dif := dift
      END;
      c^[i] := ya[i];
      d^[i] := ya[i]
   END;
   y := ya[ns];
   ns := ns-1;
   FOR m := 1 TO n-1 DO BEGIN
      FOR i := 1 TO n-m DO BEGIN
         ho := xa[i]-x;
         hp := xa[i+m]-x;
         w := c^[i+1]-d^[i];
         den := ho-hp;
         IF den = 0.0 THEN BEGIN
            writeln ('pause in routine POLINT');
            readln
         END;
         den := w/den;
         d^[i] := hp*den;
         c^[i] := ho*den
      END;
      IF 2*ns < n-m THEN
         dy := c^[ns+1]
      ELSE BEGIN
         dy := d^[ns];
         ns := ns-1
      END;
      y := y+dy
   END;
   dispose(d);
   dispose(c)
END;

Function Hod(xh,alf:real):Real;
BEGIN
     IF XH > 0 THEN Hod:=sqrt(xh/alf)
               ELSE Hod:=0.
END;

VAR
   XHOD,YHOD:RealArrayNP;
   x,y,xmin,xmax,dx,alf,dy,bodmin,bodmax:Real;
   i,j,k,mbod:Integer;
   xb:Array [1..maxbod] of Real;
   yb:Array [1..maxbod,1..3] of Real;
   n:Array [1..3] of Integer;
   f1 : text;
BEGIN
   ASSIGN(f1,'zkpol1.dat');
   REWRITE(f1);
   Writeln('Zadej xmin,xmax,konstantu alf ve fci');
   Readln(xmin,xmax,alf);
   Writeln(' Zadej stupne polynomu n1,n2,n3, pocet bodu vystupu mbod');
   Readln(n[1],n[2],n[3],mbod);
   bodmin := 2*xmin - xmax;
   bodmax := 2*xmax - xmin;
   dx := (bodmax - bodmin)/(mbod-1);
   FOR i:=1 TO mbod DO xb[i] := bodmin + (i-1)*dx;
   FOR k:=1 TO 3 DO BEGIN
       dx := (xmax-xmin)/n[k];
       FOR i:=0 TO n[k] DO BEGIN
           j := i+1;
           IF i < n[k] THEN xhod[j]:=xmin + i*dx
                       ELSE xhod[j]:=xmax;
           yhod[j] := Hod(xhod[j],alf);
      END;
      FOR i:=1 TO mbod DO
          polint(xhod,yhod,n[k]+1,xb[i],yb[i,k],dy);
      END;
   FOR i:=1 TO mbod DO
      Writeln(f1,xb[i],' ',hod(xb[i],alf),' ',yb[i,1],
                       ' ',yb[i,2],' ',yb[i,3]);
  Close(f1);
END.