PROGRAM d13r1(input,output);
(* driver for routine MDIAN1 *)
CONST
   npts=50;
TYPE
   glsarray = ARRAY [1..npts] OF real;
   narray = glsarray;
VAR
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;
   gliset : integer;
   glgset : real;
   i,j,idum : integer;
   xmed : real;
   data : glsarray;

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

(*$I GASDEV.PAS *)

(*$I SORT.PAS *)

(*$I MDIAN1.PAS *)

BEGIN
   gliset := 0;
   idum := -5;
   FOR i := 1 to npts DO BEGIN
      data[i] := gasdev(idum)
   END;
   mdian1(data,npts,xmed);
   writeln('Data drawn from a gaussian distribution');
   writeln('with zero mean and unit variance');
   writeln;
   writeln('Median of data set is',xmed:9:6);
   writeln;
   writeln('Sorted data');
   FOR i := 1 to (npts DIV 5) DO BEGIN
      FOR j := 1 to 5 DO write(data[5*i-5+j]:12:6);
      writeln
   END
END.
PROGRAM d13r2(input,output);
(* driver for routine MDIAN2 *)
CONST
   npts=50;
TYPE
   glsarray = ARRAY [1..npts] OF real;
   narray = glsarray;
VAR
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;
   gliset : integer;
   glgset : real;
   i,idum : integer;
   xmed : real;
   data : narray;

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

(*$I GASDEV.PAS *)

(*$I SORT.PAS *)

(*$I MDIAN1.PAS *)

(*$I MDIAN2.PAS *)

BEGIN
   gliset := 0;
   idum := -5;
   FOR i := 1 to npts DO BEGIN
      data[i] := gasdev(idum)
   END;
   mdian2(data,npts,xmed);
   writeln('Data drawn from a gaussian distribution');
   writeln('with zero mean, unit variance');
   writeln;
   writeln('median according to mdian2 is',xmed:9:6);
   mdian1(data,npts,xmed);
   writeln('median according to mdian1 is',xmed:9:6)
END.
PROGRAM d14r11(input,output);
(* driver for routine MEDFIT *)
CONST
   npt=100;
   spread=0.1;
   ndata=npt;
TYPE
   glndata = ARRAY [1..ndata] OF real;
   glsarray = glndata;
VAR
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;
   gliset : integer;
   glgset : real;
   a,aa,abdev,abdevt,b,chi2,q,siga,sigb : real;
   i,idum,mwt : integer;
   sig,x,y,arr : glndata;

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

(*$I GASDEV.PAS *)

(*$I GAMMLN.PAS *)

(*$I GSER.PAS *)

(*$I GCF.PAS *)

(*$I GAMMQ.PAS *)

(*$I SORT.PAS *)

(*$I FIT.PAS *)

(*$I ROFUNC.PAS *)

(*$I MEDFIT.PAS *)

BEGIN
   gliset := 0;
   idum := -1984;
   FOR i := 1 to npt DO BEGIN
      x[i] := 0.1*i;
      y[i] := -2.0*x[i]+1.0+spread*gasdev(idum);
      sig[i] := spread
   END;
   mwt := 1;
   writeln;
   fit(x,y,npt,sig,mwt,a,b,siga,sigb,chi2,q);
   writeln('According to routine FIT the result is:');
   writeln('a = ':7,a:8:4,'  uncertainty: ':13,siga:8:4);
   writeln('b = ':7,b:8:4,'  uncertainty: ':13,sigb:8:4);
   writeln('chi-squared: ':16,chi2:8:4,' for ',npt:4,' points');
   writeln('goodness-of-fit: ':20,q:8:4);
   writeln;
   writeln('According to routine MEDFIT the result is:');
   medfit(x,y,npt,a,b,abdev);
   writeln('a = ':7,a:8:4);
   writeln('b = ':7,b:8:4);
   writeln(' ':3,'absolute deviation (per data point): ',abdev:8:4);
   writeln(' ':3,'note: gaussian spread is',spread:8:4,')')
END.
PROGRAM d12r9(input,output,dfile);
(* driver for routine MEMCOF *)
CONST
   n=1000;
   m=10;
TYPE
   glnarray = ARRAY [1..n] OF real;
   glmarray = ARRAY [1..m] OF real;
VAR
   i : integer;
   pm : real;
   cof,wkm : glmarray;
   data,wk1,wk2 : glnarray;
   dfile : text;

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

