PROCEDURE mdian1(VAR x: narray; n: integer; VAR xmed: real);
(* Program using routine MDIAN1 must define the type
TYPE
   narray = ARRAY [1..n] OF real;
in the calling routine *)
VAR
   n2: integer;
BEGIN
   sort(n,x);
   n2 := n DIV 2;
   IF (2*n2 = n) THEN xmed := 0.5*(x[n2]+x[n2+1])
   ELSE xmed := x[n2+1]
END;
PROCEDURE mdian2(x: narray; n: integer; VAR xmed: real);
(* Programs using routine MDIAN2 must define the type
TYPE
   narray = ARRAY [1..n] OF real;
in the calling routine *)
LABEL 1;
CONST
   big=1.0e30;
   afac=1.5;
   amp=1.5;
VAR
   np,nm,j: integer;
   xx,xp,xm,sumx,sum,eps: real;
   stemp,dum,ap,am,aa,a: real;
BEGIN
   a := 0.5*(x[1]+x[n]);
   eps := abs(x[n]-x[1]);
   ap := big;
   am := -big;
1:   sum := 0.0; sumx := 0.0; np := 0; nm := 0; xp := big; xm := -big;
   FOR j := 1 to n DO BEGIN
      xx := x[j];
      IF (xx <> a) THEN BEGIN
         IF (xx > a) THEN BEGIN
            np := np+1;
            IF (xx < xp) THEN xp := xx END
         ELSE IF  (xx < a) THEN BEGIN
            nm := nm+1;
            IF (xx > xm) THEN xm := xx
         END;
         dum := 1.0/(eps+abs(xx-a));
         sum := sum+dum;
         sumx := sumx+xx*dum
      END
   END;
   stemp := (sumx/sum)-a;
   IF ((np-nm) >= 2) THEN BEGIN
      am := a;
      IF (stemp < 0.0) THEN aa := xp
      ELSE aa := xp+stemp*amp;
      IF (aa > ap) THEN aa := 0.5*(a+ap);
      eps := afac*abs(aa-a);
      a := aa;
      GOTO 1 END
   ELSE IF ((nm-np) >= 2) THEN BEGIN
      ap := a;
      IF (stemp > 0.0) THEN aa := xm
      ELSE aa := xm+stemp*amp;
      IF (aa < am) THEN aa := 0.5*(a+am);
      eps := afac*abs(aa-a);
      a := aa;
      GOTO 1 END
   ELSE IF (n MOD 2) = 0 THEN BEGIN
      IF (np = nm) THEN xmed := 0.5*(xp+xm)
      ELSE IF (np > nm) THEN xmed := 0.5*(a+xp)
      ELSE xmed := 0.5*(xm+a) END
   ELSE BEGIN
      IF (np = nm) THEN xmed := a
      ELSE IF (np > nm) THEN xmed := xp
      ELSE xmed := xm
   END
END;
PROCEDURE medfit(x,y: glndata; ndata: integer; VAR a,b,abdev: real);
(* Programs using routine MEDFIT must define the type
TYPE
   glndata = ARRAY [1..ndata] OF real;
in the main routine. MEDFIT also assumes that the instructions
at the beginning of ROFUNC have been carried out so that arrays
x,y and arr, and real variables aa and abdevt exist as globally
defined variables. *)
LABEL 3;
VAR
   j: integer;
   sy,sxy,sxx,sx,sigb,f2,f1,f: real;
   del,chisq,bb,b2,b1: real;
BEGIN
   sx := 0.0;
   sy := 0.0;
   sxy := 0.0;
   sxx := 0.0;
   FOR j := 1 to ndata DO BEGIN
      sx := sx+x[j];
      sy := sy+y[j];
      sxy := sxy+x[j]*y[j];
      sxx := sxx+sqr(x[j])
   END;
   del := ndata*sxx-sx*sx;
   aa := (sxx*sy-sx*sxy)/del;
   bb := (ndata*sxy-sx*sy)/del;
   chisq := 0.0;
   FOR j := 1 to ndata DO BEGIN
      chisq := chisq+sqr(y[j]-(aa+bb*x[j]))
   END;
   sigb := sqrt(chisq/del);
   b1 := bb;
   f1 := rofunc(b1);
   IF (f1 > 0.0) THEN BEGIN
      b2 := bb+abs(3.0*sigb)
   END ELSE BEGIN
      b2 := bb-abs(3.0*sigb)
   END;
   f2 := rofunc(b2);
   WHILE ((f1*f2) > 0.0) DO BEGIN
      bb := 2.0*b2-b1;
      b1 := b2;
      f1 := f2;
      b2 := bb;
      f2 := rofunc(b2)
   END;
   sigb := 0.01*sigb;
   WHILE (abs(b2-b1) > sigb) DO BEGIN
      bb := 0.5*(b1+b2);
      IF ((bb = b1) OR (bb = b2)) THEN GOTO 3;
      f := rofunc(bb);
      IF ((f*f1) >= 0.0) THEN BEGIN
         f1 := f;
         b1 := bb
      END ELSE BEGIN
         f2 := f;
         b2 := bb
      END
   END;
