PROGRAM d9r10(input,output);
(* driver for routine LAGUER *)
CONST
   m=4;
   twomp2=10;   (* twomp2=2*m+2 *)
   ntry=21;
   twontry=42;   (* twontry=2*ntry *)
   eps=1.0e-6;
TYPE
   glcarray = ARRAY [1..twomp2] OF real;
   gl2array = ARRAY [1..2] OF real;
VAR
   i,iflag,j,n : integer;
   polish : boolean;
   x : gl2array;
   a : glcarray;
   y : ARRAY [1..twontry] OF real;

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

BEGIN
   a[1] := 0.0; a[2] := 2.0;
   a[3] := 0.0; a[4] := 0.0;
   a[5] := -1.0; a[6] := -2.0;
   a[7] := 0.0; a[8] := 0.0;
   a[9] := 1.0; a[10] := 0.0;
   writeln;
   writeln('Roots of polynomial x^4-(1+2i)*x^2+2i');
   writeln('real':15,'complex':13);
   n := 0;
   polish := false;
   FOR i := 1 to ntry DO BEGIN
      x[1] := (i-11.0)/10.0;
      x[2] := (i-11.0)/10.0;
      laguer(a,m,x,eps,polish);
      IF (n = 0) THEN BEGIN
         n := 1;
         y[1] := x[1];
         y[2] := x[2];
         writeln(n:5,x[1]:12:6,x[2]:12:6)
      END ELSE BEGIN
         iflag := 0;
         FOR j := 1 to n DO BEGIN
            IF (sqrt(sqr(x[1]-y[2*j-1])+sqr(x[2]-y[2*j]))
               <= eps*sqrt(sqr(x[1])+sqr(x[2]))) THEN
               iflag := 1            
         END;
         IF (iflag = 0) THEN BEGIN
            n := n+1;
            y[2*n-1] := x[1];
            y[2*n] := x[2];
            writeln(n:5,x[1]:12:6,x[2]:12:6)
         END
      END
   END
END.
PROGRAM d14r2(input,output);
(* driver for routine LFIT *)
CONST
   npt=100;
   spread=0.1;
   nterm=3;   
TYPE
   glcovar = ARRAY [1..nterm,1..nterm] OF real;
   glnpbynp = glcovar;
   glnpbymp = ARRAY [1..nterm,1..1] OF real;
   gllista = ARRAY [1..nterm] OF integer;
   glndata = ARRAY [1..npt] OF real;
   glmma = ARRAY [1..nterm] OF real;
VAR
   gliset : integer;
   glgset : real;
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;
   chisq : real;
   i,ii,idum,j,mfit : integer;
   lista : gllista;
   a : glmma;
   covar : glcovar;
   x,y,sig : glndata;

PROCEDURE funcs(x: real; VAR afunc: glmma; mma: integer);
(* Programs using FUNCS must define the type
TYPE
   glmma = ARRAY [1..mma] OF real;
in the main routine. *)
VAR
   i : integer;
BEGIN
   afunc[1] := 1.0;
   FOR i := 2 to mma DO BEGIN
      afunc[i] := x*afunc[i-1]
   END
END;

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

(*$I GASDEV.PAS *)

(*$I GAUSSJ.PAS *)

(*$I COVSRT.PAS *)

(*$I LFIT.PAS *)

BEGIN
   gliset := 0;
   idum := -911;
   FOR i := 1 to npt DO BEGIN
      x[i] := 0.1*i;
      y[i] := nterm;
      FOR j := nterm-1 DOWNTO 1 DO BEGIN
         y[i] := j+y[i]*x[i]
      END;
      y[i] := y[i]+spread*gasdev(idum);
      sig[i] := spread
   END;
   mfit := nterm;
   FOR i := 1 to mfit DO BEGIN
      lista[i] := i
   END;
   lfit(x,y,sig,npt,a,nterm,lista,mfit,covar,nterm,chisq);
   writeln;
   writeln('parameter':9,'uncertainty':23);
   FOR i := 1 to nterm DO BEGIN
      writeln('a[':4,i:1,'] = ',a[i]:8:6,sqrt(covar[i,i]):12:6)
   END;
   writeln('chi-squared = ',chisq:12);
   writeln('full covariance matrix');
   FOR i := 1 to nterm DO BEGIN
      FOR j := 1 to nterm DO write(covar[i,j]:12);
      writeln
   END;
   writeln;
   writeln('press RETURN to continue...');
   readln;
