PROCEDURE caldat(julian: integer; VAR mm,id,iyyy: integer);
CONST
   igreg=2299161;
VAR
   je,jd,jc,jb,jalpha,ja: integer;
BEGIN
   IF (julian >= igreg) THEN BEGIN
      jalpha := trunc(((julian-1867216)-0.25)/36524.25);
      ja := julian+1+jalpha-trunc(0.25*jalpha)
   END ELSE BEGIN
      ja := julian
   END;
   jb := ja+1524;
   jc := trunc(6680.0+((jb-2439870)-122.1)/365.25);
   jd := 365*jc+trunc(0.25*jc);
   je := trunc((jb-jd)/30.6001);
   id := jb-jd-trunc(30.6001*je);
   mm := je-1;
   IF (mm > 12) THEN mm := mm-12;
   iyyy := jc-4715;
   IF (mm > 2) THEN iyyy := iyyy-1;
   IF (iyyy <= 0) THEN iyyy := iyyy-1
END;
FUNCTION cel(qqc,pp,aa,bb: real): real;
LABEL 1;
CONST
   ca=0.0003;
   pio2=1.5707963268;
VAR
   a,b,e,f,g: real;
   em,p,q,qc: real;
BEGIN
   IF (qqc = 0.0) THEN BEGIN
      writeln('pause in routine CEL'); readln END;
   qc := abs(qqc);
   a := aa;
   b := bb;
   p := pp;
   e := qc;
   em := 1.0;
   IF (p > 0.0) THEN BEGIN
      p := sqrt(p);
      b := b/p
   END ELSE BEGIN
      f := qc*qc;
      q := 1.0-f;
      g := 1.0-p;
      f := f-p;
      q := q*(b-a*p);
      p := sqrt(f/g);
      a := (a-b)/g;
      b := -q/(g*g*p)+a*p
   END;
1:   f := a;
   a := a+b/p;
   g := e/p;
   b := b+f*g;
   b := b+b;
   p := g+p;
   g := em;
   em := qc+em;
   IF (abs(g-qc) > (g*ca)) THEN BEGIN
      qc := sqrt(e);
      qc := qc+qc;
      e := qc*em;
      GOTO 1
   END;
   cel := pio2*(b+a*em)/(em*(em+p))
END;
PROCEDURE chder(a,b: real; c: glcarray; VAR cder: glcarray; n: integer);
(* Programs using routine CHDER must define the type glcarray
as in routine CHEBFT. *)
VAR
   j: integer;
   con: real;
BEGIN
   cder[n] := 0.0;
   cder[n-1] := 2*(n-1)*c[n];
   IF (n >= 3) THEN BEGIN
      FOR j := n-2 DOWNTO 1 DO BEGIN
         cder[j] := cder[j+2]+2*j*c[j+1]
      END
   END;
   con := 2.0/(b-a);
   FOR j := 1 to n DO BEGIN
      cder[j] := cder[j]*con
   END
END;
FUNCTION chebev(a,b: real; c: glcarray; m: integer; x:real): real;
(* Programs using routine CHEBEV must define the type
glcarray as in routine CHEBFT. *)
VAR
   d,dd,sv,y,y2: real;
   j: integer;
BEGIN
   IF  (((x-a)*(x-b)) > 0.0) THEN BEGIN
      writeln('pause in CHEBEV - x not in range.'); readln
   END;
   d := 0.0;
   dd := 0.0;
   y := (2.0*x-a-b)/(b-a);
   y2 := 2.0*y;
   FOR j := m DOWNTO 2 DO BEGIN
      sv := d;
      d := y2*d-dd+c[j];
      dd := sv
   END;
   chebev := y*d-dd+0.5*c[1]
END;
PROCEDURE chebft(a,b: real; VAR c: glcarray; n: integer);
(* Programs using the routine CHEBFT must define the type
TYPE
   glcarray=ARRAY [1..n] OF real;
in the calling routine, and must externally define a function
func(y:real):real which is to be analyzed. *)
CONST
   pi=3.141592653589793;
