FUNCTION ran0(VAR idum: integer): real;
(* RAN0 assumes that the system random number generator is called RANDOM,
with no arguments. Programs using RAN0 must declare the following
variables
VAR
   gly: real;
   glv: ARRAY [1..97] OF real;
in the main routine. *)
VAR
   dum: real;
   j: integer;
BEGIN
   IF (idum < 0) THEN BEGIN
      idum := 1;
      FOR j := 1 to 97 DO dum := random;
      FOR j := 1 to 97 DO glv[j] := random;
      gly := random
   END;
   j := 1+trunc(97.0*gly);
   IF ((j > 97) OR (j < 1)) THEN BEGIN
      writeln('pause in routine RAN0'); readln END;
   gly := glv[j];
   ran0 := gly;
   glv[j] := random
END;
FUNCTION ran1(VAR idum: integer): real;
(* Programs using RAN1 must declare the following variables
VAR
   glix1,glix2,glix3: integer;
   glr: ARRAY [1..97] OF real;
in the main program. *)
CONST
   m1=259200;
   ia1=7141;
   ic1=54773;
   rm1=3.8580247e-6;   (* 1.0/m1 *)
   m2=134456;
   ia2=8121;
   ic2=28411;
   rm2=7.4373773e-6;   (* 1.0/m2 *)
   m3=243000;
   ia3=4561;
   ic3=51349;
VAR
   j: integer;
BEGIN
   IF  (idum < 0)  THEN BEGIN
      glix1 := (ic1-idum) MOD m1;
      glix1 := (ia1*glix1+ic1) MOD m1;
      glix2 := glix1 MOD m2;
      glix1 := (ia1*glix1+ic1) MOD m1;
      glix3 := glix1 MOD m3;
      FOR j := 1 to 97 DO BEGIN
         glix1 := (ia1*glix1+ic1) MOD m1;
         glix2 := (ia2*glix2+ic2) MOD m2;
         glr[j] := (glix1+glix2*rm2)*rm1
      END;
      idum := 1
   END;
   glix1 := (ia1*glix1+ic1) MOD m1;
   glix2 := (ia2*glix2+ic2) MOD m2;
   glix3 := (ia3*glix3+ic3) MOD m3;
   j := 1 + (97*glix3) DIV m3;
   IF ((j > 97) OR (j < 1)) THEN BEGIN
      writeln('pause in routine RAN1'); readln
   END;
   ran1 := glr[j];
   glr[j] := (glix1+glix2*rm2)*rm1
END;
FUNCTION ran2(VAR idum: integer): real;
(* Programs using RAN2 must declare the following variables
VAR
   gliy: integer;
   glir: ARRAY [1..97] OF integer;
in the main program. *)
CONST
   m=714025;
   ia=1366;
   ic=150889;
   rm=1.400512e-6;      (* 1.0/m *)
VAR
   j: integer;
BEGIN
   IF (idum < 0) THEN BEGIN
      idum := (ic-idum) MOD m;
      FOR j := 1 to 97 DO BEGIN
         idum := (ia*idum+ic) MOD m;
         glir[j] := idum
      END;
      idum := (ia*idum+ic) MOD m;
      gliy := idum
   END;
   j := 1 + (97*gliy) DIV m;
   IF ((j > 97) OR (j < 1)) THEN BEGIN
      writeln('pause in routine RAN2'); readln
   END;
   gliy := glir[j];
   ran2 := gliy*rm;
   idum := (ia*idum+ic) MOD m;
   glir[j] := idum
END;
FUNCTION ran3(VAR idum: integer): real;
(* Programs using RAN3 must declare the following variables
VAR
   glinext,glinextp: integer;
   glma: ARRAY [1..55] OF real;
in the main routine. Machines with 4-byte integers can use the integer
implementation of this routine, substituting glma of type integer, the
commented CONST and VAR declarations, and the MOD function in the third
line after the BEGIN. *)
(* CONST
   mbig=1000000000;
   mseed=161803398;
   mz=0;
   fac=1.0e-9;
VAR
   i,ii,k,mj,mk: integer; *)