(* now test the LISTA feature *)
   FOR i := 1 to nterm DO BEGIN
      lista[i] := nterm+1-i
   END;
   lfit(x,y,sig,npt,a,nterm,lista,mfit,covar,nterm,chisq);
   writeln('parameter':9,'uncertainty':23);
   FOR i := 1 to nterm DO BEGIN
      writeln('a[':4,i:1,'] = ',a[i]:8:6,sqrt(covar[i,i]):12:6)
   END;
   writeln('chi-squared = ',chisq:12);
   writeln('full covariance matrix');
   FOR i := 1 to nterm DO BEGIN
      FOR j := 1 to nterm DO write(covar[i,j]:12);
      writeln
   END;
   writeln;
   writeln('press RETURN to continue...');
   readln;
(* now check results of restricting fit parameters *)
   ii := 1;
   FOR i := 1 to nterm DO BEGIN
      IF ((i MOD 2) = 1) THEN BEGIN
         lista[ii] := i;
         ii := ii+1
      END
   END;
   mfit := ii-1;
   lfit(x,y,sig,npt,a,nterm,lista,mfit,covar,nterm,chisq);
   writeln('parameter':9,'uncertainty':23);
   FOR i := 1 to nterm DO BEGIN
      writeln('a[':4,i:1,'] = ',a[i]:8:6,sqrt(covar[i,i]):12:6)
   END;
   writeln('chi-squared = ',chisq:12);
   writeln('full covariance matrix');
   FOR i := 1 to nterm DO BEGIN
      FOR j := 1 to nterm DO write(covar[i,j]:12);
      writeln
   END;
   writeln
END.
PROGRAM d10r7(input,output);
(* driver for routine LINMIN *)
CONST
   ndim=3;
   pio2=1.5707963;
TYPE
   glnp = ARRAY [1..ndim] OF real;
   glndim = glnp;
   glnarray = glnp;
VAR
   ncom : integer;
   pcom,xicom : glndim;
   fret,sr2,x : real;
   i,j : integer;
   p,xi : glnp;

FUNCTION fnc(x : glnp) : real;
(* calling routine must define type
TYPE
   glnp = ARRAY [1..np] OF real;
where np dimension of vector x. *)
VAR
   i : integer;
   f : real;
BEGIN
   f := 0.0;
   FOR i := 1 to 3 DO BEGIN
      f := f+sqr(x[i]-1.0)
   END;
   fnc := f
END;

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

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

(*$I MNBRAK.PAS *)

(*$I BRENT.PAS *)

(*$I LINMIN.PAS *)

BEGIN
   writeln;
   writeln('Minimum of a 3-d quadratic centered');
   writeln('at (1.0,1.0,1.0). Minimum is found');
   writeln('along a series of radials.');
   writeln;
   writeln('x':9,'y':12,'z':12,'minimum':14);
   FOR i := 0 to 10 DO BEGIN
      x := pio2*i/10.0;
      sr2 := sqrt(2.0);
      xi[1] := sr2*cos(x);
      xi[2] := sr2*sin(x);
      xi[3] := 1.0;
      p[1] := 0.0;
      p[2] := 0.0;
      p[3] := 0.0;
      linmin(p,xi,ndim,fret);
      FOR j := 1 to 3 DO write(p[j]:12:6);
      writeln(fret:12:6)
   END
END.
PROGRAM d3r5 (input,output);
(* driver for routine LOCATE *)
CONST
   n=100;
TYPE
   glnarray = ARRAY [1..n] OF real;
VAR
   i,j : integer;
   x : real;
   xx : glnarray;

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

BEGIN
(* create array to be searched *)
   FOR i := 1 to n DO BEGIN
      xx[i] := exp(i/20.0)-74.0
   END;
   writeln ('result of:  j := 0 indicates x too small');
   writeln (' ':12,'j := 100 indicates x too large');
   writeln;
   writeln ('locate ':10,'j':6,'xx(j)':11,'xx(j+1)':12);
(* do test *)
   FOR i := 1 to 19 DO BEGIN
      x := -100.0+200.0*i/20.0;
      locate(xx,n,x,j);
      IF ((j<n) AND (j>0)) THEN BEGIN
         writeln (x:10:4,j:6,xx[j]:12:6,xx[j+1]:12:6)
      END ELSE IF (j=n) THEN BEGIN
         writeln (x:10:4,j:6,xx[j]:12:6,'   upper lim')
      END ELSE BEGIN
         writeln (x:10:4,j:6,'   lower lim',xx[j+1]:12:6)
      END
   END
END.
PROGRAM d2r3(input,output,dfile);
(* driver for routine LUBKSB *)
LABEL 10,99;
CONST
   np=20;
TYPE
   glnpbynp=ARRAY [1..np,1..np] OF real;
   glnarray=ARRAY [1..np] OF real;
   glindx=ARRAY [1..np] OF integer;
