PROCEDURE laguer(a: glcarray; m: integer; VAR x: gl2array;
       eps: real; polish: boolean);
(* Programs using routine LAGUER must define the type
TYPE
   glcarray = ARRAY [1..2*m+2] OF real;
   gl2array = ARRAY [1..2] OF real;
in the main routine. *)
LABEL 99;
CONST
   pi=3.1415928;
   tiny=1.0e-15;
   mxit=100;
VAR
   j,iter,npol: integer;
   xx,yy,rr,qq,dum,angl,eps2: real;
   dxold,cdx: real;
   sq,h,gp,gm,g2,g: gl2array;
   b,d,dx,f,x1: gl2array;
BEGIN
   eps2 := eps*eps;
   IF (polish) THEN BEGIN
      dxold := sqrt(sqr(x[1])+sqr(x[2]));
      npol := 0
   END;
   FOR iter := 1 to mxit DO BEGIN
      b[1] := a[2*m+1];
      b[2] := a[2*m+2];
      d[1] := 0.0;
      d[2] := 0.0;
      f[1] := 0.0;
      f[2] := 0.0;
      FOR j := m DOWNTO 1 DO BEGIN
         dum := f[1];
         f[1] := x[1]*f[1]-x[2]*f[2]+d[1];
         f[2] := x[1]*f[2]+x[2]*dum+d[2];
         dum := d[1];
         d[1] := x[1]*d[1]-x[2]*d[2]+b[1];
         d[2] := x[1]*d[2]+x[2]*dum+b[2];
         dum := b[1];
         b[1] := x[1]*b[1]-x[2]*b[2]+a[2*j-1];
         b[2] := x[1]*b[2]+x[2]*dum+a[2*j]
      END;
      IF (sqrt(sqr(b[1])+sqr(b[2])) < tiny) THEN BEGIN
         dx[1] := 0.0;
         dx[2] := 0.0
      END ELSE IF ((sqrt(sqr(d[1])+sqr(d[2])) < tiny) AND
         (sqrt(sqr(f[1])+sqr(f[2])) < tiny)) THEN BEGIN
         dx[1] := exp((1.0/m)*ln(sqrt(sqr(b[1])+sqr(b[2]))/
             sqrt(sqr(a[2*m+1])+sqr(a[2*m+2]))));
         dx[2] := 0.0
      END ELSE BEGIN
         g[1] := (b[1]*d[1]+d[2]*b[2])/(sqr(b[1])+sqr(b[2]));
         g[2] := (b[1]*d[2]-b[2]*d[1])/(sqr(b[1])+sqr(b[2]));
         g2[1] := sqr(g[1])-sqr(g[2]);
         g2[2] := 2.0*g[1]*g[2];
         h[1] := g2[1]-2.0*(b[1]*f[1]+b[2]*f[2])/
            (sqr(b[1])+sqr(b[2]));
         h[2] := g2[2]-2.0*(b[1]*f[2]-b[2]*f[1])/
            (sqr(b[1])+sqr(b[2]));
         xx := (m-1)*(m*h[1]-g2[1]);
         yy := (m-1)*(m*h[2]-g2[2]);
         IF ((abs(xx) < tiny) AND (abs(yy) < tiny)) THEN BEGIN
            sq[1] := 0.0;
            sq[2] := 0.0
         END ELSE BEGIN
            IF (abs(xx) > abs(yy)) THEN BEGIN
               rr := abs(xx)*sqrt(1.0+sqr(yy/xx))
            END ELSE BEGIN
               rr := abs(yy)*sqrt(1.0+sqr(xx/yy))
            END;
            qq := rr+abs(xx);
            IF (xx > 0.0) THEN BEGIN
               sq[1] := sqrt(0.5*qq);
               sq[2] := abs(yy)*sqrt(0.5/qq)
            END ELSE BEGIN
               sq[1] := abs(yy)*sqrt(0.5/qq);
               sq[2] := sqrt(0.5*qq)
            END;
            IF (yy < 0.0) THEN sq[2] := -sq[2]
         END;
         gp[1] := g[1]+sq[1];
         gp[2] := g[2]+sq[2];
         gm[1] := g[1]-sq[1];
         gm[2] := g[2]-sq[2];
         IF((sqr(gp[1])+sqr(gp[2])) < (sqr(gm[1])+sqr(gm[2])))
         THEN BEGIN
            gp[1] := gm[1];
            gp[2] := gm[2]
         END;
         dx[1] := m*gp[1]/(sqr(gp[1])+sqr(gp[2]));
         dx[2] := -m*gp[2]/(sqr(gp[1])+sqr(gp[2]))
      END;
      x1[1] := x[1]-dx[1];
      x1[2] := x[2]-dx[2];
      IF ((x[1] = x1[1]) AND (x[2] = x1[2])) THEN GOTO 99;
      x[1] := x1[1];
      x[2] := x1[2];
      IF (polish) THEN BEGIN
         npol := npol+1;
         cdx := sqrt(sqr(dx[1])+sqr(dx[2]));
         IF ((npol > 9) AND (cdx >= dxold)) THEN GOTO 99;
         dxold := cdx
      END ELSE IF ((sqr(dx[1])+sqr(dx[2]))
         <= (eps2*(sqr(x[1])+sqr(x[2])))) THEN BEGIN
         GOTO 99
      END
   END;
   writeln('pause in routine LAGUER - too many iterations'); readln;
