PROGRAM d7r1 (input,output);
(* driver for routine RAN0 *)
(* calculates pi statistically using volume of unit n-sphere *)
CONST
   pi=3.1415926;
VAR
   i,j,k,idum,jpower : integer;
   x1,x2,x3,x4 : real;
   iy : ARRAY [1..3] OF integer;
   yprob : ARRAY [1..3] OF real;
   gliseed : integer;
   gly : real;
   glv : ARRAY [1..97] OF real;

FUNCTION fnc(x1,x2,x3,x4 : real) : real;
BEGIN
   fnc := sqrt(sqr(x1)+sqr(x2)+sqr(x3)+sqr(x4))
END;

FUNCTION twotoj(j : integer): integer;
BEGIN
   IF (j=0) THEN twotoj := 1
   ELSE twotoj := 2*twotoj(j-1)
END;

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

BEGIN
   idum := -1;
   FOR i := 1 to 3 DO BEGIN
      iy[i] := 0
   END;
   writeln;
   writeln ('volume of unit n-sphere, n = 2,3,4');
   writeln ('# points','     pi     ','  (4/3)*pi  ',' (1/2)*pi^2 ');
   writeln;
   FOR j := 1 to 13 DO BEGIN
      FOR k := twotoj(j-1) to twotoj(j) DO BEGIN
         x1 := ran0(idum);
         x2 := ran0(idum);
         x3 := ran0(idum);
         x4 := ran0(idum);
         IF (fnc(x1,x2,0.0,0.0) < 1.0) THEN  iy[1] := iy[1]+1;
         IF (fnc(x1,x2,x3,0.0) < 1.0) THEN  iy[2] := iy[2]+1;
         IF (fnc(x1,x2,x3,x4) < 1.0) THEN  iy[3] := iy[3]+1
      END;
      jpower := twotoj(j);
      yprob[1] := 4.0*iy[1]/jpower;
      yprob[2] := 8.0*iy[2]/jpower;
      yprob[3] := 16.0*iy[3]/jpower;
      writeln (jpower:6,yprob[1]:12:6,yprob[2]:12:6,yprob[3]:12:6)
   END;
   writeln;
   writeln ('actual',pi:12:6,(4.0*pi/3.0):12:6,(0.5*sqr(pi)):12:6)
END.

PROGRAM d7r2 (input,output);
(* driver for routine RAN1 *)
(* calculates pi statistically using volume of unit n-sphere *)
CONST
   pi=3.1415926;
VAR
   i,idum,j,k,jpower : integer;
   x1,x2,x3,x4 : real;
   iy : ARRAY [1..3] OF integer;
   yprob : ARRAY [1..3] OF real;
   glix1,glix2,glix3 : integer;
   glr : ARRAY [1..97] OF real;

FUNCTION fnc(x1,x2,x3,x4 : real) : real;
BEGIN
   fnc := sqrt(sqr(x1)+sqr(x2)+sqr(x3)+sqr(x4))
END;

FUNCTION twotoj(j : integer) : integer;
BEGIN
   IF (j=0) THEN twotoj := 1
   ELSE twotoj := 2*twotoj(j-1)
END;

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

BEGIN
   idum := -1;
   FOR i := 1 to 3 DO BEGIN
      iy[i] := 0;
   END;
   writeln;
   writeln ('volume of unit n-sphere, n = 2,3,4');
   writeln ('# points','     pi     ','  (4/3)*pi  ',' (1/2)*pi^2 ');
   writeln;
   FOR j := 1 to 13 DO BEGIN
      FOR k := twotoj(j-1) to twotoj(j) DO BEGIN
         x1 := ran1(idum);
         x2 := ran1(idum);
         x3 := ran1(idum);
         x4 := ran1(idum);
         IF (fnc(x1,x2,0.0,0.0) < 1.0) THEN  iy[1] := iy[1]+1;
         IF (fnc(x1,x2,x3,0.0) < 1.0) THEN  iy[2] := iy[2]+1;
         IF (fnc(x1,x2,x3,x4) < 1.0) THEN  iy[3] := iy[3]+1
      END;
      jpower := twotoj(j);
      yprob[1] := 4.0*iy[1]/jpower;
      yprob[2] := 8.0*iy[2]/jpower;
      yprob[3] := 16.0*iy[3]/jpower;
      writeln (jpower:6,yprob[1]:12:6,yprob[2]:12:6,yprob[3]:12:6)
   END;
   writeln;
   writeln ('actual',pi:12:6,(4.0*pi/3.0):12:6,(0.5*sqr(pi)):12:6)
