PROGRAM d10r8(input,output);
(* driver for routine F1DIM *)
CONST
   ndim=3;
TYPE
   glndim = ARRAY [1..ndim] OF real;
VAR
   i,j,ncom : integer;
   p,xi : glndim;
   pcom,xicom : glndim;

FUNCTION fnc(x : glndim) : real;
(* Programs using routine FNC must define the type
TYPE
   glndim = ARRAY [1..ndim] OF real;
in the main routine, where ndim is the dimension of vector x. *)
VAR
   f : real;
   i : integer;
BEGIN
   f := 0.0;
   FOR i := 1 to 3 DO BEGIN
      f := f+sqr(x[i]-1.0)
   END;
   fnc := f
END;

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

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

(*$I SCRSHO.PAS *)

BEGIN
   p[1] := 0.0; p[2] := 0.0; p[3] := 0.0;
   ncom := ndim;
   writeln;
   writeln('Enter vector direction along which to');
   writeln('plot the function. Minimum is in the');
   writeln('direction 1.0 1.0 1.0 - enter x y z:');
   read(xi[1],xi[2],xi[3]);
FOR i := 1 to 3 DO write(xi[i]);
writeln;
   FOR j := 1 to ndim DO BEGIN
      pcom[j] := p[j];
      xicom[j] := xi[j]
   END;
   scrsho
END.
PROGRAM d6r2(input,output,dfile);
(* driver for routine FACTRL *)
VAR
   actual : real;
   glntop,i,n,nval : integer;
   txt : string[11];
   dfile : text;
   gla : ARRAY [1..33] OF real;

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

(*$I FACTRL.PAS *)

BEGIN
   glntop := 0;   (* initialize FACTRL *)
   gla[1] := 1.0;
   glopen(dfile,'fncval.dat');
   REPEAT readln(dfile,txt) UNTIL (txt = 'N-factorial');
   readln(dfile,nval);
   writeln (txt);
   writeln ('n':6,'actual':19,'factrl(n)':21);
   FOR i := 1 to nval DO BEGIN
      readln(dfile,n,actual);
      IF (actual < (1.0e10)) THEN BEGIN
         writeln (n:6,actual:20:0,factrl(n):20:0)
      END else BEGIN
         writeln (n:6,' ':7,actual:13,
               ' ':7,factrl(n):13)
      END
   END;
   close(dfile)
END.
PROGRAM d14r10(input,output);
(* driver for routine FGAUSS *)
CONST
   npt=3;
   nlin=2;
   na=6;   (* na=3*nlin *)
TYPE
   glnparam = ARRAY [1..na] OF real;
VAR
   e1,e2,f,x,y : real;
   i,j : integer;
   a,dyda,df : glnparam;

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

BEGIN
   a[1] := 3.0; a[2] := 0.2; a[3] := 0.5;
   a[4] := 1.0; a[5] := 0.7; a[6] := 0.3;
   writeln;
   writeln('x':6,'y':8,'dyda1':9,'dyda2':8,'dyda3':8,
      'dyda4':8,'dyda5':8,'dyda6':8);
   FOR i := 1 to npt DO BEGIN
      x := 0.3*i;
      fgauss(x,a,y,dyda,na);
      e1 := exp(-sqr((x-a[2])/a[3]));
      e2 := exp(-sqr((x-a[5])/a[6]));
      f := a[1]*e1+a[4]*e2;
      df[1] := e1;
      df[4] := e2;
      df[2] := a[1]*e1*2.0*(x-a[2])/(a[3]*a[3]);
      df[5] := a[4]*e2*2.0*(x-a[5])/(a[6]*a[6]);
      df[3] := a[1]*e1*2.0*sqr(x-a[2])/(a[3]*a[3]*a[3]);
      df[6] := a[4]*e2*2.0*sqr(x-a[5])/(a[6]*a[6]*a[6]);
      writeln('from FGAUSS');
      write(x:8:4,y:8:4);
      FOR j := 1 to 6 DO write(dyda[j]:8:4);
      writeln;
      writeln('independent calc.');
      write(x:8:4,f:8:4);
      FOR j := 1 to 6 DO write(df[j]:8:4);
      writeln;
      writeln
   END
END.
PROGRAM d14r1(input,output);
(* driver for routine FIT *)
CONST
   npt=100;
   spread=0.5;
TYPE
   glndata = ARRAY [1..npt] OF real;
