PROGRAM d11r5(input,output);
(* driver for routine BALANC *)
CONST
   np=5;
TYPE
   glnp = ARRAY [1..np] OF real;
   glnpnp = ARRAY [1..np,1..np] OF real;
VAR
   i,j : integer;
   a : glnpnp;
   c,r : glnp;

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

BEGIN
   a[1,1] := 1.0; a[1,2] := 100.0; a[1,3] := 1.0;
   a[1,4] := 100.0; a[1,5] := 1.0;
   a[2,1] := 1.0; a[2,2] := 1.0; a[2,3] := 1.0;
   a[2,4] := 1.0; a[2,5] := 1.0;
   a[3,1] := 1.0; a[3,2] := 100.0; a[3,3] := 1.0;
   a[3,4] := 100.0; a[3,5] := 1.0;
   a[4,1] := 1.0; a[4,2] := 1.0; a[4,3] := 1.0;
   a[4,4] := 1.0; a[4,5] := 1.0;
   a[5,1] := 1.0; a[5,2] := 100.0; a[5,3] := 1.0;
   a[5,4] := 100.0; a[5,5] := 1.0;
(* write norms *)
   FOR i := 1 to np DO BEGIN
      r[i] := 0.0;
      c[i] := 0.0;
      FOR j := 1 to np DO BEGIN
         r[i] := r[i]+abs(a[i,j]);
         c[i] := c[i]+abs(a[j,i])
      END
   END;
   writeln('rows:');
   FOR i := 1 to np DO write(r[i]:12:2);
   writeln;
   writeln('columns:');
   FOR i := 1 to np DO write(c[i]:12:2);
   writeln;
   writeln;
   writeln('***** balancing matrix *****');
   writeln;
   balanc(a,np);
(* write norms *)
   FOR i := 1 to np DO BEGIN
      r[i] := 0.0;
      c[i] := 0.0;
      FOR j := 1 to np DO BEGIN
         r[i] := r[i]+abs(a[i,j]);
         c[i] := c[i]+abs(a[j,i])
      END
   END;
   writeln('rows:');
   FOR i := 1 to np DO write(r[i]:12:2);
   writeln;
   writeln('columns:');
   FOR i := 1 to np DO write(c[i]:12:2);
   writeln
END.
PROGRAM d3r10(input,output);
(* driver for routine BCUCOF *)
TYPE
   gl4array = ARRAY [1..4] OF real;
   gl4by4  = ARRAY [1..4,1..4] OF real;
VAR
   d1,d2,ee,x1x2 : real;
   i,j : integer;
   y,y1,y2,y12,x1,x2 : gl4array;
   c : gl4by4;
   wt : ARRAY [1..16,1..16] OF real;
   glflag : boolean;

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

BEGIN
   glflag := true;
   x1[1] := 0.0; x1[2] := 2.0; x1[3] := 2.0; x1[4] := 0.0;
   x2[1] := 0.0; x2[2] := 0.0; x2[3] := 2.0; x2[4] := 2.0;
   d1 := x1[2]-x1[1];
   d2 := x2[4]-x2[1];
   FOR i := 1 to 4 DO BEGIN
      x1x2 := x1[i]*x2[i];
      ee := exp(-x1x2);
      y[i] := x1x2*ee;
      y1[i] := x2[i]*(1.0-x1x2)*ee;
      y2[i] := x1[i]*(1.0-x1x2)*ee;
      y12[i] := (1.0-3.0*x1x2+sqr(x1x2))*ee
   END;
   bcucof(y,y1,y2,y12,d1,d2,c);
   writeln;
   writeln('Coefficients for bicubic interpolation:');
   writeln;
   FOR i := 1 to 4 DO BEGIN
      FOR j := 1 to 4 DO write(c[i,j]:12:6);
      writeln
   END
END.
PROGRAM d3r11(input,output);
(* driver for routine BCUINT *)
TYPE
   gl4array = ARRAY [1..4] OF real;
   gl4by4 = ARRAY [1..4,1..4] OF real;
VAR
   ansy,ansy1,ansy2,ey,ey1,ey2 : real;
   x1,x1l,x1u,x1x2,x2,x2l,x2u,xxyy : real;
   i : integer;
   xx,y,y1,y12,y2,yy : gl4array;
   wt : ARRAY [1..16,1..16] OF real;
   glflag : boolean;

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

(*$I BCUINT.PAS *)