99:   END;
PROCEDURE lfit(x,y,sig: glndata; ndata: integer; VAR a: glmma; mma: integer;
      lista: gllista; mfit: integer; VAR covar: glcovar;
      ncvm: integer; VAR chisq: real);
(* Programs using routine LFIT must define the types
TYPE
   glndata = ARRAY [1..ndata] OF real;
   glmma = ARRAY [1..mma] OF real;
   gllista = ARRAY [1..mma] OF integer;
   glcovar = ARRAY [1..ncvm,1..ncvm] OF real;
   glnpbymp = ARRAY [1..ncvm,1..1] OF real;
in the main routine. *)
VAR
   k,kk,j,ihit,i: integer;
   ym,wt,sum,sig2i: real;
   beta: glnpbymp;
   afunc: glmma;
BEGIN
   kk := mfit+1;
   FOR j := 1 to mma DO BEGIN
      ihit := 0;
      FOR k := 1 to mfit DO BEGIN
         IF (lista[k] = j) THEN ihit := ihit+1
      END;
      IF (ihit = 0) THEN BEGIN
         lista[kk] := j;
         kk := kk+1
      END ELSE IF (ihit > 1) THEN BEGIN
         writeln('pause in routine LFIT');
         writeln('improper permutation in LISTA'); readln
      END
   END;
   IF (kk <> (mma+1)) THEN BEGIN
      writeln('pause in routine LFIT');
      writeln('improper permutation in LISTA'); readln
   END;
   FOR j := 1 to mfit DO BEGIN
      FOR k := 1 to mfit DO BEGIN
         covar[j,k] := 0.0
      END;
      beta[j,1] := 0.0
   END;
   FOR i := 1 to ndata DO BEGIN
      funcs(x[i],afunc,mma);
      ym := y[i];
      IF (mfit < mma) THEN BEGIN
         FOR j := (mfit+1) to mma DO BEGIN
            ym := ym-a[lista[j]]*afunc[lista[j]]
         END
      END;
      sig2i := 1.0/sqr(sig[i]);
      FOR j := 1 to mfit DO BEGIN
         wt := afunc[lista[j]]*sig2i;
         FOR k := 1 to j DO BEGIN
            covar[j,k] := covar[j,k]+wt*afunc[lista[k]]
         END;
         beta[j,1] := beta[j,1]+ym*wt
      END
   END;
   IF (mfit > 1) THEN BEGIN
      FOR j := 2 to mfit DO BEGIN
         FOR k := 1 to j-1 DO BEGIN
            covar[k,j] := covar[j,k]
         END
      END
   END;
   gaussj(covar,mfit,ncvm,beta,1,1);
   FOR j := 1 to mfit DO BEGIN
      a[lista[j]] := beta[j,1]
   END;
   chisq := 0.0;
   FOR i := 1 to ndata DO BEGIN
      funcs(x[i],afunc,mma);
      sum := 0.0;
      FOR j := 1 to mma DO BEGIN
         sum := sum+a[j]*afunc[j]
      END;
      chisq := chisq+sqr((y[i]-sum)/sig[i])
   END;
   covsrt(covar,ncvm,mma,lista,mfit)
END;
PROCEDURE linmin(VAR p,xi: glnarray; n: integer; VAR fret: real);
(* Programs using routine LINMIN must define the type
TYPE
   glnarray = ARRAY [1..n] OF real;
They must also declare the variables
VAR
   ncom: integer;
   pcom,xicom: glnarray;
in the main routine. Also the function FUNC referenced by BRENT
and MNBRAK must be set to return the function F1DIM. *)
CONST
   tol=1.0e-4;
VAR
   j: integer;
   xx,xmin,fx,fb,fa,bx,ax: real;