VAR
   gliset : integer;
   glgset : real;
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;
   a,b,chi2,q,siga,sigb : real;
   i,idum,mwt : integer;
   x,y,sig : glndata;

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

(*$I GSER.PAS *)

(*$I GCF.PAS *)

(*$I GAMMQ.PAS *)

(*$I RAN3.PAS *)

(*$I GASDEV.PAS *)

(*$I FIT.PAS *)

BEGIN
   gliset := 0;
   idum := -117;
   FOR i := 1 to npt DO BEGIN
      x[i] := 0.1*i;
      y[i] := -2.0*x[i]+1.0+spread*gasdev(idum);
      sig[i] := spread
   END;
   FOR mwt := 0 to 1 DO BEGIN
      fit(x,y,npt,sig,mwt,a,b,siga,sigb,chi2,q);
      writeln;
      IF (mwt = 0) THEN BEGIN
         writeln('ignoring standard deviations')
      END ELSE BEGIN
         writeln('including standard deviation')
      END;
      writeln(' ':5,'a  =  ',a:9:6,' ':6,'uncertainty:',siga:9:6);
      writeln(' ':5,'b  =  ',b:9:6,' ':6,'uncertainty:',sigb:9:6);
      writeln(' ':5,'chi-squared: ',chi2:14:6);
      writeln(' ':5,'goodness-of-fit: ',q:10:6)
   END
END.
PROGRAM d12r11(input,output);
(* driver for routine FIXRTS *)
CONST
   npoles=6;
   twonp=12;   (* twonp=2*npoles *)
   twonp2=14;   (* twonp2=twonp+2 *)
TYPE
   gl2array = ARRAY [1..2] OF real;
   glnparray = ARRAY [1..npoles] OF real;
   glcarray = ARRAY [1..twonp2] OF real;
VAR
   i,j : integer;
   dum : real;
   polish : boolean;
   d : glnparray;
   sixth : gl2array;
   zcoef,zeros : glcarray;

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

(*$I ZROOTS.PAS *)

(*$I FIXRTS.PAS *)

BEGIN
   d[1] := 6.0; d[2] := -15.0; d[3] := 20.0;
   d[4] := -15.0; d[5] := 6.0; d[6] := 0.0;
   polish := true;
   (* finding roots of (z-1.0)^6 := 1.0 *)
   (* first write roots *)
   zcoef[2*npoles+1] := 1.0;
   zcoef[2*npoles+2] := 0.0;
   FOR i := npoles DOWNTO 1 DO BEGIN
      zcoef[2*i-1] := -d[npoles+1-i];
      zcoef[2*i] := 0.0
   END;
   zroots(zcoef,npoles,zeros,polish);
   writeln('Roots of (z-1.0)^6 = 1.0');
   writeln('Root':22,'(z-1.0)^6':27);
   FOR i := 1 to npoles DO BEGIN
      sixth[1] := 1.0;
      sixth[2] := 0.0;
      FOR j := 1 to 6 DO BEGIN
         dum := sixth[1];
         sixth[1] := sixth[1]*(zeros[2*i-1]-1.0)
            -sixth[2]*zeros[2*i];
         sixth[2] := dum*zeros[2*i]
            +sixth[2]*(zeros[2*i-1]-1.0)
      END;
      writeln(i:6,zeros[2*i-1]:12:6,zeros[2*i]:12:6,
         sixth[1]:12:6,sixth[2]:12:6)
   END;
   (* now fix them to lie within unit circle *)
   fixrts(d,npoles);
   (* check results *)
   zcoef[2*npoles+1] := 1.0;
   zcoef[2*npoles+2] := 0.0;
   FOR i := npoles DOWNTO 1 DO BEGIN
      zcoef[2*i-1] := -d[npoles+1-i];
      zcoef[2*i] := 0.0
   END;
   zroots(zcoef,npoles,zeros,polish);
   writeln;
   writeln('Roots reflected in unit circle');
   writeln('Root':22,'(z-1.0)^6':27);
   FOR i := 1 to npoles DO BEGIN
      sixth[1] := 1.0;
      sixth[2] := 0.0;
      FOR j := 1 to 6 DO BEGIN
         dum := sixth[1];
         sixth[1] := sixth[1]*(zeros[2*i-1]-1.0)
            -sixth[2]*zeros[2*i];
         sixth[2] := dum*zeros[2*i]
            +sixth[2]*(zeros[2*i-1]-1.0)
      END;
      writeln(i:6,zeros[2*i-1]:12:6,zeros[2*i]:12:6,
         sixth[1]:12:6,sixth[2]:12:6)
   END