3:   a := aa;
   b := bb;
   abdev := abdevt/ndata
END;
PROCEDURE memcof(data: glnarray; n,m: integer; VAR pm: real;
      VAR cof: glmarray; wk1,wk2: glnarray; wkm: glmarray);
(* Programs using routine MEMCOF must define the data types
TYPE
   glnarray = ARRAY [1..n] OF real;
   glmarray = ARRAY [1..m] OF real;
and must provide workspace arrays wk1,wk2,wkm with the dimensions
shown in the arguments above. *)
LABEL 99;
VAR
   k,j,i: integer;
   pneum,p,denom: real;
BEGIN
   p := 0.0;
   FOR j := 1 to n DO BEGIN
      p := p+sqr(data[j])
   END;
   pm := p/n;
   wk1[1] := data[1];
   wk2[n-1] := data[n];
   FOR j := 2 to n-1 DO BEGIN
      wk1[j] := data[j];
      wk2[j-1] := data[j]
   END;
   FOR k := 1 to m DO BEGIN
      pneum := 0.0;
      denom := 0.0;
      FOR j := 1 to n-k DO BEGIN
         pneum := pneum+wk1[j]*wk2[j];
         denom := denom+sqr(wk1[j])+sqr(wk2[j])
      END;
      cof[k] := 2.0*pneum/denom;
      pm := pm*(1.0-sqr(cof[k]));
      IF (k <> 1) THEN BEGIN
         FOR i := 1 to k-1 DO BEGIN
            cof[i] := wkm[i]-cof[k]*wkm[k-i]
         END
      END;
      IF (k = m) THEN GOTO 99;
      FOR i := 1 to k DO BEGIN
         wkm[i] := cof[i]
      END;
      FOR j := 1 to n-k-1 DO BEGIN
         wk1[j] := wk1[j]-wkm[k]*wk2[j];
         wk2[j] := wk2[j+1]-wkm[k]*wk1[j+1]
      END
   END;
99:   END;
PROCEDURE midinf(aa,bb: real; VAR s: real; n: integer);
(* Programs using MIDINF must define the function to be integrated
with the declaration FUNCTION func(x: real): real; They must also declare
the global variable
VAR
   glit: integer; *)
VAR
   j: integer;
   x,tnm,sum,del,ddel,b,a: real;
FUNCTION funk(x: real): real;
   BEGIN
      funk := func(1.0/x)/sqr(x)
   END;
BEGIN
   b := 1.0/aa;
   a := 1.0/bb;
   IF (n = 1) THEN BEGIN
      s := (b-a)*funk(0.5*(a+b));
      glit := 1
   END ELSE BEGIN
      tnm := glit;
      del := (b-a)/(3.0*tnm);
      ddel := del+del;
      x := a+0.5*del;
      sum := 0.0;
      FOR j := 1 to glit DO BEGIN
         sum := sum+funk(x);
         x := x+ddel;
         sum := sum+funk(x);
         x := x+del
      END;
      s := (s+(b-a)*sum/tnm)/3.0;
      glit := 3*glit
   END
END;
PROCEDURE midpnt(a,b: real; VAR s: real; n: integer);
(* Programs using routine MIDPNT must supply a function
func(x:real):real which is to be integrated. They must also
declare an iteration counter
VAR
   glit: integer; *)
VAR
   j: integer;
   x,tnm,sum,del,ddel: real;
BEGIN
   IF  (n = 1)  THEN BEGIN
      s := (b-a)*func(0.5*(a+b));
      glit := 1
   END ELSE BEGIN
      tnm := glit;
      del := (b-a)/(3.0*tnm);
      ddel := del+del;
      x := a+0.5*del;
      sum := 0.0;
      FOR j := 1 to glit DO BEGIN
         sum := sum+func(x);
         x := x+ddel;
         sum := sum+func(x);
         x := x+del
      END;
      s := (s+(b-a)*sum/tnm)/3.0;
      glit := 3*glit
   END
END;
PROCEDURE midsql(aa,bb: real; VAR s: real; n: integer);
(* Programs using MIDSQL must define the function to be integrated
with the declaration FUNCTION func(x: real): real; They must also declare
the global variable
VAR
   glit: integer; *)