BEGIN
   glflag := true;
   xx[1] := 0.0; xx[2] := 2.0; xx[3] := 2.0; xx[4] := 0.0;
   yy[1] := 0.0; yy[2] := 0.0; yy[3] := 2.0; yy[4] := 2.0;
   x1l := xx[1];
   x1u := xx[2];
   x2l := yy[1];
   x2u := yy[4];
   FOR i := 1 to 4 DO BEGIN
      xxyy := xx[i]*yy[i];
      y[i] := sqr(xxyy);
      y1[i] := 2.0*yy[i]*xxyy;
      y2[i] := 2.0*xx[i]*xxyy;
      y12[i] := 4.0*xxyy
   END;
   writeln;
   writeln('x1':6,'x2':8,'y':7,'expect':11,'y1':6,
      'expect':10,'y2':6,'expect':10);
   writeln;
   FOR i := 1 to 10 DO BEGIN
      x1 := 0.2*i;
      x2 := x1;
      bcuint(y,y1,y2,y12,x1l,x1u,x2l,x2u,x1,x2,ansy,ansy1,ansy2);
      x1x2 := x1*x2;
      ey := sqr(x1x2);
      ey1 := 2.0*x2*x1x2;
      ey2 := 2.0*x1*x1x2;
      writeln(x1:8:4,x2:8:4,ansy:8:4,ey:8:4,
         ansy1:8:4,ey1:8:4,ansy2:8:4,ey2:8:4)
   END
END.
PROGRAM d6r26(input,output,dfile);
(* driver for routine BESSI *)
VAR
   n,i,nval : integer;
   val,x : real;
   txt : string[27];
   dfile : text;

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

(*$I BESSI.PAS *)

BEGIN
   glopen(dfile,'fncval.dat');
   REPEAT readln(dfile,txt) UNTIL (txt = 'Modified Bessel Function In');
   readln(dfile,nval);
   writeln(txt);
   writeln;
   writeln('x':5,'actual':14,'bessi(n,x)':18);
   FOR i := 1 to nval DO BEGIN
      readln(dfile,n,x,val);
      writeln(x:6:2,'   ',val:13,'   ',bessi(n,x):13)
   END;
   close(dfile)
END.

PROGRAM d6r21(input,output,dfile);
(* driver for routine BESSI0 *)
VAR
   i,nval : integer;
   val,x : real;
   txt : string[27];
   dfile : text;

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

BEGIN
   glopen(dfile,'fncval.dat');
   REPEAT readln(dfile,txt) UNTIL (txt = 'Modified Bessel Function I0');
   readln(dfile,nval);
   writeln(txt);
   writeln;
   writeln('x':5,'actual':16,'bessi0(x)':17);
   FOR i := 1 to nval DO BEGIN
      readln(dfile,x,val);
      writeln(x:6:2,val:16:7,bessi0(x):16:7)
   END;
   close(dfile)
END.

PROGRAM d6r23(input,output,dfile);
(* driver for routine BESSI1 *)
VAR
   i,nval : integer;
   val,x : real;
   txt : string[27];
   dfile : text;

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

BEGIN
   glopen(dfile,'fncval.dat');
   REPEAT readln(dfile,txt) UNTIL (txt = 'Modified Bessel Function I1');
   readln(dfile,nval);
   writeln(txt);
   writeln;
   writeln('x':5,'actual':16,'bessi1(x)':17);
   FOR i := 1 to nval DO BEGIN
      readln(dfile,x,val);
      writeln(x:6:2,val:16:7,bessi1(x):16:7)
   END;
   close(dfile)
END.

PROGRAM d6r20(input,output,dfile);
(* driver for routine BESSJ *)
VAR
   i,n,nval : integer;
   val,x : real;
   txt : string[18];
   dfile : text;

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

(*$I BESSJ1.PAS *)

(*$I BESSJ.PAS *)

BEGIN
   glopen(dfile,'fncval.dat');
   REPEAT readln(dfile,txt) UNTIL (txt = 'Bessel Function Jn');
   readln(dfile,nval);
   writeln(txt);
   writeln('n':4,'x':7,'actual':14,'bessj(n,x)':18);
   FOR i := 1 to nval DO BEGIN
      readln(dfile,n,x,val);
      writeln(n:4,x:8:2,'   ',val:13,'   ',bessj(n,x):13)
   END;
   close(dfile)
END.

PROGRAM d6r15(input,output,dfile);
(* driver for routine BESSJ0 *)
VAR
   i,nval : integer;
   val,x : real;
   txt : string[18];
   dfile : text;

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

BEGIN
   glopen(dfile,'fncval.dat');
   REPEAT readln(dfile,txt) UNTIL (txt = 'Bessel Function J0');
   readln(dfile,nval);
   writeln (txt);
   writeln ('x':5,'actual':12,'bessj0(x)':13);
   FOR i := 1 to nval DO BEGIN
      readln(dfile,x,val);
      writeln (x:6:2,val:12:7,bessj0(x):12:7)
   END;
   close(dfile)
