FUNCTION f1dim(x: real): real;
(* Programs using F1DIM must declare the variables
TYPE
   glndim = ARRAY [1..ndim] OF real;
VAR
   ncom: integer;
   pcom,xicom: glndim;
and assign values to them externally. *)
VAR
   j: integer;
   xt: glndim;
BEGIN
   FOR j := 1 to ncom DO BEGIN
      xt[j] := pcom[j]+x*xicom[j]
   END;
   f1dim := fnc(xt)
END;
FUNCTION factln(n: integer): real;
(* Programs using routine FACTLN must declare the array
VAR
   gla: ARRAY [1..100] OF real;
and must initialize the array to the values
   FOR i := 1 to 100 DO gla[i] := -1.0;
in the main routine. *)
BEGIN
   IF  (n < 0) THEN BEGIN
      writeln ('pause in FACTLN - negative factorial'); readln END
   ELSE IF (n <= 99) THEN BEGIN
      IF  (gla[n+1] < 0.0) THEN  gla[n+1] := gammln(n+1.0);
      factln := gla[n+1] END
   ELSE BEGIN
      factln := gammln(n+1.0)
   END
END;
FUNCTION factrl(n: integer): real;
(* Programs using routing FACTRL must declare the variables
VAR
   glntop: integer;
   gla: ARRAY [1..33] OF real;
and initialize the values
   glntop := 0;
   gla[1] := 1.0;
in the main routine. *)
VAR
   j: integer;
BEGIN
   IF  (n < 0)  THEN BEGIN
      writeln('pause in FACTRL - negative factorial'); readln END
   ELSE IF (n <= glntop) THEN BEGIN
      factrl := gla[n+1] END
   ELSE IF (n <= 32) THEN BEGIN
      FOR j := glntop+1 to n DO BEGIN
         gla[j+1] := j*gla[j]
      END;
      glntop := n;
      factrl := gla[n+1]
   END ELSE BEGIN
      factrl := exp(gammln(n+1.0))
   END
END;
PROCEDURE fgauss(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;
in the main routine.   *)
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;
PROCEDURE fit(x,y: glndata; ndata: integer; sig: glndata; mwt: integer;
         VAR a,b,siga,sigb,chi2,q: real);
(* Programs using routine FIT must define the type
TYPE
   glndata = ARRAY [1..ndata] OF real;
in the main routine.   *)
VAR
   i: integer;
   wt,t,sy,sxoss,sx,st2,ss,sigdat: real;
BEGIN
   sx := 0.0;
   sy := 0.0;
   st2 := 0.0;
   b := 0.0;
   IF (mwt <> 0)THEN BEGIN
      ss := 0.0;
      FOR i := 1 to ndata DO BEGIN
         wt := 1.0/sqr(sig[i]);
         ss := ss+wt;
         sx := sx+x[i]*wt;
         sy := sy+y[i]*wt
      END
   END ELSE BEGIN
      FOR i := 1 to ndata DO BEGIN
         sx := sx+x[i];
         sy := sy+y[i]
      END;
      ss := ndata
   END;
   sxoss := sx/ss;
   IF (mwt <> 0)THEN BEGIN
      FOR i := 1 to ndata DO BEGIN
         t := (x[i]-sxoss)/sig[i];
         st2 := st2+t*t;
         b := b+t*y[i]/sig[i]
      END
   END ELSE BEGIN
      FOR i := 1 to ndata DO BEGIN
         t := x[i]-sxoss;
         st2 := st2+t*t;
         b := b+t*y[i]
      END
   END;
   b := b/st2;
   a := (sy-sx*b)/ss;
   siga := sqrt((1.0+sx*sx/(ss*st2))/ss);
   sigb := sqrt(1.0/st2);
   chi2 := 0.0;
   IF (mwt = 0)THEN BEGIN
      FOR i := 1 to ndata DO BEGIN
         chi2 := chi2+sqr(y[i]-a-b*x[i])
      END;
      q := 1.0;
      sigdat := sqrt(chi2/(ndata-2));
      siga := siga*sigdat;
      sigb := sigb*sigdat
   END ELSE BEGIN
      FOR i := 1 to ndata DO BEGIN
         chi2 := chi2+sqr((y[i]-a-b*x[i])/sig[i])
      END;
      q := gammq(0.5*(ndata-2),0.5*chi2)
   END;