VAR
   j: integer;
   x,tnm,sum,del,ddel,b,a: real;
FUNCTION funk(x: real): real;
   BEGIN
      funk := 2.0*x*func(aa+sqr(x))
   END;
BEGIN
   b := sqrt(bb-aa);
   a := 0.0;
   IF (n = 1) THEN BEGIN
      s := (b-a)*funk(0.5*(a+b));
      glit := 1
   END ELSE BEGIN
      tnm := glit;
      del := (b-a)/(3.0*tnm);
      ddel := del+del;
      x := a+0.5*del;
      sum := 0.0;
      FOR j := 1 to glit DO BEGIN
         sum := sum+funk(x);
         x := x+ddel;
         sum := sum+funk(x);
         x := x+del
      END;
      s := (s+(b-a)*sum/tnm)/3.0;
      glit := 3*glit
   END
END;
PROCEDURE midsqu(aa,bb: real; VAR s: real; n: integer);
(* Programs using MIDSQU must define the function to be integrated
with the declaration FUNCTION func(x: real): real; They must also declare
the global variable
VAR
   glit: integer; *)
VAR
   j: integer;
   x,tnm,sum,del,ddel,b,a: real;
FUNCTION funk(x: real): real;
   BEGIN
      funk := 2.0*x*func(bb-sqr(x))
   END;
BEGIN
   b := sqrt(bb-aa);
   a := 0.0;
   IF (n = 1) THEN BEGIN
      s := (b-a)*funk(0.5*(a+b));
      glit := 1
   END ELSE BEGIN
      tnm := glit;
      del := (b-a)/(3.0*tnm);
      ddel := del+del;
      x := a+0.5*del;
      sum := 0.0;
      FOR j := 1 to glit DO BEGIN
         sum := sum+funk(x);
         x := x+ddel;
         sum := sum+funk(x);
         x := x+del
      END;
      s := (s+(b-a)*sum/tnm)/3.0;
      glit := 3*glit
   END
END;
PROCEDURE mmid(y,dydx: glnarray; nvar: integer; xs,htot: real;
      nstep: integer; VAR yout: glnarray);
(* Programs using routine MMID 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
   glnarray = ARRAY [1..nvar] OF real;
in the main routine.  Note that this routine is in
single precision, unlike the FORTRAN version. *)
VAR
   n,i: integer;
   x,swap,h2,h: real;
   ym,yn: glnarray;
BEGIN
   h := htot/nstep;
   FOR i := 1 to nvar DO BEGIN
      ym[i] := y[i];
      yn[i] := y[i]+h*dydx[i]
   END;
   x := xs+h;
   derivs(x,yn,yout);
   h2 := 2.0*h;
   FOR n := 2 to nstep DO BEGIN
      FOR i := 1 to nvar DO BEGIN
         swap := ym[i]+h2*yout[i];
         ym[i] := yn[i];
         yn[i] := swap
      END;
      x := x+h;
      derivs(x,yn,yout)
   END;
   FOR i := 1 to nvar DO BEGIN
      yout[i] := 0.5*(ym[i]+yn[i]+h*yout[i])
   END
END;
PROCEDURE mnbrak(VAR ax,bx,cx,fa,fb,fc: real);
(* Programs using routine MNBRAK must supply an external
function func(x:real):real for which a minimum is to be found *)
LABEL 1;
CONST
   gold=1.618034;
   glimit=100.0;
   tiny=1.0e-20;
VAR
   ulim,u,r,q,fu,dum: real;
FUNCTION max(a,b: real): real;
   BEGIN
      IF (a > b) THEN max := a ELSE max := b
   END;
FUNCTION sign(a,b: real): real;
   BEGIN
      IF (b > 0.0) THEN sign := abs(a) ELSE sign := -abs(a)
   END;
BEGIN
   fa := func(ax);
   fb := func(bx);
   IF (fb > fa) THEN BEGIN
      dum := ax;
      ax := bx;
      bx := dum;
      dum := fb;
      fb := fa;
      fa := dum
   END;
   cx := bx+gold*(bx-ax);
   fc := func(cx);