CONST
   mbig=4.0e6;
   mseed=1618033.0;
   mz=0.0;
   fac=2.5e-7; (* 1/mbig *)
VAR
   i,ii,k: integer;
   mj,mk: real;
BEGIN
   IF (idum < 0) THEN BEGIN
      mj := mseed+idum;
      (* The following IF block is mj := mj MOD mbig; for real variables. *)
      IF mj>=0.0 THEN mj := mj-mbig*trunc(mj/mbig)
         ELSE mj := mbig-abs(mj)+mbig*trunc(abs(mj)/mbig);
      glma[55] := mj;
      mk := 1;
      FOR i := 1 to 54 DO BEGIN
         ii := 21*i MOD 55;
         glma[ii] := mk;
         mk := mj-mk;
         IF (mk < mz) THEN mk := mk+mbig;
         mj := glma[ii]
      END;
      FOR k := 1 to 4 DO BEGIN
         FOR i := 1 to 55 DO BEGIN
            glma[i] := glma[i]-glma[1+((i+30) MOD 55)];
            IF (glma[i] < mz) THEN glma[i] := glma[i]+mbig
         END
      END;
      glinext := 0;
      glinextp := 31;
      idum := 1
   END;
   glinext := glinext+1;
   IF (glinext = 56) THEN glinext := 1;
   glinextp := glinextp+1;
   IF (glinextp = 56) THEN glinextp := 1;
   mj := glma[glinext]-glma[glinextp];
   IF (mj < mz) THEN mj := mj+mbig;
   glma[glinext] := mj;
   ran3 := mj*fac
END;
FUNCTION ran4(VAR idum: integer): real;
(* Programs using routine RAN4 must define the global variables
TYPE
   gl64array = ARRAY [1..64] OF integer;
   gl65reals = ARRAY [1..65] OF real;
VAR
   glnewkey: integer;
   glinp,glkey: gl64array;
   glpow: gl65reals;
in the main routine. The initialization block, IF (idum < 0), has been
written in real arithmetic to avoid overflow on machines with 2-byte
integers. With 4-byte integers, this block can be simplified with MOD
and DIV. *)
CONST
   im=11979;
   rm=11979.0;
   a=430.0;
   c=2531.0;
   nacc=24;
VAR
   isav,j: integer;
   jot: gl64array;
   r4,dum: real;
BEGIN
   IF (idum < 0) THEN BEGIN
      dum := idum MOD im;
      IF (dum < 0.0) THEN dum := dum+rm;
      glpow[1] := 0.5;
      FOR j := 1 to 64 DO BEGIN
         dum := dum*a+c;
         dum := dum-rm*trunc(dum/rm);
         glkey[j] := trunc(2.0*dum/rm);
         glinp[j] := trunc(4.0*dum/rm) MOD 2;
         glpow[j+1] := 0.5*glpow[j]
      END;
      idum := round(dum);
      glnewkey := 1
   END;
   isav := glinp[64];
   IF (isav <> 0) THEN BEGIN
      glinp[4] := 1-glinp[4];
      glinp[3] := 1-glinp[3];
      glinp[1] := 1-glinp[1]
   END;
   FOR j := 64 DOWNTO 2 DO BEGIN
      glinp[j] := glinp[j-1]
   END;
   glinp[1] := isav;
   des(glinp,glkey,glnewkey,0,jot);
   r4 := 0;
   FOR j := 1 to nacc DO BEGIN
      IF (jot[j] <> 0) THEN r4 := r4+glpow[j]
   END;
   ran4 := r4
END;
PROCEDURE rank(n: integer; indx: gliarray; VAR irank: gliarray);
(* Programs using routine RANK must define type
TYPE
   gliarray = ARRAY [1..np] OF real;
in the main routine, with np >= n.   *)
VAR
   j: integer;