END;
PROCEDURE fixrts(VAR d: glnparray; npoles: integer);
(* Programs using routine FIXRTS must define the type
TYPE
   glnparray = ARRAY [1..npoles] OF real;
   glcarray = ARRAY [1..2*npoles+2] OF real;
in the main routine. *)
VAR
   j,i: integer;
   size: real;
   polish: boolean;
   a,roots: glcarray;
BEGIN
   a[2*npoles+1] := 1.0;
   a[2*npoles+2] := 0.0;
   FOR j := npoles DOWNTO 1 DO BEGIN
      a[2*j-1] := -d[npoles+1-j];
      a[2*j] := 0.0
   END;
   polish := true;
   zroots(a,npoles,roots,polish);
   FOR j := 1 to npoles DO BEGIN
      size := sqr(roots[2*j-1])+sqr(roots[2*j]);
      IF (size > 1.0) THEN BEGIN
         roots[2*j-1] := roots[2*j-1]/size;
         roots[2*j] := roots[2*j]/size
      END
   END;
   a[1] := -roots[1];
   a[2] := -roots[2];
   a[3] := 1.0;
   a[4] := 0.0;
   FOR j := 2 to npoles DO BEGIN
      a[2*j+1] := 1.0;
      a[2*j+2] := 0.0;
      FOR i := j DOWNTO 2 DO BEGIN
         dum := a[2*i-1];
         a[2*i-1] := a[2*i-3]-a[2*i-1]*roots[2*j-1]
               +a[2*i]*roots[2*j];
         a[2*i] := a[2*i-2]-dum*roots[2*j]
               -a[2*i]*roots[2*j-1]
      END;
      dum := a[1];
      a[1] := -a[1]*roots[2*j-1]+a[2]*roots[2*j];
      a[2] := -dum*roots[2*j]-a[2]*roots[2*j-1]
   END;
   FOR j := 1 to npoles DO BEGIN
      d[npoles+1-j] := -a[2*j-1]
   END
END;
PROCEDURE fleg(x: real; VAR pl: glnlarray; nl: integer);
(* Programs using routine FLEG must define the type
TYPE
   glnlarray = ARRAY [1..nl] OF real;
in the main routine. *)
VAR
   j: integer;
   twox,f2,f1,d: real;
BEGIN
   pl[1] := 1.0;
   pl[2] := x;
   IF (nl > 2) THEN BEGIN
      twox := 2.0*x;
      f2 := x;
      d := 1.0;
      FOR j := 3 to nl DO BEGIN
         f1 := d;
         f2 := f2+twox;
         d := d+1.0;
         pl[j] := (f2*pl[j-1]-f1*pl[j-2])/d
      END
   END
END;
PROCEDURE flmoon(n,nph: integer; VAR jd: integer; VAR frac: real);
VAR
   i: integer;
   rad,xtra,t2,t,c,as,am: real;
BEGIN
   rad := 3.14159265/180.0;
   c := n+nph/4.0;
   t := c/1236.85;
   t2 := sqr(t);
   as := 359.2242+29.105356*c;
   am := 306.0253+385.816918*c+0.010730*t2;
   jd := 2415020+28*n+7*nph;
   xtra := 0.75933+1.53058868*c+(1.178e-4-1.55e-7*t)*t2;
   IF ((nph = 0) OR (nph = 2)) THEN BEGIN
      xtra := xtra+(0.1734-3.93e-4*t)*sin(rad*as)-0.4068*sin(rad*am)
   END ELSE IF ((nph = 1) OR (nph = 3)) THEN BEGIN
      xtra := xtra+(0.1721-4.0e-4*t)*sin(rad*as)-0.6280*sin(rad*am)
   END ELSE BEGIN
      writeln('pause in FLMOON - nph is unknown.'); readln
   END;
   IF (xtra >= 0.0) THEN i := trunc(xtra) ELSE i := trunc(xtra-1.0);
   jd := jd+i;
   frac := xtra-i
END;
PROCEDURE four1(VAR data: gldarray; nn,isign: integer);
(* Programs using routine FOUR1 must define type
TYPE
   gldarray = ARRAY [1..nn2] OF real;
in the calling routine, where nn2=nn+nn. *)
VAR
   ii,jj,n,mmax,m,j,istep,i: integer;
   wtemp,wr,wpr,wpi,wi,theta: double;
   tempr,tempi: real;