END.

PROGRAM d6r17(input,output,dfile);
(* driver for routine BESSJ1 *)
VAR
   i,nval : integer;
   val,x : real;
   txt : string[18];
   dfile : text;

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

BEGIN
   glopen(dfile,'fncval.dat');
   REPEAT readln(dfile,txt) UNTIL (txt = 'Bessel Function J1');
   readln(dfile,nval);
   writeln(txt);
   writeln ('x':5,'actual':12,'bessj1(x)':13);
   FOR i := 1 to nval DO BEGIN
      readln(dfile,x,val);
      writeln(x:6:2,val:12:7,bessj1(x):12:7)
   END;
   close(dfile);
END.

PROGRAM d6r25(input,output,dfile);
(* driver for routine BESSK *)
VAR
   n,i,nval : integer;
   val,x : real;
   txt : string[27];
   dfile : text;

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

(*$I BESSI1.PAS *)

(*$I BESSK0.PAS *)

(*$I BESSK1.PAS *)

(*$I BESSK.PAS *)

BEGIN
   glopen(dfile,'fncval.dat');
   REPEAT readln(dfile,txt) UNTIL (txt = 'Modified Bessel Function Kn');
   readln(dfile,nval);
   writeln(txt);
   writeln;
   writeln('x':5,'actual':14,'bessk(n,x)':18);
   FOR i := 1 to nval DO BEGIN
      readln(dfile,n,x,val);
      writeln(x:6:2,'   ',val:13,'   ',bessk(n,x):13)
   END;
   close(dfile)
END.

PROGRAM d6r22(input,output,dfile);
(* driver for routine BESSK0 *)
VAR
   i,nval : integer;
   val,x : real;
   txt : string[27];
   dfile : text;

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

(*$I BESSK0.PAS *)

BEGIN
   glopen(dfile,'fncval.dat');
   REPEAT readln(dfile,txt) UNTIL (txt = 'Modified Bessel Function K0');
   readln(dfile,nval);
   writeln(txt);
   writeln;
   writeln('x':5,'actual':16,'bessk0(x)':17);
   FOR i := 1 to nval DO BEGIN
      readln(dfile,x,val);
      writeln(x:6:2,val:16:7,bessk0(x):16:7)
   END;
   close(dfile)
END.

PROGRAM d6r24(input,output,dfile);
(* driver for routine BESSK1 *)
VAR
   i,nval : integer;
   val,x : real;
   txt : string[27];
   dfile : text;

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

(*$I BESSK1.PAS *)

BEGIN
   glopen(dfile,'fncval.dat');
   REPEAT readln(dfile,txt) UNTIL (txt = 'Modified Bessel Function K1');
   readln(dfile,nval);
   writeln(txt);
   writeln;
   writeln('x':5,'actual':16,'bessk1(x)':17);
   FOR i := 1 to nval DO BEGIN
      readln(dfile,x,val);
      writeln(x:6:2,val:16:7,bessk1(x):16:7)
   END;
   close(dfile)
END.

PROGRAM d6r19(input,output,dfile);
(* driver for routine BESSY *)
VAR
   i,n,nval : integer;
   val,x : real;
   txt : string[18];
   dfile : text;

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

(*$I BESSJ1.PAS *)

(*$I BESSY0.PAS *)

(*$I BESSY1.PAS *)

(*$I BESSY.PAS *)

BEGIN
   glopen(dfile,'fncval.dat');
   REPEAT readln(dfile,txt) UNTIL (txt = 'Bessel Function Yn');
   readln(dfile,nval);
   writeln(txt);
   writeln('n':4,'x':7,'actual':14,'bessy(n,x)':18);
   FOR i := 1 to nval DO BEGIN
      readln(dfile,n,x,val);
      writeln (n:4,x:8:2,'   ',val:13,'   ',bessy(n,x):13)
   END;
   close(dfile)
END.

PROGRAM d6r16(input,output,dfile);
(* driver for routine BESSY0 *)
VAR
   i,nval : integer;
   val,x : real;
   txt : string[18];
   dfile : text;

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

(*$I BESSY0.PAS *)

BEGIN
   glopen(dfile,'fncval.dat');
   REPEAT readln(dfile,txt) UNTIL (txt = 'Bessel Function Y0');
   readln(dfile,nval);
   writeln (txt);
   writeln ('x':5,'actual':12,'bessy0(x)':13);
   FOR i := 1 to nval DO BEGIN
      readln(dfile,x,val);
      writeln(x:6:2,val:12:7,bessy0(x):12:7)
   END;
   close(dfile)
 END.

