PROGRAM d5r9(input,output);
(* driver for routine PCSHFT *)
LABEL 10,99;
CONST
   nval=40;
   pio2=1.5707963;
TYPE
   glcarray= ARRAY [1..nval] OF real;
VAR
   a,b,poly,x : real;
   i,j,mval : integer;
   c,d : glcarray;

FUNCTION func(x: real): real;
BEGIN
   func := sqr(x)*(sqr(x)-2.0)*sin(x)
END;

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

(*$I CHEBPC.PAS *)

(*$I PCSHFT.PAS *)

BEGIN
   a := -pio2;
   b := pio2;
   chebft(a,b,c,nval);
10:   writeln;
   writeln('How many terms in Chebyshev evaluation?');
   write('Enter n between 6 and ',nval:2,
         '. (n := 0 to end).  ');
   readln(mval);
   IF ((mval <= 0) OR (mval > nval)) THEN GOTO 99;
   chebpc(c,d,mval);
   pcshft(a,b,d,mval);
(* test shifted polynomial *)
   writeln;
   writeln('x':9,'actual':14,'polynomial':14);
   FOR i := -8 to 8 DO BEGIN
      x := i*pio2/10.0;
      poly := d[mval];
      FOR j := mval-1 DOWNTO 1 DO BEGIN
         poly := poly*x+d[j]
      END;
      writeln(x:12:6,func(x):12:6,poly:12:6)
   END;
   GOTO 10;
99:
END.
PROGRAM d13r15(input,output);
(* driver for routine PEARSN *)
CONST
   n=10;
TYPE
   nparray = ARRAY [1..n] OF real;
VAR
   prob,r,z : real;
   i : integer;
   dose,spore : nparray;

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

(*$I BETACF.PAS *)

(*$I BETAI.PAS *)

(*$I PEARSN.PAS *)

BEGIN
   dose[1] := 56.1; dose[2] := 64.1; dose[3] := 70.0;
   dose[4] := 66.6; dose[5] := 82.0; dose[6] := 91.3;
   dose[7] := 90.0; dose[8] := 99.7; dose[9] := 115.3;
   dose[10] := 110.0;
   spore[1] := 0.11; spore[2] := 0.40; spore[3] := 0.37;
   spore[4] := 0.48; spore[5] := 0.75; spore[6] := 0.66;
   spore[7] := 0.71; spore[8] := 1.20; spore[9] := 1.01;
   spore[10] := 0.95;
   writeln;
   writeln('Effect of Gamma Rays on Man-in-the-Moon Marigolds');
   writeln('Count Rate (cpm)':16,'Pollen Index':23);
   FOR i := 1 to 10 DO BEGIN
      writeln(dose[i]:10:2,spore[i]:25:2)
   END;
   pearsn(dose,spore,10,r,prob,z);
   writeln;
   writeln('PEARSN':30,'Expected':16);
   writeln('Corr. Coeff.',' ':8,r:13,'  ',0.9069586:13);
   writeln('Probability',' ':9,prob:13,'  ',0.2926505e-3:13);
   writeln('Fisher''s z',' ':10,z:13,'  ',1.510110:13);
END.
PROGRAM d6r27(input,output,dfile);
(* driver for routine PLGNDR *)
VAR
   fac,val,x : real;
   i,j,m,n,nval : integer;
   txt : string[20];
   dfile : text;

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

BEGIN
   glopen(dfile,'fncval.dat');
   REPEAT readln(dfile,txt) UNTIL (txt = 'Legendre Polynomials');
   readln(dfile,nval);
   writeln(txt);
   writeln('n':4,'m':4,'x':10,'actual':19,'plgndr(n,m,x)':20);
   FOR i := 1 to nval DO BEGIN
      readln(dfile,n,m,x,val);
      fac := 1.0;
      IF (m > 0) THEN BEGIN
         FOR j := n-m+1 to n+m DO BEGIN
            fac := fac*j
         END
      END;
      fac := 2.0*fac/(2.0*n+1.0);
      val := val*sqrt(fac);
      writeln(n:4,m:4,'   ',x:13,'   ',val:13,'   ',plgndr(n,m,x):13)
   END;
   close(dfile)
END.
PROGRAM d7r8 (input,output);
(* driver for routine POIDEV *)
LABEL 10,99;
CONST
   n = 20;
   npts = 1000;
   iscal = 200;
   llen = 50;