BEGIN
   FOR j := 1 to n DO BEGIN
      irank[indx[j]] := j
   END
END;
PROCEDURE ratint(xa,ya: glnarray; n: integer; x: real; VAR y,dy: real);
(* Programs using routine RATINT must define the type
TYPE
   glnarray = ARRAY [1..n] OF real;
in the main routine. *)
LABEL 99;
CONST
   tiny=1.0e-25;
VAR
   ns,m,i: integer;
   w,t,hh,h,dd: real;
   c,d: glnarray;
BEGIN
   ns := 1;
   hh := abs(x-xa[1]);
   FOR i := 1 to n DO BEGIN
      h := abs(x-xa[i]);
      IF  (h = 0.0) THEN BEGIN
         y := ya[i];
         dy := 0.0;
         GOTO 99 END
      ELSE IF (h < hh) THEN BEGIN
         ns := i;
         hh := h
      END;
      c[i] := ya[i];
      d[i] := ya[i]+tiny
   END;
   y := ya[ns];
   ns := ns-1;
   FOR m := 1 to n-1 DO BEGIN
      FOR i := 1 to n-m DO BEGIN
         w := c[i+1]-d[i];
         h := xa[i+m]-x;
         t := (xa[i]-x)*d[i]/h;
         dd := t-c[i+1];
         IF (dd = 0.0) THEN BEGIN
            writeln('pause in routine RATINT'); readln
         END;
         dd := w/dd;
         d[i] := c[i+1]*dd;
         c[i] := t*dd
      END;
      IF  (2*ns < n-m) THEN BEGIN
         dy := c[ns+1]
      END ELSE BEGIN
         dy := d[ns];
         ns := ns-1
      END;
      y := y+dy
   END;
99:   END;
PROCEDURE realft(VAR data: gldarray; n,isign: integer);
(* Programs using routine REALFT must define the type
TYPE
   gldarray = ARRAY [1..2*n+2] OF real;
where 2*n is the dimension of the input data array. When
routine FOUR1 is used with REALFT, its data type 'gldarray'
should be set as in this program. *)
VAR
   wr,wi,wpr,wpi,wtemp,theta: double;
   i,i1,i2,i3,i4: integer;
   c1,c2,h1r,h1i,h2r,h2i,wrs,wis: real;
BEGIN
   theta := 6.28318530717959/(2.0*n);
   wr := 1.0;
   wi := 0.0;
   c1 := 0.5;
   IF (isign = 1) THEN BEGIN
      c2 := -0.5;
      four1(data,n,1);
      data[2*n+1] := data[1];
      data[2*n+2] := data[2]
   END ELSE BEGIN
      c2 := 0.5;
      theta := -theta;
      data[2*n+1] := data[2];
      data[2*n+2] := 0.0;
      data[2] := 0.0
   END;
   wpr := -2.0*sqr(sin(0.5*theta));
   wpi := sin(theta);
   FOR i := 1 to (n DIV 2)+1 DO BEGIN
      i1 := i+i-1;
      i2 := i1+1;
      i3 := n+n+3-i2;
      i4 := i3+1;
      wrs := sngl(wr);
      wis := sngl(wi);
      h1r := c1*(data[i1]+data[i3]);
      h1i := c1*(data[i2]-data[i4]);
      h2r := -c2*(data[i2]+data[i4]);
      h2i := c2*(data[i1]-data[i3]);
      data[i1] := h1r+wrs*h2r-wis*h2i;
      data[i2] := h1i+wrs*h2i+wis*h2r;
      data[i3] := h1r-wrs*h2r+wis*h2i;
      data[i4] := -h1i+wrs*h2i+wis*h2r;
      wtemp := wr;
      wr := wr*wpr-wi*wpi+wr;
      wi := wi*wpr+wtemp*wpi+wi
   END;
   IF (isign = 1) THEN BEGIN
      data[2] := data[2*n+1]
   END ELSE BEGIN
      four1(data,n,-1)
   END