END.

PROGRAM d7r3 (input,output);
(* driver for routine RAN2 *)
(* calculates pi statistically using volume of unit n-sphere *)
CONST
   pi=3.1415926;
VAR
   i,idum,j,k,jpower : integer;
   x1,x2,x3,x4 : real;
   iy : ARRAY [1..3] OF integer;
   yprob : ARRAY [1..3] OF real;
   gliy : integer;
   glir : ARRAY [1..97] OF integer;

FUNCTION fnc(x1,x2,x3,x4 : real) : real;
BEGIN
   fnc := sqrt(sqr(x1)+sqr(x2)+sqr(x3)+sqr(x4))
END;

FUNCTION twotoj(j : integer) : integer;
BEGIN
   IF (j=0) THEN twotoj := 1
   ELSE twotoj := 2*twotoj(j-1)
END;

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

BEGIN
   idum := -1;
   FOR i := 1 to 3 DO BEGIN
      iy[i] := 0
   END;
   writeln;
   writeln ('volume of unit n-sphere, n := 2,3,4');
   writeln ('# points','     pi     ','  (4/3)*pi  ',' (1/2)*pi^2 ');
   writeln;
   FOR j := 1 to 13 DO BEGIN
      FOR k := twotoj(j-1) to twotoj(j) DO BEGIN
         x1 := ran2(idum);
         x2 := ran2(idum);
         x3 := ran2(idum);
         x4 := ran2(idum);
         IF (fnc(x1,x2,0.0,0.0) < 1.0) THEN  iy[1] := iy[1]+1;
         IF (fnc(x1,x2,x3,0.0) < 1.0) THEN  iy[2] := iy[2]+1;
         IF (fnc(x1,x2,x3,x4) < 1.0) THEN  iy[3] := iy[3]+1
      END;
      jpower := twotoj(j);
      yprob[1] := 4.0*iy[1]/jpower;
      yprob[2] := 8.0*iy[2]/jpower;
      yprob[3] := 16.0*iy[3]/jpower;
      writeln (jpower:6,yprob[1]:12:6,yprob[2]:12:6,yprob[3]:12:6)
   END;
   writeln;
   writeln ('actual',pi:12:6,(4.0*pi/3.0):12:6,(0.5*sqr(pi)):12:6)
END.

PROGRAM d7r4 (input,output);
(* driver for routine RAN3 *)
(* calculates pi statistically using volume of unit n-sphere *)
CONST
   pi=3.1415926;
VAR
   i,idum,j,k,jpower : integer;
   x1,x2,x3,x4 : real;
   iy : ARRAY [1..3] OF integer;
   yprob : ARRAY [1..3] OF real;
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;

FUNCTION fnc(x1,x2,x3,x4 : real) : real;
BEGIN
   fnc := sqrt(sqr(x1)+sqr(x2)+sqr(x3)+sqr(x4))
END;

FUNCTION twotoj(j : integer) : integer;
BEGIN
   IF (j=0) THEN twotoj := 1
   ELSE twotoj := 2*twotoj(j-1)
END;

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

BEGIN
   idum := -1;
   FOR i := 1 to 3 DO BEGIN
      iy[i] := 0
   END;
   writeln;
   writeln ('volume of unit n-sphere, n := 2,3,4');
   writeln ('# points','     pi     ','  (4/3)*pi  ',' (1/2)*pi^2 ');
   writeln;
   FOR j := 1 to 13 DO BEGIN
      FOR k := twotoj(j-1) to twotoj(j) DO BEGIN
         x1 := ran3(idum);
         x2 := ran3(idum);
         x3 := ran3(idum);
         x4 := ran3(idum);
         IF (fnc(x1,x2,0.0,0.0) < 1.0) THEN  iy[1] := iy[1]+1;
         IF (fnc(x1,x2,x3,0.0) < 1.0) THEN  iy[2] := iy[2]+1;
         IF (fnc(x1,x2,x3,x4) < 1.0) THEN  iy[3] := iy[3]+1
      END;
      FOR i := 1 to 3 DO BEGIN
         jpower := twotoj(j);
         yprob[1] := 4.0*iy[1]/jpower;
         yprob[2] := 8.0*iy[2]/jpower;
         yprob[3] := 16.0*iy[3]/jpower
      END;
      writeln (jpower:6,yprob[1]:12:6,yprob[2]:12:6,yprob[3]:12:6)
   END;
   writeln;
   writeln ('actual',pi:12:6,(4.0*pi/3.0):12:6,(0.5*sqr(pi)):12:6)