VAR
   k,j: integer;
   y,fac,bpa,bma: real;
   sum: double;
   f: glcarray;
BEGIN
   bma := 0.5*(b-a);
   bpa := 0.5*(b+a);
   FOR k := 1 to n DO BEGIN
      y := sngl(cos(pi*(k-0.5)/n));
      f[k] := func(y*bma+bpa)
   END;
   fac := 2.0/n;
   FOR j := 1 to n DO BEGIN
      sum := 0.0;
      FOR k := 1 to n DO BEGIN
         sum := sum+f[k]*cos((pi*(j-1))*((k-0.5)/n))
      END;
      c[j] := sngl(fac*sum)
   END
END;
PROCEDURE chebpc(c: glcarray;VAR d: glcarray; n: integer);
(* Programs using routine CHEBPC must define the type
glcarray as in routine CHEBFT. *)
VAR
   k,j: integer;
   sv: real;
   dd: glcarray;
BEGIN
   FOR j := 1 to n DO BEGIN
      d[j] := 0.0;
      dd[j] := 0.0
   END;
   d[1] := c[n];
   FOR j := n-1 DOWNTO 2 DO BEGIN
      FOR k := n-j+1 DOWNTO 2 DO BEGIN
         sv := d[k];
         d[k] := 2.0*d[k-1]-dd[k];
         dd[k] := sv
      END;
      sv := d[1];
      d[1] := -dd[1]+c[j];
      dd[1] := sv
   END;
   FOR j := n DOWNTO 2 DO BEGIN
      d[j] := d[j-1]-dd[j]
   END;
   d[1] := -dd[1]+0.5*c[1]
END;
PROCEDURE chint(a,b: real; c: glcarray; VAR cint: glcarray; n: integer);
(* Programs using routine CHINT must define the type glcarray as in
routine CHEBFT. *)
VAR
   j: integer;
   sum,fac,con: real;
BEGIN
   con := 0.25*(b-a);
   sum := 0.0;
   fac := 1.0;
   FOR j := 2 to n-1 DO BEGIN
      cint[j] := con*(c[j-1]-c[j+1])/(j-1);
      sum := sum+fac*cint[j];
      fac := -fac
   END;
   cint[n] := con*c[n-1]/(n-1);
   sum := sum+fac*cint[n];
   cint[1] := 2.0*sum
END;
PROCEDURE chsone(bins,ebins: barray; nbins,knstrn: integer;
       VAR df,chsq,prob: real);
(* Programs using routine CHSONE must define the types
TYPE
   barray = ARRAY [1..nbins] OF real;
in the main routine.   *)
VAR
   j: integer;
BEGIN
   df := nbins-1-knstrn;
   chsq := 0.0;
   FOR j := 1 to nbins DO BEGIN
      IF (ebins[j] <= 0.0) THEN BEGIN
         writeln('pause in CHSONE - bad expected number')
      END;
      chsq := chsq+sqr(bins[j]-ebins[j])/ebins[j]
   END;
   prob := gammq(0.5*df,0.5*chsq)
END;
PROCEDURE chstwo(bins1,bins2: barray; nbins,knstrn: integer;
       VAR df,chsq,prob: real);
(* Programs using routine CHSTWO must define type
TYPE
   barray = ARRAY [1..nbins] OF real;
in the main routine.   *)
VAR
   j: integer;
BEGIN
   df := nbins-1-knstrn;
   chsq := 0.0;
   FOR j := 1 to nbins DO BEGIN
      IF ((bins1[j] = 0.0) AND (bins2[j] = 0.0)) THEN BEGIN
         df := df-1.0
      END ELSE BEGIN
         chsq := chsq+sqr(bins1[j]-bins2[j])/
               (bins1[j]+bins2[j])
      END
   END;
   prob := gammq(0.5*df,0.5*chsq)