END;
PROCEDURE red(iz1,iz2,jz1,jz2,jm1,jm2,jmf,ic1,jc1,jcf,kc: integer;
         c: glcarray; nci,ncj,nck: integer;
         VAR s: glsarray; nsi,nsj: integer);
(* Programs using routine RED must define the types
TYPE
   glcarray = ARRAY [1..nci,1..ncj,1..nck] OF real;
   glsarray = ARRAY [1..nsi,1..nsj] OF real;
in the main routine. *)
VAR
   loff,l,j,ic,i: integer;
   vx: real;
BEGIN
   loff := jc1-jm1;
   ic := ic1;
   FOR j := jz1 to jz2 DO BEGIN
      FOR l := jm1 to jm2 DO BEGIN
         vx := c[ic,l+loff,kc];
         FOR i := iz1 to iz2 DO s[i,l] := s[i,l]-s[i,j]*vx
      END;
      vx := c[ic,jcf,kc];
      FOR i := iz1 to iz2 DO s[i,jmf] := s[i,jmf]-s[i,j]*vx;
      ic := ic+1
   END
END;
PROCEDURE rk4(y,dydx: glnarray; n: integer; x,h: real; VAR yout: glnarray);
(* Programs using routine RK4 must provide a
PROCEDURE derivs(x:real; y:glnarray; VAR dydx:glnarray);
which returns the derivatives dydx at location x, given both x and the
function values y. The calling program must also define the types
TYPE
   glnarray = ARRAY [1..nvar] OF real;
where nvar is the number of variables y. *)
VAR
   i: integer;
   xh,hh,h6: real;
   dym,dyt,yt: glnarray;
BEGIN
   hh := h*0.5;
   h6 := h/6.0;
   xh := x+hh;
   FOR i := 1 to n DO BEGIN
      yt[i] := y[i]+hh*dydx[i]
   END;
   derivs(xh,yt,dyt);
   FOR i := 1 to n DO BEGIN
      yt[i] := y[i]+hh*dyt[i]
   END;
   derivs(xh,yt,dym);
   FOR i := 1 to n DO BEGIN
      yt[i] := y[i]+h*dym[i];
      dym[i] := dyt[i]+dym[i]
   END;
   derivs(x+h,yt,dyt);
   FOR i := 1 to n DO BEGIN
      yout[i] := y[i]+h6*(dydx[i]+dyt[i]+2.0*dym[i])
   END
END;
PROCEDURE rkdumb(vstart: glnarray; nvar: integer; x1,x2: real; nstep: integer);
(* Programs using routine RKDUMB must provide a
PROCEDURE derivs(x:real; v:glnarray; VAR dvdx:glnarray);
which returns the derivatives dvdx at location x, given both x and the
function values v. They must also define the type
TYPE
   glnarray = ARRAY [1..nvar] OF real;
where nvar is the number of functions. They must also
declare the variables
VAR
   xx: ARRAY [1..200] OF real;
   y: ARRAY [1..nvar,1..200] OF real;
in the main routine. *)
VAR
   kount,k,i: integer;
   x,h: real;
   v,vout,dv: glnarray;
BEGIN
   FOR i := 1 to nvar DO BEGIN
      v[i] := vstart[i];
      y[i,1] := v[i]
   END;
   xx[1] := x1;
   x := x1;
   h := (x2-x1)/nstep;
   FOR k := 1 to nstep DO BEGIN
      derivs(x,v,dv);
      rk4(v,dv,nvar,x,h,vout);
      IF (x+h = x) THEN BEGIN
         writeln('pause in routine RKDUMB');
         writeln('stepsize to small'); readln
      END;
      x := x+h;
      xx[k+1] := x;
      FOR i := 1 to nvar DO BEGIN
         v[i] := vout[i];
         y[i,k+1] := v[i]
      END
   END