1:   IF (fb >= fc) THEN BEGIN
      r := (bx-ax)*(fb-fc);
      q := (bx-cx)*(fb-fa);
      u := bx-((bx-cx)*q-(bx-ax)*r)/
         (2.0*sign(max(abs(q-r),tiny),q-r));
      ulim := bx+glimit*(cx-bx);
      IF ((bx-u)*(u-cx) > 0.0) THEN BEGIN
         fu := func(u);
         IF (fu < fc) THEN BEGIN
            ax := bx;
            fa := fb;
            bx := u;
            fb := fu;
            GOTO 1 END
         ELSE IF (fu > fb) THEN BEGIN
            cx := u;
            fc := fu;
            GOTO 1
         END;
         u := cx+gold*(cx-bx);
         fu := func(u)
      END ELSE IF  ((cx-u)*(u-ulim) > 0.0) THEN BEGIN
         fu := func(u);
         IF (fu < fc) THEN BEGIN
            bx := cx;
            cx := u;
            u := cx+gold*(cx-bx);
            fb := fc;
            fc := fu;
            fu := func(u)
         END
      END ELSE IF  ((u-ulim)*(ulim-cx) >= 0.0) THEN BEGIN
         u := ulim;
         fu := func(u)
      END ELSE BEGIN
         u := cx+gold*(cx-bx);
         fu := func(u)
      END;
      ax := bx;
      bx := cx;
      cx := u;
      fa := fb;
      fb := fc;
      fc := fu;
      GOTO 1
   END
END;
PROCEDURE mnewt(ntrial: integer; VAR x: glnarray; n: integer;
      tolx,tolf: real);
(* Programs using routine MNEWT must define the types
TYPE
   glnarray = ARRAY [1..n] OF real;
   glnbyn = ARRAY [1..n,1..n] OF real;
   glindx = ARRAY [1..n] OF integer;
in the main routine. They must also supply a routine
PROCEDURE usrfun(x: glnarray; n: integer; VAR alpha: glnbyn;
       VAR beta: glnarray);   *)
LABEL 99;
VAR
   k,i: integer;
   errx,errf,d: real;
   beta: glnarray;
   alpha: glnbyn;
   indx: glindx;
BEGIN
   FOR  k := 1 to ntrial DO BEGIN
      usrfun(x,n,alpha,beta);
      errf := 0.0;
      FOR i := 1 to n DO errf := errf+abs(beta[i]);
      IF (errf <= tolf) THEN GOTO 99;
      ludcmp(alpha,n,np,indx,d);
      lubksb(alpha,n,np,indx,beta);
      errx := 0.0;
      FOR i := 1 to n DO BEGIN
         errx := errx+abs(beta[i]);
         x[i] := x[i]+beta[i]
      END;
      IF (errx <= tolx) THEN GOTO 99
   END;
99:   END;
    TYPE double = real; char12 = string[12];
    FUNCTION sngl(x:real):real; BEGIN sngl := x END;
    PROCEDURE glopen(VAR infile:text; filename:char12);
        BEGIN assign(infile,filename); reset(infile) END;
PROCEDURE moment(data: narray; n: integer;
       VAR ave,adev,sdev,svar,skew,curt: real);
(* Programs using routine MOMENT must define the type
TYPE
   narray = ARRAY [1..n] OF real;
in the calling routine *)
VAR
   j: integer;
   s,p: real;
BEGIN
   IF (n <= 1) THEN BEGIN
      writeln('pause in MOMENT - n must be at least 2'); readln
   END;
   s := 0.0;
   FOR j := 1 to n DO s := s+data[j];
   ave := s/n;
   adev := 0.0;
   svar := 0.0;
   skew := 0.0;
   curt := 0.0;
   FOR j := 1 to n DO BEGIN
      s := data[j]-ave;
      adev := adev+abs(s);
      p := s*s;
      svar := svar+p;
      p := p*s;
      skew := skew+p;
      p := p*s;
      curt := curt+p
   END;
   adev := adev/n;
   svar := svar/(n-1);
   sdev := sqrt(svar);
   IF (svar <> 0.0) THEN BEGIN
      skew := skew/(n*sdev*sdev*sdev);
      curt := curt/(n*sqr(svar))-3.0
   END ELSE BEGIN
      writeln('pause in MOMENT - no skew/kurtosis when variance = 0'); readln
   END
END;
PROCEDURE mprove(a,alud: glnpbynp; n,np: integer; indx: glindx;
       b: glnarray; VAR x: glnarray);