END.

PROGRAM d7r12(input,output);
(* driver for routine RAN4 *)
CONST
   npt=50;
TYPE
   gl64array = ARRAY [1..64] OF integer;
   gl56array = ARRAY [1..56] OF integer;
   gl48array = ARRAY [1..48] OF integer;
   gl32array = ARRAY [1..32] OF integer;
   glnparray = ARRAY [1..npt] OF real;
   gl65reals = ARRAY [1..65] OF real;
VAR
   glnewkey : integer;
   glinp,glkey,ip,ipm : gl64array;
   glpow : gl65reals;
   glicd,ipc1 : gl56array;
   ave,vrnce : real;
   idum,j : integer;
   y : glnparray;
   ipc2,iet : gl48array;
   ipp : gl32array;
   is : ARRAY [1..16,1..4,1..8] OF integer;
   ibin : ARRAY [1..4,1..16] OF integer;
   ksflg,cyflg,desflg : boolean;

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

(*$I DESKS.PAS *)

(*$I DES.PAS *)

(*$I RAN4.PAS *)

BEGIN
   ksflg := true;
   cyflg := true;
   desflg := true;
   idum := -123;
   ave := 0.0;
   writeln('First 10 Random Numbers with idum = ',idum:5);
   writeln;
   writeln('#':4,'RAN4':11);
   FOR j := 1 to 10 DO y[j] := ran4(idum);
   FOR j := 1 to 10 DO writeln(j:4,y[j]:12:6);
   writeln;
   writeln('Average and Variance of Next ',npt:3);
   FOR j := 1 to npt DO BEGIN
      y[j] := ran4(idum)
   END;
   avevar(y,npt,ave,vrnce);
   writeln;
   writeln('Average: ',ave:10:4);
   writeln('Variance:',vrnce:10:4);
   writeln;
   writeln('Expected Result for an Infinite Sample:');
   writeln;
   writeln('Average: ',0.5:10:4);
   writeln('Variance:',1.0/12.0:10:4)
END.
PROGRAM d8r7(input,output,dfile);
(* driver for routine RANK *)
CONST
   np=100;
TYPE
   gliarray = ARRAY [1..np] OF integer;
   glsarray = ARRAY [1..np] OF real;
VAR
   i,j,k,l : integer;
   a : glsarray;
   b : ARRAY [1..10] OF real;
   indx,irank : gliarray;
   dfile : text;

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

(*$I RANK.PAS *)

BEGIN
   glopen(dfile,'tarray.dat');
   readln(dfile);
   FOR i := 1 to 100 DO BEGIN
      read(dfile,a[i])
   END;
   close(dfile);
   indexx(np,a,indx);
   rank(np,indx,irank);
   writeln('original array is:');
   FOR i := 1 to 10 DO BEGIN
      FOR j := 1 to 10 DO BEGIN
         write (a[10*(i-1)+j]:6:2)
      END;
      writeln
   END;
   writeln('table of ranks is:');
   FOR i := 1 to 10 DO BEGIN
      FOR j := 1 to 10 DO BEGIN
         write(irank[10*(i-1)+j]:6)
      END;
      writeln
   END;
   writeln('press return to continue...');
   readln;
   writeln('array sorted according to rank table:');
   FOR i := 1 to 10 DO BEGIN
      FOR j := 1 to 10 DO BEGIN
         k := 10*(i-1)+j;
         FOR l := 1 to 100 DO BEGIN
            IF (irank[l] = k) THEN b[j] := a[l]
         END
      END;
      FOR j := 1 to 10 DO BEGIN
         write(b[j]:6:2)
      END;
   writeln
   END
END.
PROGRAM d3r2(input,output);
(* driver for routine RATINT *)
CONST
   npt=6;
   eps=1.0;
TYPE
   glnarray = ARRAY [1..npt] OF real;