BEGIN
   glopen(dfile,'spctrl.dat');
   FOR i := 1 to n DO BEGIN
      read(dfile,data[i])
   END;
   close(dfile);
   memcof(data,n,m,pm,cof,wk1,wk2,wkm);
   writeln('Coefficients for spectral estimation of SPCTRL.DAT');
   writeln;
   FOR i := 1 to m DO writeln('a[',i:2,'] =',cof[i]:12:6);
   writeln;
   writeln('b0 =',pm:12:6)
END.
PROGRAM d4r5(input,output);
(* driver for routine MIDPNT *)
CONST
   nmax=10;
VAR
   a,b,s : real;
   i,glit : integer;

FUNCTION func(x: real): real;
(* Function for testing integration *)
BEGIN
   func := 1.0/sqrt(x)
END;

FUNCTION fint(x: real): real;
(* Integral of 'func' *)
BEGIN
   fint := 2.0*sqrt(x)
END;

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

BEGIN
   a := 0.0;
   b := 1.0;
   writeln;
   writeln ('Integral of func computed with MIDPNT');
   writeln ('Actual value of integral is',(fint(b)-fint(a)):7:4);
   writeln ('n':6,'Approx. integral':29);
   FOR i := 1 to nmax DO BEGIN
      midpnt(a,b,s,i);
      writeln (i:6,s:24:6)
   END
END.
PROGRAM d15r5(input,output);
(* driver for routine MMID *)
CONST
   nvar=4;
   x1=1.0;
   htot=0.5;
TYPE
   glnarray = ARRAY [1..nvar] OF real;
VAR
   b1,b2,b3,b4,xf : real;
   i,ii : integer;
   y,yout,dydx : glnarray;

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

(*$I BESSJ1.PAS *)

(*$I BESSJ.PAS *)

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

BEGIN
   y[1] := bessj0(x1);
   y[2] := bessj1(x1);
   y[3] := bessj(2,x1);
   y[4] := bessj(3,x1);
   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];
   xf := x1+htot;
   b1 := bessj0(xf);
   b2 := bessj1(xf);
   b3 := bessj(2,xf);
   b4 := bessj(3,xf);
   writeln('First four Bessel functions:');
   FOR ii := 1 to 10 DO BEGIN
      i := 5*ii;
      mmid(y,dydx,nvar,x1,htot,i,yout);
      writeln;
      writeln('x := ',x1:5:2,
         ' to ',x1+htot:5:2,' in ',i:2,' steps');
      writeln('integration':14,'bessj':9);
      writeln(yout[1]:12:6,b1:12:6);
      writeln(yout[2]:12:6,b2:12:6);
      writeln(yout[3]:12:6,b3:12:6);
      writeln(yout[4]:12:6,b4:12:6);
      writeln('press return to continue...');
      readln
   END
END.
PROGRAM d10r1(input,output);
(* driver for routine MNBRAK *)
VAR
   ax,bx,cx,fa,fb,fc : real;
   i : integer;

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

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

(*$I MNBRAK.PAS *)

BEGIN
   FOR i := 1 to 10 DO BEGIN
      ax := i*0.5;
      bx := (i+1.0)*0.5;
      mnbrak(ax,bx,cx,fa,fb,fc);
      writeln('a':14,'b':12,'c':12);
      writeln('x':3,'  ',ax:12:6,bx:12:6,cx:12:6);
      writeln('f':3,'  ',fa:12:6,fb:12:6,fc:12:6)
   END
END.
PROGRAM d9r13(input,output);
(* driver for routine MNEWT *)
CONST
   ntrial=5;
   tolx=1.0e-6;
   n=4;
   np=n;
   tolf=1.0e-6;
TYPE
   glnarray = ARRAY [1..n] OF real;
   glnbyn = ARRAY [1..n,1..n] OF real;
   glindx = ARRAY [1..n] OF integer;
   glnpbynp = glnbyn;
VAR
   i,j,k,kk : integer;
   xx : real;
   x,beta : glnarray;
   alpha : glnbyn;

PROCEDURE usrfun(x: glnarray; n: integer; VAR alpha: glnbyn;
       VAR beta: glnarray);