BEGIN
   n := 2*nn;
   j := 1;
   FOR ii := 1 to nn DO BEGIN
      i := 2*ii-1;
      IF (j > i) THEN BEGIN
         tempr := data[j];
         tempi := data[j+1];
         data[j] := data[i];
         data[j+1] := data[i+1];
         data[i] := tempr;
         data[i+1] := tempi
      END;
      m := n DIV 2;
      WHILE ((m >= 2) AND (j > m))  DO BEGIN
         j := j-m;
         m := m DIV 2
      END;
      j := j+m
   END;
   mmax := 2;
   WHILE (n > mmax) DO BEGIN
      istep := 2*mmax;
      theta := 6.28318530717959/(isign*mmax);
      wpr := -2.0*sqr(sin(0.5*theta));
      wpi := sin(theta);
      wr := 1.0;
      wi := 0.0;
      FOR ii := 1 to (mmax DIV 2) DO BEGIN
         m := 2*ii-1;
         FOR jj := 0 to ((n-m) DIV istep) DO BEGIN
            i := m + jj*istep;
            j := i+mmax;
            tempr := sngl(wr)*data[j]-sngl(wi)*data[j+1];
            tempi := sngl(wr)*data[j+1]+sngl(wi)*data[j];
            data[j] := data[i]-tempr;
            data[j+1] := data[i+1]-tempi;
            data[i] := data[i]+tempr;
            data[i+1] := data[i+1]+tempi
         END;
         wtemp := wr;
         wr := wr*wpr-wi*wpi+wr;
         wi := wi*wpr+wtemp*wpi+wi
      END;
      mmax := istep
   END
END;
PROCEDURE fourn(VAR data: gldarray; nn: glnnarray; ndim,isign: integer);
(* Programs using routine FOURN must define the types
TYPE
   gldarray = ARRAY [1..ndat2] OF real;
   glnnarray = ARRAY [1..ndim] OF integer;
where ndat2 is twice the product of the transform lengths nn(i). *)
VAR
   i1,i2,i3,i2rev,i3rev,ibit,idim: integer;
   ip1,ip2,ip3,ifp1,ifp2,k1,k2,n: integer;
   ii1,ii2,ii3: integer;
   nprev,nrem,ntot: integer;
   tempi,tempr: real;
   theta,wi,wpi,wpr,wr,wtemp: double;
BEGIN
   ntot := 1;
   FOR idim := 1 to ndim DO BEGIN
      ntot := ntot*nn[idim]
   END;
   nprev := 1;
   FOR idim := 1 to ndim DO BEGIN
      n := nn[idim];
      nrem := ntot DIV (n*nprev);
      ip1 := 2*nprev;
      ip2 := ip1*n;
      ip3 := ip2*nrem;
      i2rev := 1;
      FOR ii2 := 0 to ((ip2-1) DIV ip1) DO BEGIN
         i2 := 1+ii2*ip1;
         IF (i2 < i2rev) THEN BEGIN
            FOR ii1 := 0 to ((ip1-2) DIV 2) DO BEGIN
               i1 := i2+ii1*2;
               FOR ii3 := 0 to ((ip3-i1) DIV ip2) DO BEGIN
                  i3 := i1+ii3*ip2;
                  i3rev := i2rev+i3-i2;
                  tempr := data[i3];
                  tempi := data[i3+1];
                  data[i3] := data[i3rev];
                  data[i3+1] := data[i3rev+1];
                  data[i3rev] := tempr;
                  data[i3rev+1] := tempi
               END
            END
         END;
         ibit := ip2 DIV 2;
         WHILE ((ibit >= ip1) AND (i2rev > ibit)) DO BEGIN
            i2rev := i2rev-ibit;
            ibit := ibit DIV 2
         END;
         i2rev := i2rev+ibit
      END;
      ifp1 := ip1;
      WHILE (ifp1 < ip2) DO BEGIN
         ifp2 := 2*ifp1;
         theta := isign*6.28318530717959/(ifp2 DIV ip1);
         wpr := -2.0*sqr(sin(0.5*theta));
         wpi := sin(theta);
         wr := 1.0;
         wi := 0.0;
         FOR ii3 := 0 to ((ifp1-1) DIV ip1) DO BEGIN
            i3 := 1+ii3*ip1;
            FOR ii1 := 0 to ((ip1-2) DIV 2) DO BEGIN
               i1 := i3+ii1*2;
               FOR ii2 := 0 to ((ip3-i1) DIV ifp2) DO BEGIN
                  i2 := i1+ii2*ifp2;
                  k1 := i2;
                  k2 := k1+ifp1;
                  tempr := sngl(wr)*data[k2]
                     -sngl(wi)*data[k2+1];
                  tempi := sngl(wr)*data[k2+1]
                     +sngl(wi)*data[k2];
                  data[k2] := data[k1]-tempr;
                  data[k2+1] := data[k1+1]-tempi;
                  data[k1] := data[k1]+tempr;
                  data[k1+1] := data[k1+1]+tempi
               END
            END;
            wtemp := wr;
            wr := wr*wpr-wi*wpi+wr;
            wi := wi*wpr+wtemp*wpi+wi
         END;
         ifp1 := ifp2
      END;
      nprev := n*nprev
   END