BEGIN
   ncom := n;
   FOR j := 1 to n DO BEGIN
      pcom[j] := p[j];
      xicom[j] := xi[j]
   END;
   ax := 0.0;
   xx := 1.0;
   bx := 2.0;
   mnbrak(ax,xx,bx,fa,fx,fb);
   fret := brent(ax,xx,bx,tol,xmin);
   FOR j := 1 to n DO BEGIN
      xi[j] := xmin*xi[j];
      p[j] := p[j]+xi[j]
   END
END;
PROCEDURE locate(xx: glnarray; n: integer;
         x: real; VAR j: integer);
(* Programs which use routine LOCATE must define the type
TYPE
   glnarray = ARRAY [1..n] OF real;
in the main routine. *)
VAR
   ju,jm,jl: integer;
BEGIN
   jl := 0;
   ju := n+1;
   WHILE (ju-jl > 1) DO BEGIN
      jm := (ju+jl) DIV 2;
      IF ((xx[n] > xx[1]) = (x > xx[jm])) THEN jl := jm
      ELSE ju := jm
   END;
   j := jl
END;
PROCEDURE lubksb(a: glnpbynp; n,np: integer; indx: glindx; VAR b: glnarray);
(* Programs using LUBKSB must define the types
TYPE
   glnarray = ARRAY [1..n] OF real;
   glindx = ARRAY [1..n] OF integer;
   glnpbynp = ARRAY [1..np,1..np] OF real;
in the main routine *)
VAR
   j,ip,ii,i: integer;
   sum: real;
BEGIN
   ii := 0;
   FOR i := 1 to n DO BEGIN
      ip := indx[i];
      sum := b[ip];
      b[ip] := b[i];
      IF  (ii <> 0) THEN BEGIN
         FOR j := ii to i-1 DO BEGIN
            sum := sum-a[i,j]*b[j]
         END
      END ELSE IF (sum <> 0.0) THEN BEGIN
         ii := i
      END;
      b[i] := sum
   END;
   FOR i := n DOWNTO 1 DO BEGIN
      sum := b[i];
      IF (i < n) THEN BEGIN
         FOR j := i+1 to n DO BEGIN
            sum := sum-a[i,j]*b[j]
         END
      END;
      b[i] := sum/a[i,i]
   END
END;
PROCEDURE ludcmp(VAR a: glnpbynp; n,np: integer;
       VAR indx: glindx; VAR d: real);
(* Programs using LUDCMP must define the types
TYPE
   glnpbynp = ARRAY [1..np,1..np] OF real;
   glnarray = ARRAY [1..n] OF real;
   glindx = ARRAY [1..n] OF integer;
in the main routine. *)
CONST
   tiny=1.0e-20;
VAR
   k,j,imax,i: integer;
   sum,dum,big: real;
   vv: glnarray;
BEGIN
   d := 1.0;
   FOR i := 1 to n DO BEGIN
      big := 0.0;
      FOR j := 1 to n DO IF (abs(a[i,j]) > big) THEN big := abs(a[i,j]);
      IF (big = 0.0) THEN BEGIN
         writeln('pause in LUDCMP - singular matrix'); readln
      END;
      vv[i] := 1.0/big
   END;
   FOR j := 1 to n DO BEGIN
      IF (j > 1) THEN BEGIN
         FOR i := 1 to j-1 DO BEGIN
            sum := a[i,j];
            IF (i > 1) THEN BEGIN
               FOR k := 1 to i-1 DO BEGIN
                  sum := sum-a[i,k]*a[k,j]
               END;
               a[i,j] := sum
            END
         END
      END;
      big := 0.0;
      FOR i := j to n DO BEGIN
         sum := a[i,j];
         IF (j > 1) THEN BEGIN
            FOR k := 1 to j-1 DO BEGIN
               sum := sum-a[i,k]*a[k,j]
            END;
            a[i,j] := sum
         END;
         dum := vv[i]*abs(sum);
         IF (dum > big) THEN BEGIN
            big := dum;
            imax := i
         END
      END;
      IF (j <> imax) THEN BEGIN
         FOR k := 1 to n DO BEGIN
            dum := a[imax,k];
            a[imax,k] := a[j,k];
            a[j,k] := dum
         END;
         d := -d;
         vv[imax] := vv[j]
      END;
      indx[j] := imax;
      IF (j <> n) THEN BEGIN
         IF (a[j,j] = 0.0) THEN a[j,j] := tiny;
         dum := 1.0/a[j,j];
         FOR i := j+1 to n DO BEGIN
            a[i,j] := a[i,j]*dum
         END
      END
   END;
   IF (a[n,n] = 0.0) THEN a[n,n] := tiny
END;