(* Programs using routine USRFUN must define the types
TYPE
   glnarray = ARRAY [1..n] OF real;
   glnbyn = ARRAY [1..n,1..n] OF real;
in the main routine. *)
BEGIN
   alpha[1,1] := -2.0*x[1];
   alpha[1,2] := -2.0*x[2];
   alpha[1,3] := -2.0*x[3];
   alpha[1,4] := 1.0;
   alpha[2,1] := 2.0*x[1];
   alpha[2,2] := 2.0*x[2];
   alpha[2,3] := 2.0*x[3];
   alpha[2,4] := 2.0*x[4];
   alpha[3,1] := 1.0;
   alpha[3,2] := -1.0;
   alpha[3,3] := 0.0;
   alpha[3,4] := 0.0;
   alpha[4,1] := 0.0;
   alpha[4,2] := 1.0;
   alpha[4,3] := -1.0;
   alpha[4,4] := 0.0;
   beta[1] := sqr(x[1])+sqr(x[2])+sqr(x[3])-x[4];
   beta[2] := -sqr(x[1])-sqr(x[2])-sqr(x[3])-sqr(x[4])+1.0;
   beta[3] := -x[1]+x[2];
   beta[4] := -x[2]+x[3]
END;

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

(*$I LUDCMP.PAS *)

(*$I MNEWT.PAS *)

BEGIN
   FOR kk := 1 to 2 DO BEGIN
      FOR k := 1 to 3 DO BEGIN
         xx := 0.2*k*(2*kk-3);
         writeln('Starting vector number',k:2);
         FOR i := 1 to 4 DO BEGIN
            x[i] := xx+0.2*i;
            writeln('x[':7,i:1,']  :=  ',x[i]:5:2)
         END;
         writeln;
         FOR j := 1 to ntrial DO BEGIN
            mnewt(1,x,n,tolx,tolf);
            usrfun(x,n,alpha,beta);
            writeln('i':5,'x[i]':13,'f':13);
            FOR i := 1 to n DO BEGIN
               writeln(i:5,x[i]:14:6,-beta[i]:15:6)
            END;
            writeln;
            writeln('press RETURN to continue...');
            readln
         END
      END
   END
END.
PROGRAM d13r0(input,output);
(* driver for routine MOMENT *)
CONST
   pi=3.14159265;
   npts=5000;
   nbin=100;
   nppnb=5100;   (* nppnb=npts+nbin *)
TYPE
   narray = ARRAY [1..nppnb] OF real;
VAR
   adev,ave,curt,sdev,skew : real;
   vrnce,x : real;
   i,j,k,nlim : integer;
   data : narray;

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

BEGIN
   i := 1;
   FOR j := 1 to nbin DO BEGIN
      x := pi*j/nbin;
      nlim := round(sin(x)*pi/2.0*npts/nbin);
      FOR k := 1 to nlim DO BEGIN
         data[i] := x;
         i := i+1
      END
   END;
   writeln('moments of a sinusoidal distribution');
   writeln;
   moment(data,npts,ave,adev,sdev,vrnce,skew,curt);
   writeln('calculated':39,'expected':11);
   writeln;
   writeln('Mean :',' ':19,ave:12:4,pi/2.0:12:4);
   writeln('Average Deviation :',' ':6,adev:12:4,(pi/2.0)-1.0:12:4);
   writeln('Standard Deviation :',' ':5,sdev:12:4,0.683667:12:4);
   writeln('Variance :',' ':15,vrnce:12:4,0.467401:12:4);
   writeln('Skewness :',' ':15,skew:12:4,0.0:12:4);
   writeln('Kurtosis :',' ':15,curt:12:4,-0.806249:12:4)
END.
PROGRAM d2r5(input,output);
(* driver for routine MPROVE *)
CONST
   n=5;
   np=n;
TYPE
   glnarray = ARRAY [1..n] OF real;
   glindx = ARRAY [1..n] OF integer;
   glnpbynp = ARRAY [1..np,1..np] OF real;
VAR
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;
   d : real;
   i,idum,j : integer;
   a,aa : glnpbynp;
   b,x : glnarray;
   indx : glindx;

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

(*$I LUDCMP.PAS *)

(*$I LUBKSB.PAS *)

(*$I MPROVE.PAS *)

BEGIN
   a[1,1] := 1.0; a[1,2] := 2.0; a[1,3] := 3.0; a[1,4] := 4.0; 
   a[1,5] := 5.0; a[2,1] := 2.0; a[2,2] := 3.0; a[2,3] := 4.0;
   a[2,4] := 5.0; a[2,5] := 1.0; a[3,1] := 1.0; a[3,2] := 1.0;
   a[3,3] := 1.0; a[3,4] := 1.0; a[3,5] := 1.0; a[4,1] := 4.0;
   a[4,2] := 5.0; a[4,3] := 1.0; a[4,4] := 2.0; a[4,5] := 3.0;
   a[5,1] := 5.0; a[5,2] := 1.0; a[5,3] := 2.0; a[5,4] := 3.0;
   a[5,5] := 4.0;
   b[1] := 1.0; b[2] := 1.0; b[3] := 1.0; b[4] := 1.0; b[5] := 1.0;
   FOR i := 1 to n DO BEGIN
      x[i] := b[i];
      FOR j := 1 to n DO BEGIN
         aa[i,j] := a[i,j]
      END
   END;
   ludcmp(aa,n,np,indx,d);
   lubksb(aa,n,np,indx,x);
   writeln;
   writeln('Solution vector for the equations:');
   FOR i := 1 to n DO write(x[i]:12:6);
   writeln;