VAR
   dyy,xx,yexp,yy : real;
   i : integer;
   x,y : glnarray;

FUNCTION f(x,eps: real): real;
BEGIN
   f := x*exp(-x)/(sqr(x-1.0)+eps*eps)
END;

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

BEGIN
   FOR i := 1 to npt DO BEGIN
      x[i] := i*2.0/npt;
      y[i] := f(x[i],eps)
   END;
   writeln('Diagonal rational function interpolation');
   writeln;
   writeln('x':5,'interp.':13,'accuracy':14,'actual':12);
   FOR i := 1 to 10 DO BEGIN
      xx := 0.2*i;
      ratint(x,y,npt,xx,yy,dyy);
      yexp := f(xx,eps);
      writeln(xx:6:2,yy:12:6,' ':4,dyy:11,yexp:12:6)
   END
END.
PROGRAM d12r3(input,output);
(* driver for routine REALFT *)
LABEL 1,99;
CONST
   eps=1.0e-3;
   np=32;
   np2=34;   (* np+2 *)
   width=50.0;
   pi=3.1415926;
TYPE
   gldarray = ARRAY [1..np2] OF real;
VAR
   big,per,scal,small : real;
   i,j,n,nlim : integer;
   data,size : gldarray;

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

(*$I REALFT.PAS *)

BEGIN
   n := np DIV 2;
1:   writeln('Period of sinusoid in channels (2-',np:2,')');
   readln(per);
   IF (per <= 0.0) THEN GOTO 99;
   FOR i := 1 to np DO BEGIN
      data[i] := cos(2.0*pi*(i-1)/per)
   END;
   realft(data,n,+1);
   big := -1.0e10;
   FOR i := 1 to n DO BEGIN
      size[i] := sqrt(sqr(data[2*i-1])+sqr(data[2*i]));
      IF (size[i] > big) THEN big := size[i]
   END;
   size[1] := data[1];
   IF (size[1] > big) THEN big := size[1];
   scal := width/big;
   FOR i := 1 to n DO BEGIN
      nlim := round(scal*size[i]+eps);
      write(i:4,' ');
      FOR j := 1 to nlim+1 DO write('*');
      writeln
   END;
   writeln('press RETURN to continue ...');
   readln;
   realft(data,n,-1);
   big := -1.0e10;
   small := 1.0e10;
   FOR i := 1 to np DO BEGIN
      IF (data[i] < small) THEN small := data[i];
      IF (data[i] > big) THEN big := data[i]
   END;
   scal := width/(big-small);
   FOR i := 1 to np DO BEGIN
      nlim := round(scal*(data[i]-small)+eps);
      write(i:4,' ');
      FOR j := 1 to nlim+1 DO write('*');
      writeln
   END;
   GOTO 1;
99:
END.
PROGRAM d15r1(input,output);
(* driver for routine RK4 *)
CONST
   n=4;
TYPE
   glnarray = ARRAY [1..n] OF real;
VAR
   h,x : real;
   i,j : integer;
   y,dydx,yout : glnarray;

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

(*$I BESSJ1.PAS *)

(*$I BESSJ.PAS *)

PROCEDURE derivs(x : real; y : glnarray;
       VAR dydx : glnarray);
(* The arrays y and dydx must carry the dimension given to
them in the calling 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 RK4.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];
   writeln;
   writeln('Bessel function:':16,'j0':5,
         'j1':12,'j3':12,'j4':12);
   FOR i := 1 to 5 DO BEGIN
      h := 0.2*i;
      rk4(y,dydx,n,x,h,yout);
      writeln;
      writeln('for a step size of:',h:6:2);
      write('rk4: ':11);
      FOR j := 1 to 4 DO write(yout[j]:12:6);
      writeln;
      writeln('actual: ':11,bessj0(x+h):12:6,
         bessj1(x+h):12:6,bessj(2,x+h):12:6,bessj(3,x+h):12:6)
   END
END.
PROGRAM d15r2(input,output);
(* driver for routine RKDUMB *)
CONST
   nvar=4;
   nstep=150;
TYPE
   glnarray = ARRAY [1..nvar] OF real;
VAR
   i,j : integer;
   x1,x2 : real;
   vstart : glnarray;
   xx : ARRAY [1..200] OF real;
   y : ARRAY [1..nvar,1..200] OF real;

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

