PROGRAM d3r4 (input,output);
(* driver for routine SPLINT *)
LABEL 1;
CONST
   np=10;
   pi=3.1415926;
TYPE
   glnarray = ARRAY [1..np] OF real;
VAR
   i,nfunc : integer;
   f,x,y,yp1,yp2,ypn : real;
   xa,ya,y2 : glnarray;

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

(*$I SPLINT.PAS *)

BEGIN
   FOR nfunc := 1 to 2 DO BEGIN
      IF (nfunc = 1) THEN BEGIN
         writeln;
         writeln ('sine function from 0 to pi');
         FOR i := 1 to np DO BEGIN
            xa[i] := i*pi/np;
            ya[i] := sin(xa[i])
         END;
         yp1 := cos(xa[1]);
         ypn := cos(xa[np])
      END ELSE IF (nfunc = 2) THEN BEGIN
         writeln;
         writeln ('exponential function from 0 to 1');
         FOR i := 1 to np DO BEGIN
            xa[i] := 1.0*i/np;
            ya[i] := exp(xa[i])
         END;
         yp1 := exp(xa[1]);
         ypn := exp(xa[np])
      END ELSE BEGIN
         GOTO 1
      END;
(* call spline to get second derivatives *)
      spline(xa,ya,np,yp1,yp2,y2);
(* call splint FOR interpolations *)
      writeln;
      writeln ('x':9,'f(x)':13,'interpolation':17);
      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;
         splint(xa,ya,y2,np,x,y);
         writeln (x:12:6,f:12:6,y:12:6)
      END;
      writeln;
      writeln ('***********************************');
      writeln ('press RETURN');
      readln
   END;
1:
END.