END;
PROCEDURE fpoly(x: real; VAR p: glnparray; np: integer);
(* Programs using routine FPOLY must define the type
TYPE
   glnparray = ARRAY [1..np] OF real;
in the main routine.   *)
VAR
   j: integer;
BEGIN
   p[1] := 1.0;
   FOR j := 2 to np DO BEGIN
      p[j] := p[j-1]*x
   END
END;
PROCEDURE frprmn(VAR p: glnarray; n: integer; ftol: real;
       VAR iter: integer; VAR fret: real);
(* Programs using routine FRPRMN must supply a
FUNCTION fnc(p: glnarray):real; and a
PROCEDURE dfnc(p: glnarray; VAR g: glnarray);
which evaluate a function and its gradient. They must
also define the type
TYPE
   glnarray = ARRAY [1..n] OF real;
in the main routine. *)
LABEL 99;
CONST
   itmax=200;
   eps=1.0e-10;
VAR
   j,its: integer;
   gg,gam,fp,dgg: real;
   g,h,xi: glnarray;
BEGIN
   fp := fnc(p);
   dfnc(p,xi);
    FOR j := 1 to n DO BEGIN
      g[j] := -xi[j];
      h[j] := g[j];
      xi[j] := h[j]
   END;
   FOR its := 1 to itmax DO BEGIN
      iter := its;
      linmin(p,xi,n,fret);
      IF ((2.0*abs(fret-fp)) <= (ftol*(abs(fret)+abs(fp)+eps)))
         THEN GOTO 99;
      fp := fnc(p);
      dfnc(p,xi);
      gg := 0.0;
      dgg := 0.0;
      FOR j := 1 to n DO BEGIN
         gg := gg+sqr(g[j]);
(*         dgg := dgg+sqr(xi[j])   *)
         dgg := dgg+(xi[j]+g[j])*xi[j]
      END;
      IF (gg = 0.0) THEN GOTO 99;
      gam := dgg/gg;
      FOR j := 1 to n DO BEGIN
         g[j] := -xi[j];
         h[j] := g[j]+gam*h[j];
         xi[j] := h[j]
      END
   END;
   writeln('pause in routine FRPRMN');
   writeln('too many iterations'); readln;
99:   END;

PROCEDURE ftest(data1: gln1array; n1: integer;
      data2: gln2array; n2: integer; VAR f,prob: real);
(* Programs using routine FTEST must define the types
TYPE
   gln1array = ARRAY [1..n1] OF real;
   gln2array = ARRAY [1..n2] OF real;
in the main routine, with n2 less than or equal to n1. *)
VAR
   i: integer;
   var2,var1,df2,df1,ave2,ave1: real;
BEGIN
   IF (n2 > n1) THEN BEGIN
      writeln('pause in routine FTEST');
      writeln('first array must be the larger'); readln
   END;
   avevar(data1,n1,ave1,var1);
   FOR i := 1 to n2 DO data1[i] := data2[i];
   avevar(data1,n2,ave2,var2);
   IF (var1 > var2) THEN BEGIN
      f := var1/var2;
      df1 := n1-1;
      df2 := n2-1
   END ELSE BEGIN
      f := var2/var1;
      df1 := n2-1;
      df2 := n1-1
   END;
   prob  :=  betai(0.5*df2,0.5*df1,df2/(df2+df1*f))
         +(1.0-betai(0.5*df1,0.5*df2,df1/(df1+df2/f)))
END;

