PROGRAM d7r7 (input,output);
(* driver for routine GAMDEV *)
LABEL 10,99;
CONST
   n=20;
   npts=1000;
   iscal=200;
   llen=50;
VAR
   i,ia,idum,j,k,klim : integer;
   words : ARRAY [1..50] OF char;
   dist : ARRAY [0..20] OF real;
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;

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

(*$I GAMDEV.PAS *)

BEGIN
   idum := -13;
10:   FOR j := 0 to 20 DO BEGIN
      dist[j] := 0.0
   END;
   writeln('Select order of Gamma distribution (n=1..20), -1 to end --  ');
   readln (ia);
   IF (ia < 0) THEN GOTO 99;
   IF (ia > 20) THEN GOTO 10;
   FOR i := 1 to npts DO BEGIN
      j := trunc(gamdev(ia,idum));
      IF ((j >= 0) AND (j <= 20)) THEN dist[j] := dist[j]+1
   END;
   writeln;
   writeln ('gamma-distribution deviate, order ',
         ia:2,' of ',npts:6,' points');
   writeln ('x':6,'p(x)':7,'graph:':9);
   FOR j := 0 to 19 DO BEGIN
      dist[j] := dist[j]/npts;
      FOR k := 1 to 50 DO BEGIN
         words[k] := ' '
      END;
      klim := trunc(iscal*dist[j]);
      IF (klim > llen) THEN klim := llen;
      FOR k := 1 to klim DO BEGIN
         words[k] := '*'
      END;
      write (j:6,dist[j]:8:4,'  ');
      FOR k := 1 to klim DO BEGIN
         write (words[k])
      END;
      writeln
   END;
   GOTO 10;
99:
END.
PROGRAM d6r1(input,output,dfile);
(* driver for routine GAMMLN *)
CONST
   pi = 3.1415926;
VAR
   i,nval : integer;
   actual,calc,x : real;
   txt : string[14];
   dfile : text;

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

BEGIN
   glopen(dfile,'fncval.dat');
   REPEAT readln(dfile,txt) UNTIL (txt = 'Gamma Function');
   readln(dfile,nval);
   writeln ('log of gamma function:');
   writeln ('x':10,'actual':21,'gammln(x)':22);
   FOR i := 1 to nval DO BEGIN
      readln(dfile,x,actual);
      IF (x > 0.0) THEN BEGIN
         IF (x >= 1.0) THEN BEGIN
            calc := gammln(x)
         END ELSE BEGIN
            calc := gammln(x+1.0)-ln(x)
         END;
         writeln (x:12:2,ln(actual):20:6,calc:20:6)
      END
   END;
   close(dfile)
END.

PROGRAM d6r4(input,output,dfile);
(* driver for routine FACTLN *)
VAR
   i,n,nval : integer;
   val : real;
   txt : string[11];
   dfile : text;
   gla : ARRAY [1..100] OF real;

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

(*$I FACTLN.PAS *)

BEGIN
   FOR i := 1 to 100 DO gla[i] := -1.0;  (* initialize FACTLN *)
   glopen(dfile,'fncval.dat');
   REPEAT readln(dfile,txt) UNTIL (txt = 'N-factorial');
   readln(dfile,nval);
   writeln ('log of n-factorial');
   writeln ('n':6,'actual':19,'factln(n)':21);
   FOR i := 1 to nval DO BEGIN
      readln(dfile,n,val);
      writeln (n:6,ln(val):20:7,factln(n):20:7)
   END
END.

PROGRAM d6r3(input,output,dfile);
(* driver for routine BICO *)
VAR
   binco : real;
   i,k,n,nval : integer;
   txt : string[21];
   dfile : text;
   gla : ARRAY [1..100] OF real;

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

(*$I FACTLN.PAS *)

(*$I BICO.PAS *)

BEGIN
   FOR i := 1 to 100 DO gla[i] := -1.0; (* initialize FACTLN *)
   glopen(dfile,'fncval.dat');
   REPEAT readln(dfile,txt) UNTIL (txt = 'Binomial Coefficients');
   readln(dfile,nval);
   writeln (txt);
   writeln ('n':6,'k':6,'actual':12,'bico(n,k)':12);
   FOR i := 1 to nval DO BEGIN
      readln(dfile,n,k,binco);
      writeln (n:6,k:6,binco:12:0,bico(n,k):12:0)
   END