END.
PROGRAM d14r7(input,output);
(* driver for routine FLEG *)
CONST
   nval=5;
   dx=0.2;
   npoly=5;
TYPE
   glnlarray = ARRAY [1..npoly] OF real;
VAR
   i,j : integer;
   x : real;
   afunc : glnlarray;

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

(*$I FLEG.PAS *)

BEGIN
   writeln;
   writeln('legendre polynomials':43);
   writeln('n=1':9,'n=2':10,'n=3':10,'n=4':10,'n=5':10);
   FOR i := 1 to nval DO BEGIN
      x := i*dx;
      fleg(x,afunc,npoly);
      writeln('x  := ',x:6:2);
      FOR j := 1 to npoly DO write(afunc[j]:10:4);
      writeln('  routine FLEG');
      FOR j := 1 to npoly DO write(plgndr(j-1,0,x):10:4);
      writeln('  routine PLGNDR');
      writeln
   END
END.
PROGRAM d1r1(input,output);
(* driver for routine FLMOON *)
CONST
   zon=-5.0;
TYPE
   name = PACKED ARRAY [1..13] OF char;
VAR
   timzon,frac,secs : real;
   i,i1,i2,i3,id,im,iy : integer;
   j1,j2,n,nph : integer;
   phase : ARRAY [0..3] OF name;

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

(*$I CALDAT.PAS *)

(*$I FLMOON.PAS *)

BEGIN
   timzon := zon/24.0;
   phase[0] := 'new moon     ';
   phase[1] := 'first quarter';
   phase[2] := 'full moon    ';
   phase[3] := 'last quarter ';
   writeln('date of the next few phases of the moon');
   writeln('enter today''s date (e.g. 1 31 1982)  :  ');
   readln(im,id,iy);
(* approximate number of full moons since january 1900 *)
   n := trunc(12.37*(iy-1900+trunc((im-0.5)/12.0)));
   nph := 2;
   j1 := julday(im,id,iy);
   flmoon(n,nph,j2,frac);
   n := n+trunc((j1-j2)/28.0);
   writeln;
   writeln('date':10,'time(est)':19,'phase':9);
   FOR i := 1 to 20 DO BEGIN
      flmoon(n,nph,j2,frac);
      frac := 24.0*(frac+timzon);
      IF  (frac < 0.0)  THEN BEGIN
         j2 := j2-1;
         frac := frac+24.0
      END;
      IF  (frac > 12.0)  THEN BEGIN
         j2 := j2+1;
         frac := frac-12.0
      END ELSE BEGIN
         frac := frac+12.0
      END;
      i1 := trunc(frac);
      secs := 3600.0*(frac-i1);
      i2 := trunc(secs/60.0);
      i3 := trunc(secs-60*i2);
      caldat(j2,im,id,iy);
      writeln(im:5,id:3,iy:5,
         i1:9,':',i2:2,':',i3:2,' ':5,phase[nph]);
      IF  (nph = 3)  THEN BEGIN
         nph := 0;
         n := n+1
      END ELSE BEGIN
         nph := nph+1
      END
   END
END.
PROGRAM d12r1(input,output);
(* driver for routine FOUR1 *)
CONST
   nn=32;
   nn2=64; (* 2*nn *)
TYPE
   gldarray = ARRAY [1..nn2] OF real;
VAR
   ii,i,isign : integer;
   data,dcmp : gldarray;

PROCEDURE prntft(data : gldarray; nn : integer);
VAR
   ii,mm,n : integer;
BEGIN
   writeln('n':4,'real(n)':13,'imag.(n)':13,'real(N-n)':12,'imag.(N-n)':13);
   writeln (0:4,data[1]:14:6,data[2]:12:6,data[1]:12:6,data[2]:12:6);
   mm := nn DIV 2;
   FOR ii := 1 to mm DO BEGIN
      n := 2*ii+1;
      writeln (((n-1) DIV 2):4,data[n]:14:6,data[n+1]:12:6,
         data[2*nn+2-n]:12:6,data[2*nn+3-n]:12:6)
   END;
   writeln (' press return to continue ...');
   readln
END;

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