END;
PROCEDURE rkqc(VAR y,dydx: glarray; n: integer; VAR x: real;
      htry,eps: real; yscal: glarray; VAR hdid,hnext: real);
(* Programs using routine RKQC must provide a
PROCEDURE derivs(x:real; y:glnarray; VAR dydx:glnarray);
which returns the derivatives dydx at location x, given both x and the
function values y. They must also define the type
TYPE
   glarray = ARRAY [1..n] OF real;
in the main routine.   *)
LABEL 1;
CONST
   pgrow=-0.20;
   pshrnk=-0.25;
   fcor=0.06666666;   (* 1.0/15.0 *)
   one=1.0;
   safety=0.9;
   errcon=6.0e-4;
VAR
   i: integer;
   xsav,hh,h,temp,errmax: real;
   dysav,ysav,ytemp: glarray;
BEGIN
   xsav := x;
   FOR i := 1 to n DO BEGIN
      ysav[i] := y[i];
      dysav[i] := dydx[i]
   END;
   h := htry;
1:   hh := 0.5*h;
   rk4(ysav,dysav,n,xsav,hh,ytemp);
   x := xsav+hh;
   derivs(x,ytemp,dydx);
   rk4(ytemp,dydx,n,x,hh,y);
   x := xsav+h;
   IF (x = xsav) THEN BEGIN
      writeln('pause in routine RKQC');
      writeln('stepsize too small'); readln
   END;
   rk4(ysav,dysav,n,xsav,h,ytemp);
   errmax := 0.0;
   FOR i := 1 to n DO BEGIN
      ytemp[i] := y[i]-ytemp[i];
      temp := abs(ytemp[i]/yscal[i]);
      IF (errmax < temp) THEN errmax := temp
   END;
   errmax := errmax/eps;
   IF (errmax > one) THEN BEGIN
      h := safety*h*exp(pshrnk*ln(errmax));
      GOTO 1 END
   ELSE BEGIN
      hdid := h;
      IF (errmax > errcon) THEN BEGIN
         hnext := safety*h*exp(pgrow*ln(errmax))
      END ELSE BEGIN
         hnext := 4.0*h
      END
   END;
   FOR i := 1 to n DO BEGIN
      y[i] := y[i]+ytemp[i]*fcor
   END
END;
FUNCTION rofunc(b: real): real;
(* Programs using ROFUNC must declare
CONST
   ndata = [the number of data points]
VAR
   aa,abdevt: real;
   x,y,arr: ARRAY [1..ndata] OF real;
in the main routine. *)
VAR
   d,sum: real;
   j,n1,nmh,nml: integer;
BEGIN
   n1 := ndata+1;
   nml := n1 DIV 2;
   nmh := n1-nml;
   FOR j := 1 to ndata DO arr[j] := y[j]-b*x[j];
   sort(ndata,arr);
   aa := 0.5*(arr[nml]+arr[nmh]);
   sum := 0.0;
   abdevt := 0.0;
   FOR j := 1 to ndata DO BEGIN
      d := y[j]-(b*x[j]+aa);
      abdevt := abdevt+abs(d);
      IF (d > 0.0) THEN sum := sum+x[j] ELSE sum := sum-x[j]
   END;
   rofunc := sum
END;
FUNCTION rtbis(x1,x2,xacc: real): real;
(* Programs using routine RTBIS must externally define a function
fx(x:real):real which is to be analyzed for roots. *)
LABEL 99;
CONST
   jmax=40;
VAR
   dx,f,fmid,xmid,rtb: real;
   j: integer;