END.

PROGRAM d6r6(input,output,dfile);
(* driver for routine GAMMP *)
VAR
   a,val,x : real;
   i,nval : integer;
   txt : string[25];
   dfile : text;

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

(*$I GCF.PAS *)

(*$I GSER.PAS *)

(*$I GAMMP.PAS *)

BEGIN
   glopen(dfile,'fncval.dat');
   REPEAT readln(dfile,txt) UNTIL (txt = 'Incomplete Gamma Function');
   readln(dfile,nval);
   writeln(txt);
   writeln('a':4,'x':11,'actual':14,'gammp(a,x)':14);
   FOR i := 1 to nval DO BEGIN
      readln(dfile,a,x,val);
      writeln(a:6:2,x:12:6,val:12:6,gammp(a,x):12:6)
   END;
   close(dfile)
END.

PROGRAM d6r7(input,output,dfile);
(* driver for routine GAMMQ *)
VAR
   a,val,x : real;
   i,nval : integer;
   txt : string[25];
   dfile : text;

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

(*$I GSER.PAS *)

(*$I GCF.PAS *)

(*$I GAMMQ.PAS *)

BEGIN
   glopen(dfile,'fncval.dat');
   REPEAT readln(dfile,txt) UNTIL (txt = 'Incomplete Gamma Function');
   readln(dfile,nval);
   writeln(txt);
   writeln('a':4,'x':11,'actual':14,'gammq(a,x)':14);
   FOR i := 1 to nval DO BEGIN
      readln(dfile,a,x,val);
      writeln(a:6:2,x:12:6,(1.0-val):12:6,gammq(a,x):12:6)
   END;
   close(dfile)
END.

PROGRAM d6r5(input,output,dfile);
(* driver for routine BETA *)
VAR
   i,nval : integer;
   val,w,z : real;
   txt : string[13];
   dfile : text;

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

(*$I BETA.PAS *)

BEGIN
   glopen(dfile,'fncval.dat');
   REPEAT readln(dfile,txt) UNTIL (txt = 'Beta Function');
   readln(dfile,nval);
   writeln(txt);
   writeln('w':5,'z':6,'actual':16,'beta(w,z)':20);
   FOR i := 1 to nval DO BEGIN
      readln(dfile,w,z,val);
      writeln(w:6:2,z:6:2,' ':5,val:13,' ':5,beta(w,z):13)
   END;
   close(dfile)
 END.

PROGRAM d6r13(input,output,dfile);
(* driver for routines BETAI and BETACF *)
VAR
   a,b,val,x : real;
   i,nval : integer;
   txt : string[24];
   dfile : text;

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

(*$I BETACF.PAS *)

(*$I BETAI.PAS *)

BEGIN
   glopen(dfile,'fncval.dat');
   REPEAT readln(dfile,txt) UNTIL (txt = 'Incomplete Beta Function');
   readln(dfile,nval);
   writeln(txt);
   writeln('a':5,'b':10,'x':12,'actual':14,'betai(x)':13);
   FOR i := 1 to nval DO BEGIN
      readln(dfile,a,b,x,val);
      writeln(a:6:2,b:12:6,x:12:6,val:12:6,betai(a,b,x):12:6)
   END;
   close(dfile)
END.
PROGRAM d7r6 (input,output);
(* driver for routine GASDEV *)
CONST
   n=20;
   nover2=10;   (* n/2 *)
   npts=1000;
   iscal=400;
   llen=50;
VAR
   i,idum,iset,j,k,klim : integer;
   gset : real;
   words : ARRAY [1..50] OF char;
   dist : ARRAY [-nover2..nover2] OF real;
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;
   gliset : integer;
   glgset : real;

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

(*$I GASDEV.PAS *)