(* now phoney up x and let mprove fit it *)
   idum := -13;
   FOR i := 1 to n DO BEGIN
      x[i] := x[i]*(1.0+0.2*ran3(idum))
   END;
   writeln;
   writeln('Solution vector with noise added:');
   FOR i := 1 to n DO write(x[i]:12:6);
   writeln;
   mprove(a,aa,n,np,indx,b,x);
   writeln;
   writeln('Solution vector recovered by mprove:');
   FOR i := 1 to n DO write(x[i]:12:6);
   writeln
END.
PROGRAM d14r9(input,output);
(* driver for routine MRQCOF *)
CONST
   npt=100;
   mma=6;
   spread=0.1;
TYPE
   glndata = ARRAY [1..npt] OF real;
   glmma = ARRAY [1..mma] OF real;
   glnparam=glmma;
   gllista = ARRAY [1..mma] OF integer;
   glnalbynal = ARRAY [1..mma,1..mma] OF real;
VAR
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;
   gliset : integer;
   glgset : real;
   chisq : real;
   i,j,idum,mfit : integer;
   x,y,sig : glndata;
   a,beta,gues : glmma;
   lista : gllista;
   covar,alpha : glnalbynal;

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

(*$I GASDEV.PAS *)

PROCEDURE funcs(x: real; a: glnparam; VAR y: real;
       VAR dyda: glnparam; na: integer);
(* Programs using routine FGAUSS must define the type
TYPE
   glnparam = ARRAY [1..na] OF real;
where na is the number of parameters. *)
VAR
   i,ii : integer;
   fac,ex,arg : real;
BEGIN
   y := 0.0;
   FOR ii := 1 to (na DIV 3) DO BEGIN
      i := 3*ii-2;
      arg := (x-a[i+1])/a[i+2];
      ex := exp(-sqr(arg));
      fac := a[i]*ex*2.0*arg;
      y := y+a[i]*ex;
      dyda[i] := ex;
      dyda[i+1] := fac/a[i+2];
      dyda[i+2] := fac*arg/a[i+2]
   END
END;

(*$I MRQCOF.PAS *)

BEGIN
   gliset := 0;
   a[1] := 5.0; a[2] := 2.0; a[3] := 3.0;
   a[4] := 2.0; a[5] := 5.0; a[6] := 3.0;
   gues[1] := 4.9; gues[2] := 2.1; gues[3] := 2.9;
   gues[4] := 2.1; gues[5] := 4.9; gues[6] := 3.1;
   idum := -911;
(* first try sum of two gaussians *)
   FOR i := 1 to 100 DO BEGIN
      x[i] := 0.1*i;
      y[i] := 0.0;
      y[i] := y[i]+a[1]*exp(-sqr((x[i]-a[2])/a[3]));
      y[i] := y[i]+a[4]*exp(-sqr((x[i]-a[5])/a[6]));
      y[i] := y[i]*(1.0+spread*gasdev(idum));
      sig[i] := spread*y[i]
   END;
   mfit := mma;
   FOR i := 1 to mfit DO BEGIN
      lista[i] := i
   END;
   FOR i := 1 to mma DO BEGIN
      a[i] := gues[i]
   END;
   mrqcof(x,y,sig,npt,a,mma,lista,mfit,alpha,beta,mma,chisq);
   writeln;
   writeln('matrix alpha');
   FOR i := 1 to mma DO BEGIN
      FOR j := 1 to mma DO write(alpha[i,j]:12:4);
      writeln
   END;
   writeln('vector beta');
   FOR i := 1 to mma DO write(beta[i]:12:4);
   writeln;
   writeln('chi-squared:',chisq:12:4);
   writeln;
(* next fix one line and improve the other *)
   FOR i := 1 to 3 DO BEGIN
      lista[i] := i+3
   END;
   mfit := 3;
   FOR i := 1 to mma DO BEGIN
      a[i] := gues[i]
   END;
   mrqcof(x,y,sig,npt,a,mma,lista,mfit,alpha,beta,mma,chisq);
   writeln('matrix alpha');
   FOR i := 1 to mfit DO BEGIN
      FOR j := 1 to mfit DO write(alpha[i,j]:12:4);
      writeln
   END;
   writeln('vector beta');
   FOR i := 1 to mfit DO write(beta[i]:12:4);
   writeln;
   writeln('chi-squared:',chisq:12:4);
   writeln