END;
PROCEDURE cntab1(nn: narray; ni,nj: integer;
       VAR chisq,df,prob,cramrv,ccc: real);
(* Programs using routine CNTAB1 must define type
TYPE
   narray = ARRAY [1..ni,1..nj] OF integer;
in the calling routine. *)
CONST
   maxi=100;
   maxj=100;
   tiny=1.0e-30;
VAR
   nnj,nni,j,i,min: integer;
   sum,expctd: real;
   sumi: ARRAY[1..maxi] OF real;
   sumj: ARRAY[1..maxj] OF real;
BEGIN
   sum := 0;
   nni := ni;
   nnj := nj;
   FOR i := 1 to ni DO BEGIN
      sumi[i] := 0.0;
      FOR j := 1 to nj DO BEGIN
         sumi[i] := sumi[i]+nn[i,j];
         sum := sum+nn[i,j];
      END;
      IF (sumi[i] = 0.0) THEN nni := nni-1;
   END;
   FOR j := 1 to nj DO BEGIN
      sumj[j] := 0.0;
      FOR i := 1 to ni DO BEGIN
         sumj[j] := sumj[j]+nn[i,j];
      END;
      IF (sumj[j] = 0.0) THEN nnj := nj-1;
   END;
   df := nni*nnj-nni-nnj+1;
   chisq := 0.0;
   FOR i := 1 to ni DO BEGIN
      FOR j := 1 to nj DO BEGIN
         expctd := sumj[j]*sumi[i]/sum;
         chisq := chisq+sqr(nn[i,j]-expctd)/(expctd+tiny)
      END
   END;
   prob := gammq(0.5*df,0.5*chisq);
   IF ((nni-1) < (nnj-1)) THEN BEGIN
      min := nni-1
   END ELSE BEGIN
      min := nnj-1
   END;
   cramrv := sqrt(chisq/(sum*min));
   ccc := sqrt(chisq/(chisq+sum))
END;
PROCEDURE cntab2(nn: narray; ni,nj: integer;
      VAR h,hx,hy,hygx,hxgy,uygx,uxgy,uxy: real);
(* Programs using routine CNTAB2 must define type
TYPE
   narray = ARRAY [1..ni,1..nj] OF integer;
in the calling routine. *)
CONST
   maxi=100;
   maxj=100;
   tiny=1.0e-30;
VAR
   j,i: integer;
   sum,p: real;
   sumi: ARRAY[1..maxi] OF real;
   sumj: ARRAY[1..maxj] OF real;
BEGIN
   sum := 0;
   FOR i := 1 to ni DO BEGIN
      sumi[i] := 0.0;
      FOR j := 1 to nj DO BEGIN
         sumi[i] := sumi[i]+nn[i,j];
         sum := sum+nn[i,j]
      END
   END;
   FOR j := 1 to nj DO BEGIN
      sumj[j] := 0.0;
      FOR i := 1 to ni DO BEGIN
         sumj[j] := sumj[j]+nn[i,j]
      END
   END;
   hx := 0.0;
   FOR i := 1 to ni DO BEGIN
      IF (sumi[i] <> 0.0) THEN BEGIN
         p := sumi[i]/sum;
         hx := hx-p*ln(p)
      END
   END;
   hy := 0.0;
   FOR j := 1 to nj DO BEGIN
      IF (sumj[j] <> 0.0) THEN BEGIN
         p := sumj[j]/sum;
         hy := hy-p*ln(p)
      END
   END;
   h := 0.0;
   FOR i := 1 to ni DO BEGIN
      FOR j := 1 to nj DO BEGIN
         IF (nn[i,j] <> 0) THEN BEGIN
            p := nn[i,j]/sum;
            h := h-p*ln(p)
         END
      END
   END;
   hygx := h-hx;
   hxgy := h-hy;
   uygx := (hy-hygx)/(hy+tiny);
   uxgy := (hx-hxgy)/(hx+tiny);
   uxy := 2.0*(hx+hy-h)/(hx+hy+tiny)