BEGIN
   fmid := fx(x2);
   f := fx(x1);
   IF ((f*fmid) >= 0.0) THEN BEGIN
      writeln('pause in RTBIS');
      writeln('Root must be bracketed for bisection.'); readln
   END;
   IF (f < 0.0) THEN BEGIN
      rtb := x1;
      dx := x2-x1
   END ELSE BEGIN
      rtb := x2;
      dx := x1-x2
   END;
   FOR j := 1 to jmax DO BEGIN
      dx := dx*0.5;
      xmid := rtb+dx;
      fmid := fx(xmid);
      IF (fmid < 0.0) THEN rtb := xmid;
      IF ((abs(dx) < xacc) OR (fmid = 0.0)) THEN GOTO 99
   END;
   writeln('pause in RTBIS - too many bisections'); readln;
99:   rtbis := rtb
END;
FUNCTION rtflsp(x1,x2,xacc: real): real;
(* Programs using routine RTFLSP must externally define a function
fx(x:real):real which is to be analyzed for roots. *)
LABEL 99;
CONST
   maxit=30;
VAR
   xl,xh,swap,fl: real;
   dx,del,f,fh,rtf: real;
   j: integer;
BEGIN
   fl := fx(x1);
   fh := fx(x2);
   IF (fl*fh > 0.0) THEN BEGIN
      writeln('pause in routine RTFLSP');
      writeln('Root must be bracketed for false position'); readln
   END;
   IF (fl < 0.0) THEN BEGIN
      xl := x1;
      xh := x2
   END ELSE BEGIN
      xl := x2;
      xh := x1;
      swap := fl;
      fl := fh;
      fh := swap
   END;
   dx := xh-xl;
   FOR j := 1 to maxit DO BEGIN
      rtf := xl+dx*fl/(fl-fh);
      f := fx(rtf);
      IF (f < 0.0)  THEN BEGIN
         del := xl-rtf;
         xl := rtf;
         fl := f
      END ELSE BEGIN
         del := xh-rtf;
         xh := rtf;
         fh := f
      END;
      dx := xh-xl;
      IF ((abs(del) < xacc) OR (f = 0.0)) THEN GOTO 99
   END;
   writeln('pause in routine RTFLSP');
   writeln('maximum number of iterations exceeded'); readln;
99:   rtflsp := rtf
END;
FUNCTION rtnewt(x1,x2,xacc: real): real;
(* Programs using routine RTNEWT must externally define procedure
funcd(x,f,df:real) which returns the function value f and its
derivative df at the point x. *)
LABEL 99;
CONST
   jmax=20;
VAR
   df,dx,f,rtn: real;
   j: integer;
BEGIN
   rtn := 0.5*(x1+x2);
   FOR j := 1 to jmax DO BEGIN
      funcd(rtn,f,df);
      dx := f/df;
      rtn := rtn-dx;
      IF ((x1-rtn)*(rtn-x2) < 0.0) THEN BEGIN
         writeln('pause in routine RTNEWT');
         writeln('jumped out of brackets'); readln
      END;
      IF (abs(dx) < xacc) THEN GOTO 99
   END;
   writeln('pause in routine RTNEWT');
   writeln('maximum number of iterations exceeded'); readln;
99:   rtnewt := rtn
END;
FUNCTION rtsafe(x1,x2,xacc: real): real;
(* Programs using routine RTSAFE must externally define procedure
funcd(x,f,df:real) which returns the function value f and its
derivative df at the point x. *)
LABEL 99;
CONST
   maxit=100;
VAR
   df,dx,dxold,f,fh,fl: real;
   swap,temp,xh,xl,rts: real;
   j: integer;