(*$I BESSJ1.PAS *)

(*$I BESSJ.PAS *)

PROCEDURE derivs(x : real; y : glnarray;
       VAR dydx : glnarray);
(* The arrays y and dydx must carry the dimension given to
them in the calling 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 RK4.PAS *)

(*$I RKDUMB.PAS *)

BEGIN
   x1 := 1.0;
   vstart[1] := bessj0(x1);
   vstart[2] := bessj1(x1);
   vstart[3] := bessj(2,x1);
   vstart[4] := bessj(3,x1);
   x2 := 20.0;
   rkdumb(vstart,nvar,x1,x2,nstep);
   writeln('x':8,'integrated':17,'bessj3':10);
   FOR i := 1 to (nstep DIV 10) DO BEGIN
      j := 10*i;
      writeln(xx[j]:10:4,'  ',y[4,j]:12:6,bessj(3,xx[j]):12:6)
   END
END.
PROGRAM d15r3(input,output);
(* driver for routine RKQC *)
CONST
   n=4;
TYPE
   glarray = ARRAY [1..n] OF real;
   glnarray = glarray;
VAR
   eps,hdid,hnext,htry,x : real;
   i : integer;
   y,dydx,yscal : glarray;

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

(*$I BESSJ1.PAS *)

(*$I BESSJ.PAS *)

PROCEDURE derivs(x: real; y: glarray; VAR dydx: glarray);
(* Programs using derivs must define type
TYPE
   glarray = ARRAY [1..4] OF real; *)
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 RK4.PAS *)