END;
PROCEDURE convlv(data: glnarray; n: integer; respns: glnarray; m: integer;
       isign: integer; VAR ans: gln2array);
(* Programs using routine CONVLV must define the types
TYPE
   glnarray = ARRAY [-m..npm] OF real;
   gln2array = ARRAY [-m..n2pm] OF real;
where n and m are the dimension of the data and the response function
respectively, and npm=n+m, n2pm=2*n+m. NOTE: when used with CONVLV, the data
dimension in FOUR1 and in TWOFFT must be the same as gln2array here.
i.e. TYPE  gldarray = gln2array; gl2narray = gln2array *)
VAR
   no2,i,ii: integer;
   dum,mag2: real;
   fft: gln2array;
BEGIN
   FOR i := 1 to ((m-1) DIV 2) DO BEGIN
      respns[n+1-i] := respns[m+1-i]
   END;
   FOR i := ((m+3) DIV 2) to (n-((m-1) DIV 2)) DO BEGIN
      respns[i] := 0.0
   END;
   twofft(data,respns,fft,ans,n);
   no2 := n DIV 2;
   FOR i := 1 to (no2+1) DO BEGIN
      ii := 2*i;
      IF (isign = 1) THEN BEGIN
         dum := ans[ii-1];
         ans[ii-1] := (fft[ii-1]*ans[ii-1]-fft[ii]*ans[ii])/no2;
         ans[ii] := (fft[ii]*dum+fft[ii-1]*ans[ii])/no2
      END ELSE IF (isign = -1) THEN BEGIN
         IF ((sqr(ans[ii-1])+sqr(ans[ii])) = 0.0) THEN BEGIN
            writeln('pause in routine CONVLV');
            writeln('deconvolving at response zero'); readln
         END;
         dum := ans[ii-1];
         mag2 := sqr(ans[ii-1])+sqr(ans[ii]);
         ans[ii-1] := fft[ii-1]*ans[ii-1]+fft[ii]*ans[ii]
               /mag2/no2;
         ans[ii] := fft[ii]*dum-fft[ii-1]*ans[ii]
               /mag2/no2
      END ELSE BEGIN
         writeln('pause in routine CONVLV');
         writeln('no meaning for ISIGN'); readln
      END
   END;
   ans[2] := ans[n+1];
   realft(ans,no2,-1)
END;
PROCEDURE correl(data1,data2: glnarray; n: integer; VAR ans: gl2narray);
(* Programs using routine CORREL must define the type
TYPE
   glnarray = ARRAY [1..n] OF real;
   gl2narray = ARRAY [1..2*n] OF real;
in the main routine. *)
VAR
   no2,i,ii: integer;
   dum: real;
   fft: gl2narray;
BEGIN
   twofft(data1,data2,fft,ans,n);
   no2 := n DIV 2;
   FOR i := 1 to (no2+1) DO BEGIN
      ii := 2*i;
      dum := ans[ii-1];
      ans[ii-1] := (fft[ii-1]*ans[ii-1]+fft[ii]*ans[ii])/no2;
      ans[ii] := (fft[ii]*dum-fft[ii-1]*ans[ii])/no2
   END;
   ans[2] := ans[n+1];
   realft(ans,no2,-1)
END;
PROCEDURE cosft(VAR y: glyarray; n,isign: integer);
(* Programs using routine COSFT must define the type
TYPE
   glyarray = ARRAY [1..n] OF real;
where n is the dimension of the input data array. *)
VAR
   enf0,even,odd,sum,sume,sumo,y1,y2: real;
   theta,wi,wr,wpi,wpr,wtemp: double;
   jj,j,m,n2: integer;