BEGIN
   gliset := 0;   (* initializes routine gasdev *)
   idum := -13;   (* initializes ran3 *)
   FOR j := -nover2 to nover2 DO BEGIN
      dist[j] := 0.0
   END;
   FOR i := 1 to npts DO BEGIN
      j := round(0.25*n*gasdev(idum));
      IF ((j >= -nover2) AND (j <= nover2)) THEN dist[j] := dist[j]+1
   END;
   writeln ('normally distributed deviate of ',npts:6,' points');
   writeln ('x':5,'p(x)':10,'graph:':9);
   FOR j := -nover2 to nover2 DO BEGIN
      dist[j] := dist[j]/npts;
      FOR k := 1 to 50 DO BEGIN
         words[k] := ' '
      END;
      klim := trunc(iscal*dist[j]);
      IF  (klim > llen) THEN  klim := llen;
      FOR k := 1 to klim DO BEGIN
         words[k] := '*'
      END;
      write (j/(0.25*n):8:4,dist[j]:8:4,'  ');
      FOR k := 1 to 50 DO BEGIN
         write (words[k])
      END;
      writeln
   END
END.
PROGRAM d4r8(input,output);
(* driver for routine GAULEG *)
CONST
   npoint=10;
   x1=0.0;
   x2=1.0;
   x3=10.0;
   nval=10;
TYPE
   double = real;
   darray=ARRAY [1..npoint] OF double;
VAR
   i,j : integer;
   xx : real;
   x,w : darray;

FUNCTION func(x: real): real;
BEGIN
   func := x*exp(-x)
END;

FUNCTION sngl(x : real): real;
BEGIN
   sngl := x
END;

(*$I GAULEG.PAS *)

BEGIN
   gauleg(x1,x2,x,w,npoint);
   writeln;
   writeln('#':2,'x[i]':11,'w[i]':12);
   FOR i := 1 to npoint DO BEGIN
      writeln(i:2,x[i]:12:6,w[i]:12:6)
   END;
(* demonstrate the use of gauleg for an integral *)
   gauleg(x1,x3,x,w,npoint);
   xx := 0.0;
   FOR i := 1 to npoint DO BEGIN
      xx := xx+sngl(w[i])*func(sngl(x[i]))
   END;
   writeln;
   writeln('Integral from GAULEG: ',xx:12:6);
   writeln('Actual value: ',((1.0+x1)*exp(-x1)-(1.0+x3)*exp(-x3)):12:6)
END.
PROGRAM d2r1 (input,output,dfile);
(* driver program for subroutine GAUSSJ *)
LABEL 10,99;
CONST
   np=20;
   mp=20;
TYPE
   glnpbynp = ARRAY [1..np,1..np] OF real;
   glnpbymp = ARRAY [1..np,1..mp] OF real;
   glnp = ARRAY [1..np] of integer;
VAR
   j,k,l,m,n : integer;
   a,ai,u : glnpbynp;
   b,x,t : glnpbymp;
   dfile : text;

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

BEGIN
   glopen(dfile,'matrx1.dat');
10:   readln(dfile);
   readln(dfile);
   readln(dfile,n,m);
   readln(dfile);
   FOR k := 1 to n DO BEGIN
      FOR l := 1 to n-1 DO read(dfile,a[k,l]);
      readln(dfile,a[k,n])
   END;
   readln(dfile);
   FOR l := 1 to m DO BEGIN
      FOR k := 1 to n-1 DO read(dfile,b[k,l]);
      readln(dfile,b[n,l])
   END;
(* save matrices for later testing of results *)
   FOR l := 1 to n DO BEGIN
      FOR k := 1 to n DO BEGIN
         ai[k,l] := a[k,l]
      END;
      FOR k := 1 to m DO BEGIN
         x[l,k] := b[l,k]
      END
   END;
(* invert matrix *)
   gaussj(ai,n,np,x,m,mp);
   writeln;
   writeln('Inverse of matrix a : ');
   FOR k := 1 to n DO BEGIN
      FOR l := 1 to n-1 DO write(ai[k,l]:12:6);
      writeln(ai[k,n]:12:6)
   END;
(* test results *)
(* check inverse *)
   writeln('a times a-inverse (compare with unit matrix)');
   FOR k := 1 to n DO BEGIN
      FOR l := 1 to n DO BEGIN
         u[k,l] := 0.0;
         FOR j := 1 to n DO BEGIN
            u[k,l] := u[k,l]+a[k,j]*ai[j,l]
         END
      END;
      FOR l := 1 to n-1 DO write(u[k,l]:12:6);
      writeln(u[k,n]:12:6)
   END;