(*$I RKQC.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 := 0.1;
   writeln('eps':8,'htry':13,'hdid':12,'hnext':13);
   FOR i := 1 to 15 DO BEGIN
      eps := exp(-i);
      rkqc(y,dydx,n,x,htry,eps,yscal,hdid,hnext);
      writeln(eps:13,htry:8:2,hdid:14:6,hnext:12:6)
   END
END.
PROGRAM d14r12(input,output);
(* driver for routine ROFUNC *)
CONST
   spread=0.05;
   ndata=100;
TYPE
   glsarray = ARRAY [1..ndata] OF real;
VAR
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;
   gliset : integer;
   glgset : real;
   aa,abdevt,b,rf : real;
   i,idum,npt : integer;
   arr,x,y : glsarray;

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

(*$I GASDEV.PAS *)

(*$I SORT.PAS *)

(*$I ROFUNC.PAS *)

BEGIN
   gliset := 0;
   idum := -11;
   FOR i := 1 to ndata DO BEGIN
      x[i] := 0.1*i;
      y[i] := -2.0*x[i]+1.0+spread*gasdev(idum)
   END;
   writeln('b':9,'a':9,'ROFUNC':12,'ABDEVT':10);
   FOR i := -5 to 5 DO BEGIN
      b := -2.0+0.02*i;
      rf := rofunc(b);
      writeln(b:10:2,aa:9:2,rofunc(b):11:2,abdevt:10:2)
   END
END.
PROGRAM d9r4(input,output);
(* driver for routine RTBIS *)
CONST
   n=100;
   nbmax=20;
   x1=1.0;
   x2=50.0;
TYPE
   glnbmax = ARRAY [1..nbmax] OF real;
VAR
   i,nb : integer;
   xacc,root : real;
   xb1,xb2 : glnbmax;

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

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

(*$I ZBRAK.PAS *)

(*$I RTBIS.PAS *)

BEGIN
   nb := nbmax;
   zbrak(x1,x2,n,xb1,xb2,nb);
   writeln;
   writeln('roots of bessj0:');
   writeln('x':20,'f(x)':15);
   FOR i := 1 to nb DO BEGIN
      xacc := (1.0e-6)*(xb1[i]+xb2[i])/2.0;
      root := rtbis(xb1[i],xb2[i],xacc);
      writeln('  root ',i:2,'  ',root:12:6,fx(root):14:6)
   END
END.
PROGRAM d9r5(input,output);
(* driver for routine RTFLSP *)
CONST
   n=100;
   nbmax=20;
   x1=1.0;
   x2=50.0;
TYPE
   glnbmax = ARRAY [1..nbmax] OF real;
VAR
   i,nb : integer;
   root,xacc : real;
   xb1,xb2 : glnbmax;

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

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

(*$I ZBRAK.PAS *)

(*$I RTFLSP.PAS *)

BEGIN
   nb := nbmax;
   zbrak(x1,x2,n,xb1,xb2,nb);
   writeln;
   writeln('roots of bessj0:');
   writeln('x':20,'f(x)':15);
   FOR i := 1 to nb DO BEGIN
      xacc := (1.0e-6)*(xb1[i]+xb2[i])/2.0;
      root := rtflsp(xb1[i],xb2[i],xacc);
      writeln('  root ',i:2,'  ',root:12:6,fx(root):14:6)
   END
END.
PROGRAM d9r8(input,output);
(* driver for routine RTNEWT *)
CONST
   n=100;
   nbmax=20;
   x1=1.0;
   x2=50.0;
TYPE
   glnbmax = ARRAY [1..nbmax] OF real;
VAR
   i,nb : integer;
   root,xacc : real;
   xb1,xb2 : glnbmax;

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

(*$I BESSJ1.PAS *)

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

PROCEDURE funcd(x: real; VAR fn,df: real);
BEGIN
   fn := fx(x);
   df := -bessj1(x)
END;

(*$I ZBRAK.PAS *)

(*$I RTNEWT.PAS *)

BEGIN
   nb := nbmax;
   zbrak(x1,x2,n,xb1,xb2,nb);
   writeln;
   writeln('roots of bessj0:');
   writeln('x':20,'f(x)':15);
   FOR i := 1 to nb DO BEGIN
      xacc := (1.0e-6)*(xb1[i]+xb2[i])/2.0;
      root := rtnewt(xb1[i],xb2[i],xacc);
      writeln('  root ',i:2,'  ',root:12:6,fx(root):14:6)
   END
END.
PROGRAM d9r9(input,output);
(* driver for routine RTSAFE *)
CONST
   n=100;
   nbmax=20;
   x1=1.0;
   x2=50.0;
TYPE
   glnbmax = ARRAY [1..nbmax] OF real;
VAR
   i,nb : integer;
   root,xacc : real;
   xb1,xb2 : glnbmax;

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

(*$I BESSJ1.PAS *)

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

PROCEDURE funcd(x: real; VAR fn,df: real);
BEGIN
   fn := fx(x);
   df := -bessj1(x)
END;

(*$I ZBRAK.PAS *)

(*$I RTSAFE.PAS *)

BEGIN
   nb := nbmax;
   zbrak(x1,x2,n,xb1,xb2,nb);
   writeln;
   writeln('roots of bessj0:');
   writeln('x':19,'f(x)':16);
   FOR i := 1 to nb DO BEGIN
      xacc := (1.0e-6)*(xb1[i]+xb2[i])/2.0;
      root := rtsafe(xb1[i],xb2[i],xacc);
      writeln(' root ',i:2,'  ',root:12:6,fx(root):14:6)
   END
END.
PROGRAM d9r6(input,output);
(* driver for routine RTSEC *)
CONST
   n=100;
   nbmax=20;
   x1=1.0;
   x2=50.0;
TYPE
   glnbmax = ARRAY [1..nbmax] OF real;
VAR
   i,nb : integer;
   root,xacc : real;
   xb1,xb2 : glnbmax;

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

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

(*$I ZBRAK.PAS *)

(*$I RTSEC.PAS *)

BEGIN
   nb := nbmax;
   zbrak(x1,x2,n,xb1,xb2,nb);
   writeln;
   writeln('roots of bessj0:');
   writeln('x':20,'f(x)':15);
   FOR i := 1 to nb DO BEGIN
      xacc := (1.0e-6)*(xb1[i]+xb2[i])/2.0;
      root := rtsec(xb1[i],xb2[i],xacc);
      writeln('  root ',i:2,'  ',root:12:6,fx(root):14:6)
   END
END.
PROGRAM d15r7(input,output);
(* driver for routine RZEXTR *)
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;
   gld : ARRAY [1..glnmax,1..glncol] OF real;
   dum,xest : real;
   i,iest,j : integer;
   dy,yest,yz : glyarray;

(*$I MODFILE.PAS *)
(*$I RZEXTR.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;
      rzextr(iest,xest,yest,yz,dy,nv,nuse);
      writeln;
      writeln('iest  :=  ',i:2,'   xest =',xest:8:4);
      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.
