PROGRAM d3r1 (input,output);
(* driver for routine POLINT *)
LABEL 1;
CONST
   np=10;   (* maximum value for n *)
   pi=3.1415926;
TYPE
   glnarray = ARRAY [1..np] OF real;
VAR
   i,n,nfunc : integer;
   dy,f,x,y : real;
   xa,ya : glnarray;

(*$I MODFILE.PAS *)
(*$I POLINT.PAS *)

BEGIN
   writeln ('generation of interpolation tables');
   writeln (' ... sin(x)   0<x<pi');
   writeln (' ... exp(x)   0<x<1 ');
   writeln ('how many entries go in these tables? (note: n<10)');
   readln (n);
   FOR nfunc := 1 to 2 DO BEGIN
      IF (nfunc = 1) THEN BEGIN
         writeln;
         writeln ('sine function from 0 to pi');
         FOR i := 1 to n DO BEGIN
            xa[i] := i*pi/n;
            ya[i] := sin(xa[i])
         END;
      END ELSE IF (nfunc = 2) THEN BEGIN
         writeln;
         writeln('exponential function from 0 to 1');
         FOR i := 1 to n DO BEGIN
            xa[i] := i*1.0/n;
            ya[i] := exp(xa[i])
         END;
      END ELSE BEGIN
         GOTO 1
      END;
      writeln;
      writeln ('x':9,'f(x)':13,'interpolated':16,'error':11);
      FOR i := 1 to 10 DO BEGIN
         IF (nfunc = 1) THEN BEGIN
            x := (-0.05+i/10.0)*pi;
            f := sin(x)
         END ELSE IF  (nfunc = 2) THEN BEGIN
            x := (-0.05+i/10.0);
            f := exp(x)
         END;
         polint(xa,ya,n,x,y,dy);
         writeln (x:12:6,f:12:6,y:12:6,' ':4,dy:11)
      END;
      writeln;
      writeln ('***********************************');
      writeln ('press RETURN');
      readln
   END;
1:
END.