BEGIN
   writeln('h(t) := real-valued even-function');
   writeln('h(n) := h(N-n) and real?');
   FOR ii := 1 to nn DO BEGIN
      i := 2*ii-1;
      data[i] := 1.0/(sqr((i-nn-1.0)/nn)+1.0);
      data[i+1] := 0.0
   END;
   isign := 1;
   four1(data,nn,isign);
   prntft(data,nn);
   writeln('h(t) := imaginary-valued even-function');
   writeln('h(n) := h(N-n) and imaginary?');
   FOR ii := 1 to nn DO BEGIN
      i := 2*ii-1;
      data[i+1] := 1.0/(sqr((i-nn-1.0)/nn)+1.0);
      data[i] := 0.0
   END;
   isign := 1;
   four1(data,nn,isign);
   prntft(data,nn);
   writeln('h(t) := real-valued odd-function');
   writeln('h(n) := -h(N-n) and imaginary?');
   FOR ii := 1 to nn DO BEGIN
      i := 2*ii-1;
      data[i] := ((i-nn-1.0)/nn)/(sqr((i-nn-1.0)/nn)+1.0);
      data[i+1] := 0.0
   END;
   data[1] := 0.0;
   isign := 1;
   four1(data,nn,isign);
   prntft(data,nn);
   writeln('h(t) := imaginary-valued odd-function');
   writeln('h(n) := -h(N-n) and real?');
   FOR ii := 1 to nn DO BEGIN
      i := 2*ii-1;
      data[i+1] := ((i-nn-1.0)/nn)/(sqr((i-nn-1.0)/nn)+1.0);
      data[i] := 0.0
   END;
   data[2] := 0.0;
   isign := 1;
   four1(data,nn,isign);
   prntft(data,nn);
(* transform, inverse-transform test *)
   FOR ii := 1 to nn DO BEGIN
      i := 2*ii-1;
      data[i] := 1.0/(sqr(0.5*(i-nn-1.0)/nn)+1.0);
      dcmp[i] := data[i];
      data[i+1] := (0.25*(i-nn-1.0)/nn)
            *exp(-sqr(0.5*(i-nn-1.0)/nn));
      dcmp[i+1] := data[i+1]
   END;
   isign := 1;
   four1(data,nn,isign);
   isign := -1;
   four1(data,nn,isign);
   writeln;
   writeln('double fourier transform:':33,'original data:':23);
   writeln;
   writeln('k':4,'real h(k)':14,'imag h(k)':13,
      'real h(k)':17,'imag h(k)':12);
   FOR ii := 1 to (nn DIV 2) DO BEGIN
      i := 2*ii-1;
      writeln(((i+1) DIV 2):4,dcmp[i]:14:6,dcmp[i+1]:12:6,
         data[i]/nn:17:6,data[i+1]/nn:12:6)
   END
END.
PROGRAM d12r13(input,output);
(* driver for routine FOURN *)
CONST
   ndim=3;
   ndat2=1024;
TYPE
   gldarray = ARRAY [1..ndat2] OF real;
   glnnarray = ARRAY [1..ndim] OF integer;
VAR
   i,isign,j,k,l,ll,ndum : integer;
   data : gldarray;
   nn : glnnarray;

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

BEGIN
   ndum := 2;
   FOR i := 1 to ndim DO BEGIN
      ndum := ndum*2;
      nn[i] := ndum
   END;
   FOR i := 1 to nn[3] DO BEGIN
      FOR j := 1 to nn[2] DO BEGIN
         FOR k := 1 to nn[1] DO BEGIN
            l := k+(j-1)*nn[1]+(i-1)*nn[2]*nn[1];
            ll := 2*l-1;
            data[ll] := ll;
            data[ll+1] := ll+1
         END
      END
   END;
   isign := +1;
   fourn(data,nn,ndim,isign);
   isign := -1;
   writeln('Double 3-dimensional transform');
   writeln;
   writeln('Double transf.':22,'Original data':24,'Ratio':20);
   writeln('real':10,'imag.':13,'real':12,'imag.':13,'real':11,'imag.':13);
   writeln;
   fourn(data,nn,ndim,isign);
   FOR i := 1 to 4 DO BEGIN
      j := 2*i;
      k := 2*j;
      l := k+(j-1)*nn[1]+(i-1)*nn[2]*nn[1];
      ll := 2*l-1;
      writeln(data[ll]:12:2,data[ll+1]:12:2,ll:10,ll+1:12,
         data[ll]/ll:14:2,data[ll+1]/(ll+1):12:2)
   END;
   writeln;
   writeln('The product of transform lengths is:',nn[1]*nn[2]*nn[3]:4)