PROGRAM d6r18(input,output,dfile);
(* driver for routine BESSY1 *)
VAR
   i,nval : integer;
   val,x : real;
   txt : string[18];
   dfile : text;

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

(*$I BESSY1.PAS *)

BEGIN
   glopen(dfile,'fncval.dat');
   REPEAT readln(dfile,txt) UNTIL (txt = 'Bessel Function Y1');
   readln(dfile,nval);
   writeln(txt);
   writeln('x':5,'actual':12,'bessy1(x)':13);
   FOR i := 1 to nval DO BEGIN
      readln(dfile,x,val);
      writeln(x:6:2,val:12:7,bessy1(x):12:7)
   END;
   close(dfile)
END.
PROGRAM d7r9(input,output);
(* driver for routine BNLDEV *)
LABEL 10,99;
CONST
   n=20;
   npts=1000;
   iscal=200;
   nn=100;
VAR
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;
   i,j,k,idum,klim,llen : integer;
   pp,xm : real;
   dist : ARRAY [0..20] OF real;
   txt : PACKED ARRAY [1..50] OF char;

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

(*$I GAMMLN.PAS *)

(*$I BNLDEV.PAS *)

BEGIN
   idum := -133;
   llen := 50;
10:   FOR j := 0 to 20 DO BEGIN
      dist[j] := 0.0
   END;
   writeln('Mean of binomial distribution (0 to 20) (neg to end)');
   read(xm);
   IF (xm < 0) THEN GOTO 99;
   pp := xm/nn;
   FOR i := 1 to npts DO BEGIN
      j := round(bnldev(pp,nn,idum));
      IF ((j >= 0) AND (j <= 20)) THEN dist[j] := dist[j]+1
   END;
   writeln('x':4,'p(x)':8,'graph:':10);
   FOR j := 0 to 19 DO BEGIN
      FOR k := 1 to llen DO txt[k] := ' ';
      dist[j] := dist[j]/npts;
      klim := round(iscal*dist[j])+1;
      IF (klim > llen) THEN klim := llen;
      FOR k := 1 to klim DO BEGIN
         txt[k] := '*'
      END;
      writeln(j:4,dist[j]:9:4,'   ',txt)
   END;
   GOTO 10;
99:   
END.
PROGRAM d10r3(input,output);
(* driver for routine BRENT *)
CONST
   tol=1.0e-6;
   eql=1.0e-4;
VAR
   ax,bx,cx,fa,fb,fc,xmin,bren : 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 BRENT.PAS *)

BEGIN
   nmin := 0;
   writeln;
   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);
      bren := brent(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 d15r6(input,output);
(* driver for routine BSSTEP *)
CONST
   n=4;
   glimax=11;
   glnmax=10;
   glncol=7;
TYPE
   glyarray = ARRAY [1..n] OF real;
   glnarray = glyarray;
VAR
   glx : ARRAY [1..glimax] OF real;
   gld : ARRAY [1..glnmax,1..glncol] OF real;
   eps,hdid,hnext,htry,x : real;
   i : integer;
   y,dydx,yscal : glyarray;

PROCEDURE derivs(x: real; y: glyarray; VAR dydx: glyarray);
(* Programs using DERIVS must define the type
TYPE
   glyarray = ARRAY [1..4] OF real;
in the main routine. *)
BEGIN
   dydx[1] := -y[2];
   dydx[2] := y[1]-(1.0/x)*y[2];
   dydx[3] := y[2]-(2.0/x)*y[3];
   dydx[4] := y[3]-(3.0/x)*y[4]
END;

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

(*$I BESSJ1.PAS *)

(*$I BESSJ.PAS *)

(*$I MMID.PAS *)

(*$I RZEXTR.PAS *)

(*$I BSSTEP.PAS *)

BEGIN
   x := 1.0;
   y[1] := bessj0(x);
   y[2] := bessj1(x);
   y[3] := bessj(2,x);
   y[4] := bessj(3,x);
   dydx[1] := -y[2];
   dydx[2] := y[1]-y[2];
   dydx[3] := y[2]-2.0*y[3];
   dydx[4] := y[3]-3.0*y[4];
   FOR i := 1 to n DO BEGIN
      yscal[i] := 1.0
   END;
   htry := 1.0;
   writeln;
   writeln('eps':10,'htry':12,'hdid':12,'hnext':12);
   FOR i := 1 to 15 DO BEGIN
      eps := exp(-i);
      bsstep(y,dydx,n,x,htry,eps,yscal,hdid,hnext);
      writeln('   ',eps:11,htry:8:2,hdid:14:6,hnext:12:6)
   END
END.