BEGIN
   funcd(x1,fl,df);
   funcd(x2,fh,df);
   IF (fl*fh >= 0.0) THEN BEGIN
      writeln('pause in routine RTSAFE');
      writeln('root must be bracketed'); readln
   END;
   IF (fl < 0.0) THEN BEGIN
      xl := x1;
      xh := x2
   END ELSE BEGIN
      xh := x1;
      xl := x2;
      swap := fl;
      fl := fh;
      fh := swap
   END;
   rts := 0.5*(x1+x2);
   dxold := abs(x2-x1);
   dx := dxold;
   funcd(rts,f,df);
   FOR j := 1 to maxit DO BEGIN
      IF((((rts-xh)*df-f)*((rts-xl)*df-f) >= 0.0)
      OR (abs(2.0*f) > abs(dxold*df))) THEN BEGIN
         dxold := dx;
         dx := 0.5*(xh-xl);
         rts := xl+dx;
         IF (xl = rts) THEN GOTO 99 END
      ELSE BEGIN
         dxold := dx;
         dx := f/df;
         temp := rts;
         rts := rts-dx;
         IF (temp = rts) THEN GOTO 99
      END;
      IF (abs(dx) < xacc) THEN GOTO 99;
      funcd(rts,f,df);
      IF (f < 0.0)  THEN BEGIN
         xl := rts;
         fl := f
      END ELSE BEGIN
         xh := rts;
         fh := f
      END
   END;
   writeln('pause in RTSAFE');
   writeln('maximum number of iterations exceeded'); readln;
99:   rtsafe := rts
END;
FUNCTION rtsec(x1,x2,xacc: real): real;
(* Programs using routine RTSEC must externally define a
function fx(x:real):real which is to be analyzed for roots. *)
LABEL 99;
CONST
   maxit=30;
VAR
   dx,f,fl,swap,xl,rts: real;
   j: integer;
BEGIN
   fl := fx(x1);
   f := fx(x2);
   IF (abs(fl) < abs(f)) THEN BEGIN
      rts := x1;
      xl := x2;
      swap := fl;
      fl := f;
      f := swap
   END ELSE BEGIN
      xl := x1;
      rts := x2
   END;
   FOR j := 1 to maxit DO BEGIN
      dx := (xl-rts)*f/(f-fl);
      xl := rts;
      fl := f;
      rts := rts+dx;
      f := fx(rts);
      IF ((abs(dx) < xacc) OR (f = 0.0)) THEN GOTO 99
   END;
   writeln('pause in routine RTSEC');
   writeln('maximum number of iterations exceeded'); readln;
99:   rtsec := rts
END;
PROCEDURE rzextr(iest: integer; xest: real; yest: glyarray;
       VAR yz,dy: glyarray; nv,nuse: integer);
(* Programs using routine RZEXTR must declare
TYPE
   glyarray = ARRAY [1..nv] OF real;
CONST
   glimax=11;
   glnmax=10;
   glncol=7;
VAR
   glx: ARRAY [1..glimax] OF real;
   gld: ARRAY [1..glnmax,1..glncol] OF real;
in the main routine. *)
CONST
   ncol=7;
VAR
   m1,k,j: integer;
   yy,v,ddy,c,b1,b: real;
   fx: ARRAY [1..ncol] OF real;
BEGIN
   glx[iest] := xest;
   IF (iest = 1) THEN BEGIN
      FOR j := 1 to nv DO BEGIN
         yz[j] := yest[j];
         gld[j,1] := yest[j];
         dy[j] := yest[j]
      END
   END ELSE BEGIN
      IF (iest < nuse) THEN m1 := iest ELSE m1 := nuse;
      FOR k := 1 to m1-1 DO BEGIN
         fx[k+1] := glx[iest-k]/xest
      END;
      FOR j := 1 to nv DO BEGIN
         yy := yest[j];
         v := gld[j,1];
         c := yy;
         gld[j,1] := yy;
         FOR k := 2 to m1 DO BEGIN
            b1 := fx[k]*v;
            b := b1-c;
            IF (b <> 0.0) THEN BEGIN
               b := (c-v)/b;
               ddy := c*b;
               c := b1*b
            END ELSE BEGIN
               ddy := v
            END;
            v := gld[j,k];
            gld[j,k] := ddy;
            yy := yy+ddy
         END;
         dy[j] := ddy;
         yz[j] := yy
      END
   END
END;