END.
PROGRAM d14r6(input,output);
(* driver for routine FPOLY *)
CONST
   nval=15;
   dx=0.1;
   npoly=5;
TYPE
   glnparray = ARRAY [1..npoly] OF real;
VAR
   i,j : integer;
   x : real;
   afunc : glnparray;

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

BEGIN
   writeln;
   writeln('powers of x':38);
   writeln('x':8,'x**0':11,'x**1':10,'x**2':10,'x**3':10,'x**4':10);
   FOR i := 1 to nval DO BEGIN
      x := i*dx;
      fpoly(x,afunc,npoly);
      write(x:10:4);
      FOR j := 1 to npoly DO write(afunc[j]:10:4);
      writeln
   END
END.
PROGRAM d10r9(input,output);
(* driver for routine FRPRMN *)
CONST
   ndim=3;
   ftol=1.0e-6;
   pio2=1.5707963;
TYPE
   glnarray = ARRAY [1..ndim] OF real;
   glndim=glnarray;
VAR
   ncom : integer;
   pcom,xicom : glnarray; 
   angl,fret : real;
   iter,k : integer;
   p : glnarray;

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

(*$I BESSJ1.PAS *)

FUNCTION fnc(x: glnarray): real;
BEGIN
   fnc := 1.0-bessj0(x[1]-0.5)*bessj0(x[2]-0.5)*bessj0(x[3]-0.5)
END;

PROCEDURE dfnc(x: glnarray; VAR df: glnarray);
BEGIN
   df[1] := bessj1(x[1]-0.5)*bessj0(x[2]-0.5)*bessj0(x[3]-0.5);
   df[2] := bessj0(x[1]-0.5)*bessj1(x[2]-0.5)*bessj0(x[3]-0.5);
   df[3] := bessj0(x[1]-0.5)*bessj0(x[2]-0.5)*bessj1(x[3]-0.5)
END;

(*$I F1DIM.PAS *)

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

(*$I MNBRAK.PAS *)

(*$I BRENT.PAS *)

(*$I LINMIN.PAS *)

(*$I FRPRMN.PAS *)

BEGIN
   writeln('Program finds the minimum of a function');
   writeln('with different trial starting vectors.');
   writeln('True minimum is (0.5,0.5,0.5)');
   FOR k := 0 to 4 DO BEGIN
      angl := pio2*k/4.0;
      p[1] := 2.0*cos(angl);
      p[2] := 2.0*sin(angl);
      p[3] := 0.0;
      writeln;
      writeln('Starting vector: (',
         p[1]:6:4,',',p[2]:6:4,',',p[3]:6:4,')');
      frprmn(p,ndim,ftol,iter,fret);
      writeln('Iterations:',iter:3);
      writeln('Solution vector: (',
         p[1]:6:4,',',p[2]:6:4,',',p[3]:6:4,')');
      writeln('Func. value at solution',fret:14)
   END
END.
PROGRAM d13r7(input,output);
(* driver for routine FTEST *)
CONST
   npts=1000;   (* npts must be > mpts for FTEST *)
   mpts=500;
   eps=0.04;
   nval=10;
TYPE
   gln1array = ARRAY [1..npts] OF real;
   gln2array = ARRAY [1..mpts] OF real;
   glnparray = gln1array;
VAR
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;
   gliset : integer;
   glgset : real;
   f,factor,prob,vrnce : real;
   i,idum,j : integer;
   data1 : gln1array;
   data2 : gln2array;

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

(*$I GASDEV.PAS *)

(*$I GAMMLN.PAS *)

(*$I BETACF.PAS *)

(*$I BETAI.PAS *)

(*$I AVEVAR.PAS *)

(*$I FTEST.PAS *)

BEGIN
(* generate two gaussian distributions with 
different variances *)
   gliset := 0;
   idum := -144;
   writeln;
   writeln('Variance 1 = ':16,1.0:5:2);
   writeln('Variance 2':13,'Ratio':11,'Probability':16);
   FOR i := 1 to (nval+1) DO BEGIN
      FOR j := 1 to npts DO BEGIN
         data1[j] := gasdev(idum)
      END;
      vrnce := 1.0+(i-1)*eps;
      factor := sqrt(vrnce);
      FOR j := 1 to mpts DO BEGIN
         data2[j] := factor*gasdev(idum)
      END;
      ftest(data1,npts,data2,mpts,f,prob);
      writeln(vrnce:11:4,f:13:4,prob:13:4)
   END
END.