BEGIN
   theta := 3.14159265358979/n;
   wr := 1.0;
   wi := 0.0;
   wpr := -2.0*sqr(sin(0.5*theta));
   wpi := sin(theta);
   sum := y[1];
   m := n DIV 2;
   n2 := n+2;
   FOR j := 2 to (m+1) DO BEGIN
      wtemp := wr;
      wr := wr*wpr-wi*wpi+wr;
      wi := wi*wpr+wtemp*wpi+wi;
      y1 := 0.5*(y[j]+y[n2-j]);
      y2 := (y[j]-y[n2-j]);
      y[j] := y1-sngl(wi)*y2;
      y[n2-j] := y1+sngl(wi)*y2;
      sum := sum+sngl(wr)*y2
   END;
   realft(y,m,+1);
   y[2] := sum;
   FOR jj := 2 to m DO BEGIN
      j := 2*jj;
      sum := sum+y[j];
      y[j] := sum
   END;
   IF (isign = -1) THEN BEGIN
      even := y[1];
      odd := y[2];
      FOR jj := 1 to (m-1) DO BEGIN
         j := 2*jj+1;
         even := even+y[j];
         odd := odd+y[j+1]
      END;
      enf0 := 2.0*(even-odd);
      sumo := y[1]-enf0;
      sume := (2.0*odd/n)-sumo;
      y[1] := 0.5*enf0;
      y[2] := y[2]-sume;
      FOR jj := 1 to (m-1) DO BEGIN
         j := 2*jj+1;
         y[j] := y[j]-sumo;
         y[j+1] := y[j+1]-sume
      END
   END
END;
PROCEDURE covsrt(VAR covar: glcovar; ncvm: integer; ma: integer;
       lista: gllista; mfit: integer);
(* Programs using routine COVSRT must define the types
TYPE
   glcovar = ARRAY [1..ncvm,1..ncvm] OF real;
   gllista = ARRAY [1..mfit] OF integer;
in the calling program. *)
VAR
   j,i: integer;
   swap: real;
BEGIN
   FOR j := 1 to ma-1 DO BEGIN
      FOR i := j+1 to ma DO BEGIN
         covar[i,j] := 0.0
      END
   END;
   FOR i := 1 to mfit-1 DO BEGIN
      FOR j := i+1 to mfit DO BEGIN
         IF (lista[j] > lista[i])  THEN BEGIN
            covar[lista[j],lista[i]] := covar[i,j]
         END ELSE BEGIN
            covar[lista[i],lista[j]] := covar[i,j]
         END
      END
   END;
   swap := covar[1,1];
   FOR j := 1 to ma DO BEGIN
      covar[1,j] := covar[j,j];
      covar[j,j] := 0.0
   END;
   covar[lista[1],lista[1]] := swap;
   FOR j := 2 to mfit DO BEGIN
      covar[lista[j],lista[j]] := covar[1,j]
   END;
   FOR j := 2 to ma DO BEGIN
      FOR i := 1 to j-1 DO BEGIN
         covar[i,j] := covar[j,i]
      END
   END
END;
PROCEDURE crank(n: integer; VAR w: narray; VAR s: real);
(* Programs using routine CRANK must define type
TYPE
   narray = ARRAY [1..n] OF real;
in the calling routine *)
LABEL 2;
VAR
   j,ji,jt,lbl1,lbl2: integer;
   t,rank: real;
BEGIN
   s := 0.0;
   j := 1;
   WHILE (j < n) DO BEGIN
      IF (w[j+1] <> w[j]) THEN BEGIN
         w[j] := j;
         j := j+1
      END ELSE BEGIN
         FOR jt := j+1 to n DO BEGIN
            IF (w[jt] <> w[J]) THEN GOTO 2;
         END;
         jt := n+1;
2:         rank := 0.5*(j+jt-1);
         FOR ji := j to jt-1 DO W[ji] := rank;
         t := jt-j;
         s := s+t*t*t-t;
         j := jt
      END
   END;
   IF (j = n) THEN w[n] := n
END;