VAR
   i,idum,j,k,klim : integer;
   xm : real;
   dist : ARRAY [0..20] OF real;
   txt : PACKED ARRAY [1..50] OF char;
   gloldm,glsq,glalxm,glg : real;
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;

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

(*$I GAMMLN.PAS *)

(*$I POIDEV.PAS *)

BEGIN
   gloldm := -1.0;   (* initializes routine poidev *)
   idum := -13;
10:   FOR j := 0 to 20 DO BEGIN
      dist[j] := 0.0
   END;
   writeln('Mean of Poisson distribution (x := 0.0 to 20.0); neg. to end');
   readln(xm);
   IF (xm < 0.0) THEN GOTO 99;
   IF (xm > 20.0) THEN GOTO 10;
   FOR i := 1 to npts DO BEGIN
      j := trunc(poidev(xm,idum));
      IF ((j >= 0) AND (j <= 20)) THEN dist[j] := dist[j]+1
   END;
   writeln('Poisson-distributed deviate, mean ',xm:5:2,
         ' of ',npts:6,' points');
   writeln('x':5,'p(x)':8,'graph:':10);
   FOR j := 0 to 19 DO BEGIN
      dist[j] := dist[j]/npts;
      FOR k := 1 to 50 DO BEGIN
         txt[k] := ' '
      END;
      klim := trunc(iscal*dist[j]);
      IF (klim > llen) THEN klim := llen;
      FOR k := 1 to klim DO BEGIN
         txt[k] := '*'
      END;
      writeln (1.0*j:6:2,dist[j]:8:4,'   ',txt)
   END;
   GOTO 10;
99:
END.
PROGRAM d3r7 (input,output);
(* driver for routine POLCOE *)
LABEL 1;
CONST
   np=5;
   pi=3.1415926;
TYPE
   glnarray = ARRAY [1..np] OF real;
VAR
   i,j,nfunc : integer;
   f,sum,x : real;
   coeff,xa,ya : glnarray;

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

BEGIN
   FOR nfunc := 1 to 2 DO BEGIN
      IF (nfunc = 1) THEN BEGIN
         writeln ('sine function from 0 to pi');
         writeln;
         FOR i := 1 to np DO BEGIN
            xa[i] := i*pi/np;
            ya[i] := sin(xa[i])
         END
      END ELSE IF (nfunc = 2) THEN BEGIN
         writeln ('exponential function from 0 to 1');
         writeln;
         FOR i := 1 to np DO BEGIN
            xa[i] := 1.0*i/np;
            ya[i] := exp(xa[i])
         END
      END ELSE BEGIN
         GOTO 1
      END;
      polcoe(xa,ya,np,coeff);
      writeln (' ':2,'coefficients');
      FOR i := 1 to np DO BEGIN
         writeln (coeff[i]:12:6)
      END;
      writeln;
      writeln('x':9,'f(x)':13,'polynomial':15);
      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;
         sum := coeff[np];
         FOR j := np-1 DOWNTO 1 DO BEGIN
            sum := coeff[j]+sum*x
         END;
         writeln (x:12:6,f:12:6,sum:12:6)
      END;
      writeln;
      writeln ('************************************');
      writeln ('press RETURN');
      readln
   END;
1:
END.
PROGRAM d3r8 (input,output);
(* driver for routine POLCOF *)
LABEL 1;
CONST
   np=5;
   pi=3.1415926;
TYPE
   glnarray = ARRAY [1..np] OF real;
VAR
   i,j,nfunc : integer;
   f,sum,x : real;
   coeff,xa,ya : glnarray;

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