(* Programs using routine MPROVE 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,i: integer;
   sdp: double;
   r: glnarray;
BEGIN
   FOR i := 1 to n DO BEGIN
      sdp := -b[i];
      FOR j := 1 to n DO BEGIN
         sdp := sdp+a[i,j]*x[j]
      END;
      r[i] := sngl(sdp)
   END;
   lubksb(alud,n,np,indx,r);
   FOR i := 1 to n DO BEGIN
      x[i] := x[i]-r[i]
   END
END;
PROCEDURE mrqcof(x,y,sig: glndata; ndata: integer;
       VAR a: glmma; mma: integer; lista: gllista;
       mfit: integer; VAR alpha: glnalbynal;
       VAR beta: glmma; nalp: integer; VAR chisq: real);
(* Programs using routine MRQMIN must provide a
PROCEDURE funcs(xx:real; a:glmma; yfit:real; dyda:glmma; mma:integer);
that evaluates the fitting function yfit and its derivatives dyda
with respect to the parameters a at point xx.  Also they
must define the types
TYPE
   glndata = ARRAY [1..ndata] OF real;
   glmma = ARRAY [1..mma] OF real;
   gllista = ARRAY [1..mma] OF integer;
   glnalbynal = ARRAY [1..nalp,1..nalp] OF real;
in the main routine *)
VAR
   k,j,i: integer;
   ymod,wt,sig2i,dy: real;
   dyda: glmma;
BEGIN
   FOR j := 1 to mfit DO BEGIN
      FOR k := 1 to j DO BEGIN
         alpha[j,k] := 0.0
      END;
      beta[j] := 0.0
   END;
   chisq := 0.0;
   FOR i := 1 to ndata DO BEGIN
      funcs(x[i],a,ymod,dyda,mma);
      sig2i := 1.0/(sig[i]*sig[i]);
      dy := y[i]-ymod;
      FOR j := 1 to mfit DO BEGIN
         wt := dyda[lista[j]]*sig2i;
         FOR k := 1 to j DO BEGIN
            alpha[j,k] := alpha[j,k]+wt*dyda[lista[k]]
         END;
         beta[j] := beta[j]+dy*wt
      END;
      chisq := chisq+dy*dy*sig2i
   END;
   FOR j := 2 to mfit DO BEGIN
      FOR k := 1 to j-1 DO BEGIN
         alpha[k,j] := alpha[j,k]
      END
   END
END;
PROCEDURE mrqmin(x,y,sig: glndata; ndata: integer;
       VAR a: glmma; mma: integer; lista: gllista;
       mfit: integer; VAR covar,alpha: glncabynca;
       nca: integer; VAR chisq,alamda: real);
(* Programs using routine MRQMIN must define the types
TYPE
   glndata = ARRAY [1..ndata] OF real;
   glmma = ARRAY [1..mma] OF real;
   gllista = ARRAY [1..mma] OF integer;
   glncabynca = ARRAY [1..nca,1..nca] OF real;
and the variables
VAR
   glochisq: real;
   glatry,glbeta: glmma;
in the main routine. Also note that this routine calls MRQCOF, which
requires a user-defined procedure FUNCS, described in that routine.   *)
LABEL 99;
VAR
   k,kk,j,ihit: integer;
   da: glmma;
   oneda: glncabynca;
BEGIN
   IF (alamda < 0.0) THEN 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 1 in routine MRQMIN');
            writeln('Improper permutation in LISTA'); readln
         END
      END;
      IF (kk <> (mma+1)) THEN BEGIN
         writeln('pause 2 in routine MRQMIN');
         writeln('Improper permutation in LISTA'); readln
      END;
      alamda := 0.001;
      mrqcof(x,y,sig,ndata,a,mma,lista,mfit,alpha,glbeta,nca,chisq);
      glochisq := chisq;
      FOR j := 1 to mma DO BEGIN
         glatry[j] := a[j]
      END
   END;
   FOR j := 1 to mfit DO BEGIN
      FOR k := 1 to mfit DO BEGIN
         covar[j,k] := alpha[j,k]
      END;
      covar[j,j] := alpha[j,j]*(1.0+alamda);
      oneda[j,1] := glbeta[j]
   END;
   gaussj(covar,mfit,nca,oneda,1,1);
   FOR j := 1 to mfit DO da[j] := oneda[j,1];
   IF (alamda = 0.0) THEN BEGIN
      covsrt(covar,nca,mma,lista,mfit);
      GOTO 99
   END;
   FOR j := 1 to mfit DO BEGIN
      glatry[lista[j]] := glatry[lista[j]]+da[j]
   END;
   mrqcof(x,y,sig,ndata,glatry,mma,lista,mfit,covar,da,nca,chisq);
   IF (chisq < glochisq) THEN BEGIN
      alamda := 0.1*alamda;
      glochisq := chisq;
      FOR j := 1 to mfit DO BEGIN
         FOR k := 1 to mfit DO BEGIN
            alpha[j,k] := covar[j,k]
         END;
         glbeta[j] := da[j];
         a[lista[j]] := glatry[lista[j]]
      END
   END ELSE BEGIN
      alamda := 10.0*alamda;
      chisq := glochisq
   END;
99:   END;