(* check vector solutions *)
   writeln;
   writeln('Check the following vectors for equality:');
   writeln('original':20,'matrix*sol''n':15);
   FOR l := 1 to m DO BEGIN
   writeln('vector ',l:2,':');
      FOR k := 1 to n DO BEGIN
         t[k,l] := 0.0;
         FOR j := 1 to n DO BEGIN
            t[k,l] := t[k,l]+a[k,j]*x[j,l]
         END;
      writeln(' ':8,b[k,l]:12:6,t[k,l]:12:6);
      END
   END;
   writeln('***********************************');
   IF eof(dfile) THEN GOTO 99;
   writeln('press RETURN for next problem:');
   readln;
   GOTO 10;
99:   close(dfile)
END.
PROGRAM d6r9(input,output,dfile);
(* driver for routine GCF *)
VAR
   a,gammcf,gln,val,x : real;
   i,nval : integer;
   txt : string[25];
   dfile : text;

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

(*$I GCF.PAS *)

BEGIN
   glopen(dfile,'fncval.dat');
   REPEAT readln(dfile,txt) UNTIL (txt = 'Incomplete Gamma Function');
   readln(dfile,nval);
   writeln(txt);
   writeln('a':4,'x':11,'actual':14,'gcf(a,x)':13,
      'gammln(a)':13,'gln':8);
   FOR i := 1 to nval DO BEGIN
      readln(dfile,a,x,val);
      IF (x >= (a+1.0)) THEN BEGIN
         gcf(a,x,gammcf,gln);
         writeln(a:6:2,x:12:6,(1.0-val):12:6,
            gammcf:12:6,gammln(a):12:6,gln:12:6)
      END
   END;
   close(dfile)
END.
PROGRAM d10r2(input,output);
(* driver for routine GOLDEN *)
CONST
   tol=1.0e-6;
   eql=1.0e-3;
VAR
   ax,bx,cx,fa,fb,fc,xmin,gold : real;
   i,iflag,j,nmin : integer;
   amin : ARRAY [1..20] OF real;

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

(*$I BESSJ1.PAS *)

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

(*$I MNBRAK.PAS *)

(*$I GOLDEN.PAS *)

BEGIN
   nmin := 0;
   writeln('minima of the function bessj0');
   writeln('min. #':10,'x':8,'bessj0(x)':17,'bessj1(x)':12);
   FOR i := 1 to 100 DO BEGIN
      ax := i;
      bx := i+1.0;
      mnbrak(ax,bx,cx,fa,fb,fc);
      gold := golden(ax,bx,cx,tol,xmin);
      IF (nmin = 0) THEN BEGIN
         amin[1] := xmin;
         nmin := 1;
         writeln(nmin:7,xmin:15:6,bessj0(xmin):12:6,
            bessj1(xmin):12:6)
      END ELSE BEGIN
         iflag := 0;
         FOR j := 1 to nmin DO BEGIN
            IF  (abs(xmin-amin[j]) <= eql*xmin)
               THEN iflag := 1;
         END;
         IF (iflag = 0) THEN BEGIN
            nmin := nmin+1;
            amin[nmin] := xmin;
            writeln(nmin:7,xmin:15:6,bessj0(xmin):12:6,
               bessj1(xmin):12:6)
         END
      END
   END
END.
PROGRAM d6r8(input,output,dfile);
(* driver for routine GSER *)
VAR
   a,gamser,gln,val,x : real;
   i,nval : integer;
   txt : string[25];
   dfile : text;

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

(*$I GSER.PAS *)

BEGIN
   glopen(dfile,'fncval.dat');
   REPEAT readln(dfile,txt) UNTIL (txt = 'Incomplete Gamma Function');
   readln(dfile,nval);
   writeln(txt);
   writeln('a':4,'x':11,'actual':14,'gser(a,x)':14,
         'gammln(a)':12,'gln':8);
   FOR i := 1 to nval DO BEGIN
      readln(dfile,a,x,val);
      gser(a,x,gamser,gln);
      writeln(a:6:2,x:12:6,val:12:6,gamser:12:6,
            gammln(a):12:6,gln:12:6)
   END;
   close(dfile)
END.