END.
PROGRAM D14R8(input,output);
(* driver for routine MRQMIN *)
LABEL 1;
CONST
   npt=100;
   mma=6;
   spread=0.001;
TYPE
   glmma = ARRAY [1..mma] OF real;
   glnparam = glmma;
   gllista = ARRAY [1..mma] OF integer;
   glcovar = ARRAY [1..mma,1..mma] OF real;
   glnalbynal = glcovar;
   glncabynca = glcovar;
   glnpbynp = glcovar;
   glnpbymp = glcovar;
   glndata = ARRAY [1..npt] OF real;
VAR
   glochisq : real;
   glatry,glbeta : glmma;
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;
   gliset : integer;
   glgset : real;
   alamda,chisq,ochisq : real;
   i,idum,itst,j,jj,k,mfit : integer;
   x,y,sig : glndata;
   lista : gllista;
   a,gues : glmma;
   covar,alpha : glcovar;

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

(*$I RAN3.PAS *)

(*$I GASDEV.PAS *)

PROCEDURE funcs(x: real; a: glnparam; VAR y: real;
       VAR dyda: glnparam; na: integer);
(* Programs using routine FGAUSS must define the type
TYPE
   glnparam = ARRAY [1..na] OF real;
where na is the number of parameters. *)
VAR
   i,ii : integer;
   fac,ex,arg : real;
BEGIN
   y := 0.0;
   FOR ii := 1 to (na DIV 3) DO BEGIN
      i := 3*ii-2;
      arg := (x-a[i+1])/a[i+2];
      ex := exp(-sqr(arg));
      fac := a[i]*ex*2.0*arg;
      y := y+a[i]*ex;
      dyda[i] := ex;
      dyda[i+1] := fac/a[i+2];
      dyda[i+2] := fac*arg/a[i+2]
   END
END;

(*$I GAUSSJ.PAS *)

(*$I MRQCOF.PAS *)

(*$I MRQMIN.PAS *)

BEGIN
   gliset := 0;
   a[1] := 5.0; a[2] := 2.0; a[3] := 3.0;
   a[4] := 2.0; a[5] := 5.0; a[6] := 3.0;
   gues[1] := 4.5; gues[2] := 2.2; gues[3] := 2.8;
   gues[4] := 2.5; gues[5] := 4.9; gues[6] := 2.8;
   idum := -911;
   FOR i := 1 to 100 DO BEGIN
      x[i] := 0.1*i;
      y[i] := 0.0;
      FOR jj := 1 to 2 DO BEGIN
         j := 3*jj-2;
         y[i] := y[i]+a[j]*exp(-sqr((x[i]-a[j+1])/a[j+2]))
      END;
      y[i] := y[i]*(1.0+spread*gasdev(idum));
      sig[i] := spread*y[i]
   END;
   mfit := 6;
   FOR i := 1 to mfit DO lista[i] := i;
   alamda := -1;
   FOR i := 1 to mma DO a[i] := gues[i];
   mrqmin(x,y,sig,npt,a,mma,lista,mfit,covar,alpha,
      mma,chisq,alamda);
   k := 1;
   itst := 0;
1:   writeln;
   writeln('Iteration #',k:2,'chi-squared:':17,chisq:10:4,
      'alamda:':10,alamda:9);
   writeln('a[1]':7,'a[2]':8,'a[3]':8,'a[4]':8,'a[5]':8,'a[6]':8);
   FOR i := 1 to 6 DO write(a[i]:8:4);
   writeln;
   k := k+1;
   ochisq := chisq;
   mrqmin(x,y,sig,npt,a,mma,lista,mfit,covar,alpha,
      mma,chisq,alamda);
   IF (chisq > ochisq) THEN BEGIN
      itst := 0
   END ELSE IF (abs(ochisq-chisq) < 0.1) THEN BEGIN
      itst := itst+1
   END;
   IF (itst < 2) THEN GOTO 1;
   alamda := 0.0;
   mrqmin(x,y,sig,npt,a,mma,lista,mfit,covar,alpha,
      mma,chisq,alamda);
   writeln('Uncertainties:');
   FOR i := 1 to 6 DO write(sqrt(covar[i,i]):8:4);
   writeln
END.