(*$I POLCOF.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
      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
      END ELSE BEGIN
         GOTO 1
      END;
      polcof(xa,ya,np,coeff);
      writeln;
      writeln ('  coefficients');
      FOR i := 1 to np DO BEGIN
         writeln (coeff[i]:12:6)
      END;
      writeln;
      writeln ('x':9,'f(x)':13,'polynomial':15);
      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;
         sum := coeff[np];
         FOR j := np-1 DOWNTO 1 DO BEGIN
            sum := coeff[j]+sum*x
         END;
         writeln (x:12:6,f:12:6,sum:12:6)
      END;
      writeln ('************************************');
      writeln ('press RETURN');
      readln
   END;
1:
END.
PROGRAM d5r3(input,output);
(* driver for routine POLDIV *)
(* (x-1)**5/(x+1)**3 *)
CONST
   n=6;
   nv=4;
TYPE
   glnarray = ARRAY [1..n] OF real;
   glnvarray = ARRAY [1..nv] OF real;
VAR
   i : integer;
   u,q,r : glnarray;
   v : glnvarray;

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

BEGIN
   u[1] := -1.0; u[2] := 5.0; u[3] := -10.0;
   u[4] := 10.0; u[5] := -5.0; u[6] := 1.0;
   v[1] := 1.0; v[2] := 3.0; v[3] := 3.0; v[4] := 1.0;
   poldiv(u,n,v,nv,q,r);
   writeln;
   writeln('x^0':10,'x^1':10,'x^2':10,'x^3':10,'x^4':10,'x^5':10);
   writeln;
   writeln('quotient polynomial coefficients:');
   FOR i := 1 to 6 DO write(q[i]:10:2);
   writeln;
   writeln('expected quotient coefficients:');
   writeln(31.0:10:2,-8.0:10:2,1.0:10:2,0.0:10:2,0.0:10:2,0.0:10:2);
   writeln;
   writeln('remainder polynomial coefficients:');
   FOR i := 1 to 4 DO write(r[i]:10:2);
   writeln;
   writeln('expected remainder coefficients:');
   writeln(-32.0:10:2,-80.0:10:2,-80.0:10:2,0.0:10:2)
END.
PROGRAM d3r9 (input,output);
(* driver for routine POLIN2 *)
CONST
   n=5;
   pi=3.1415926;
TYPE
   glmarray = ARRAY [1..n] OF real;
   glnarray = glmarray;
   glmbyn = ARRAY [1..n,1..n] OF real;
VAR
   i,j : integer;
   dy,f,x1,x2,y : real;
   x1a,x2a : glmarray;
   ya : glmbyn;

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

(*$I POLIN2.PAS *)

BEGIN
   FOR i := 1 to n DO BEGIN
      x1a[i] := i*pi/n;
      FOR j := 1 to n DO BEGIN
         x2a[j] := 1.0*j/n;
         ya[i,j] := sin(x1a[i])*exp(x2a[j])
      END
   END;
(* test 2-dimensional interpolation *)
   writeln;
   writeln ('Two dimensional interpolation of sin(x1)exp(x2)');
   writeln;
   writeln ('x1':9,'x2':12,'f(x)':13,'interpolated':16,'error':11);
   FOR i := 1 to 4 DO BEGIN
      x1 := (-0.1+i/5.0)*pi;
      FOR j := 1 to 4 DO BEGIN
         x2 := -0.1+j/5.0;
         f := sin(x1)*exp(x2);
         polin2(x1a,x2a,ya,n,n,x1,x2,y,dy);
         writeln (x1:12:6,x2:12:6,f:12:6,y:12:6,dy:15:6)
      END;
      writeln ('***********************************');
   END
END.
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.
PROGRAM d10r6(input,output);
(* driver for routine POWELL *)
CONST
   ndim=3;
   ftol=1.0e-6;
TYPE
   glnarray = ARRAY [1..ndim] OF real;
   glndim = glnarray;
   glnpbynp = ARRAY [1..ndim,1..ndim] OF real;
   gl3array = glnarray;
VAR
   ncom : integer;
   pcom,xicom : glnarray; 
   fret : real;
   i,iter,j,np : integer;
   p : glnarray;
   xi : glnpbynp;

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

FUNCTION fnc(x: gl3array): real;
(* Programs using FNC must define the type
TYPE
   gl3array = ARRAY [1..3] OF real;
in the main routine. *)
BEGIN
   fnc := 0.5-bessj0(sqr(x[1]-1.0)+sqr(x[2]-2.0)+sqr(x[3]-3.0))
END;

(*$I F1DIM.PAS *)

FUNCTION func(x: real): real;
BEGIN
   func := f1dim(x)
END;

(*$I MNBRAK.PAS *)

(*$I BRENT.PAS *)

(*$I LINMIN.PAS *)

(*$I POWELL.PAS *)

BEGIN
   np := ndim;
   xi[1,1] := 1.0; xi[1,2] := 0.0; xi[1,3] := 0.0;
   xi[2,1] := 0.0; xi[2,2] := 1.0; xi[2,3] := 0.0;
   xi[3,1] := 0.0; xi[3,2] := 0.0; xi[3,3] := 1.0;
   p[1] := 1.5; p[2] := 1.5; p[3] := 2.5;
   powell(p,xi,ndim,np,ftol,iter,fret);
   writeln('Iterations:',iter:3);
   writeln;
   writeln('Minimum found at: ');
   FOR i := 1 to ndim DO write(p[i]:12:6);
   writeln;
   writeln;
   writeln('Minimum function value =',fret:12:6);
   writeln;
   writeln('True minimum of function is at:');
   writeln(1.0:12:6,2.0:12:6,3.0:12:6)
END.
PROGRAM d12r12(input,output);
(* driver for routine PREDIC *)
CONST
   npts=300;
   npoles=10;
   twonp2=22;   (* twonp2=2*npoles+2 *)
   nfut=20;
   pi=3.1415926;
TYPE
   glnarray = ARRAY [1..npts] OF real;
   gldarray = glnarray;
   glmarray = ARRAY [1..npoles] OF real;
   glnparray = glmarray;
   glnfarray = ARRAY [1..nfut] OF real;      
   glcarray = ARRAY [1..twonp2] OF real;
   gl2array = ARRAY [1..2] OF real;
VAR
   i : integer;
   dum : real;
   d,wkm : glmarray;
   data,wk1,wk2 : glnarray;
   future : glnfarray;

FUNCTION f(n,npts: integer): real;
CONST
   pi=3.1415926;
BEGIN
   f := exp(-1.0*n/npts)*sin(2.0*pi*n/50.0)
      +exp(-2.0*n/npts)*sin(2.2*pi*n/50.0)
END;

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

(*$I ZROOTS.PAS *)

(*$I MEMCOF.PAS *)

(*$I FIXRTS.PAS *)

(*$I PREDIC.PAS *)
   
BEGIN
   FOR i := 1 to npts DO BEGIN
      data[i] := f(i,npts)
   END;
   memcof(data,npts,npoles,dum,d,wk1,wk2,wkm);
   fixrts(d,npoles);
   predic(data,npts,d,npoles,future,nfut);
   writeln('I':6,'Actual':11,'PREDIC':12);
   FOR i := 1 to nfut DO BEGIN
      writeln (i:6,f(i+npts,npts):12:6,future[i]:12:6)
   END
END.
PROGRAM d13r12(input,output);
(* driver for routine PROBKS *)
CONST
   npts=20;
   eps=0.1;
   iscal=40;
VAR
   alam,aval : real;
   i,j : integer;
   txt : PACKED ARRAY [1..iscal] OF char;

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

BEGIN
   writeln('probability function for kolmogorov-smirnov statistic');
   writeln;
   writeln('lambda':7,'value:':10,'graph:':13);
   FOR i := 1 to npts DO BEGIN
      alam := i*eps;
      aval := probks(alam);
      FOR j := 1 to iscal DO BEGIN
         IF (j <= round(iscal*aval)) THEN BEGIN
            txt[j] := '*'
         END ELSE BEGIN
            txt[j] := ' '
         END
      END;
      writeln(alam:8:6,aval:10:6,' ':5,txt)
   END
END.
PROGRAM d15r8(input,output);
(* driver for routine PZEXTR *)
CONST
   glimax=11;
   glnmax=10;
   glncol=7;
   nv=4;
   nuse=5;
TYPE
   glyarray = ARRAY [1..nv] OF real;
VAR
   glx : ARRAY [1..glimax] OF real;
   glqcol : ARRAY [1..glnmax,1..glncol] OF real;
   dum,xest : real;
   i,iest,j : integer;
   dy,yest,yz : glyarray;

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

BEGIN
(* feed values from a rational function *)
(* fn(x) := (1-x+x**3)/(x+1)**n *)
   FOR i := 1 to 10 DO BEGIN
      iest := i;
      xest := 1.0/i;
      dum := 1.0-xest+xest*xest*xest;
      FOR j := 1 to nv DO BEGIN
         dum := dum/(xest+1.0);
         yest[j] := dum
      END;
      pzextr(iest,xest,yest,yz,dy,nv,nuse);
      writeln;
      writeln('i  :=  ',i:2);
      write('Extrap. function:');
      FOR j := 1 to nv DO write(yz[j]:12:6);
      writeln;
      write('Estimated error: ');
      FOR j := 1 to nv DO write(dy[j]:12:6);
      writeln
   END;
   writeln;
   writeln('actual values:   ',1.0:12:6,1.0:12:6,1.0:12:6,1.0:12:6)
END.