VAR
   j,k,l,m,n : integer;
   p : real;
   a,b,c : glnpbynp;
   indx : glindx;
   x : glnarray;
   dfile : text;

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

(*$I LUBKSB.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 matrix a for later testing *)
   FOR l := 1 to n DO BEGIN
      FOR k := 1 to n DO BEGIN
         c[k,l] := a[k,l]
      END
   END;
(* do lu decomposition *)
   ludcmp(c,n,np,indx,p);
(* solve equations for each right-hand vector *)
   FOR k := 1 to m DO BEGIN
      FOR l := 1 to n DO BEGIN
         x[l] := b[l,k]
      END;
      lubksb(c,n,np,indx,x);
(* test results with original matrix *)
      writeln('right-hand side vector:');
      FOR l := 1 to n-1 DO write(b[l,k]:12:6);
      writeln(b[n,k]:12:6); 
      writeln ('result of matrix applied to sol''n vector');
      FOR l := 1 to n DO BEGIN
         b[l,k] := 0.0;
         FOR j := 1 to n DO BEGIN
            b[l,k] := b[l,k]+a[l,j]*x[j]
         END
      END;
      FOR l := 1 to n-1 DO write(b[l,k]:12:6);
      writeln(b[n,k]:12:6);
      writeln('***********************************')
   END;
   IF eof(dfile) THEN GOTO 99;
   writeln('press RETURN for next problem:');
   readln;
   GOTO 10;
99:   close(dfile)
END.
PROGRAM d2r2 (input,output,dfile);
(* driver for routine LUDCMP *)
LABEL 10,99;
CONST
   np=20;
TYPE
   glnpbynp = ARRAY [1..np,1..np] OF real;
   glnarray = ARRAY [1..np] OF real;
   glindx = ARRAY [1..np] OF integer;
VAR
   j,k,l,m,n,dum : integer;
   d : real;
   a,xl,xu,x : glnpbynp;
   indx,jndx : glindx;
   dfile : text;

(*$I MODFILE.PAS *)
(*$I LUDCMP.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,x[k,l]);
      readln(dfile,x[k,l])
   END;
(* print out a-matrix for comparison with product of lower *)
(* and upper decomposition matrices *)
   writeln('original matrix:');
   FOR k := 1 to n DO BEGIN
      FOR l := 1 to n-1 DO write(a[k,l]:12:6);
      writeln(a[k,n]:12:6)
   END;
(* perform the decomposition *)
   ludcmp(a,n,np,indx,d);
(* compose separately the lower and upper matrices *)
   FOR k := 1 to n DO BEGIN
      FOR l := 1 to n DO BEGIN
         IF (l > k) THEN BEGIN
            xu[k,l] := a[k,l];
            xl[k,l] := 0.0
         END ELSE IF (l < k) THEN BEGIN
            xu[k,l] := 0.0;
            xl[k,l] := a[k,l]
         END ELSE BEGIN
            xu[k,l] := a[k,l];
            xl[k,l] := 1.0
         END
      END
   END;
(* compute product of lower and upper matrices for *)
(* comparison with original matrix *)
   FOR k := 1 to n DO BEGIN
      jndx[k] := k;
      FOR l := 1 to n DO BEGIN
         x[k,l] := 0.0;
         FOR j := 1 to n DO BEGIN
            x[k,l] := x[k,l]+xl[k,j]*xu[j,l]
         END
      END
   END;
   writeln('product of lower and upper matrices (rows unscrambled):');
   FOR k := 1 to n DO BEGIN
      dum := jndx[indx[k]];
      jndx[indx[k]] := jndx[k];
      jndx[k] := dum
   END;
   FOR k := 1 to n DO BEGIN
      FOR j := 1 to n DO BEGIN
         IF (jndx[j] = k) THEN BEGIN
            FOR l := 1 to n-1 DO write(x[j,l]:12:6);
            writeln(x[j,n]:12:6)
         END
      END
   END;
   writeln('lower matrix of the decomposition:');
   FOR k := 1 to n DO BEGIN
      FOR l := 1 to n-1 DO write(xl[k,l]:12:6);
      writeln(xl[k,n]:12:6)
   END;
   writeln('upper matrix of the decomposition:');
   FOR k := 1 to n DO BEGIN
      FOR l := 1 to n-1 DO write(xu[k,l]:12:6);
      writeln(xu[k,n]:12:6)
   END;
   writeln('***********************************');
   IF eof(dfile) THEN GOTO 99;
   writeln('press return FOR next problem:');
   readln;
   GOTO 10;
99:   close(dfile)
END.
