PROCEDURE pcshft(a,b: real; VAR d: glcarray; n: integer);
(* Programs using routine PCSHFT must define the type
glcarray as in routine CHEBFT. *)
VAR
   k,j: integer;
   fac,cnst: real;
BEGIN
   cnst := 2.0/(b-a);
   fac := cnst;
   FOR j := 2 to n DO BEGIN
      d[j] := d[j]*fac;
      fac := fac*cnst
   END;
   cnst := 0.5*(a+b);
   FOR j := 1 to n-1 DO BEGIN
      FOR k := n-1 DOWNTO j DO BEGIN
         d[k] := d[k]-cnst*d[k+1]
      END
   END
END;
PROCEDURE pearsn(x,y: nparray; n: integer; VAR r,prob,z: real);
(* Programs using routine PEARSN must define type
TYPE
   nparray = ARRAY [1..n] OF real;
in the main routine. *)
CONST
   tiny=1.0e-20;
VAR
   j: integer;
   yt,xt,t,syy,sxy,sxx,df,ay,ax: real;
BEGIN
   ax := 0.0;
   ay := 0.0;
   FOR j := 1 to n DO BEGIN
      ax := ax+x[j];
      ay := ay+y[j]
   END;
   ax := ax/n;
   ay := ay/n;
   sxx := 0.0;
   syy := 0.0;
   sxy := 0.0;
   FOR j := 1 to n DO BEGIN
      xt := x[j]-ax;
      yt := y[j]-ay;
      sxx := sxx+sqr(xt);
      syy := syy+sqr(yt);
      sxy := sxy+xt*yt;
   END;
   r := sxy/sqrt(sxx*syy);
   z := 0.5*ln(((1.0+r)+tiny)/((1.0-r)+tiny));
   df := n-2;
   t := r*sqrt(df/(((1.0-r)+tiny)*((1.0+r)+tiny)));
   prob := betai(0.5*df,0.5,df/(df+sqr(t)))
(*   prob := erfcc(abs(z*sqrt(n-1.0))/1.4142136)   *)
END;
PROCEDURE piksr2(n: integer; VAR arr,brr: glsarray);
(* Programs using routine PIKSR2 must define the type
TYPE
   glsarray = ARRAY [1..np] OF real;
in the main routine, with np >= n.   *)
LABEL 10;
VAR
   j,i: integer;
   b,a: real;
BEGIN
   FOR j := 2 to n DO BEGIN
      a := arr[j];
      b := brr[j];
      FOR i := j-1 DOWNTO 1 DO BEGIN
         IF (arr[i] <= a) THEN GOTO 10;
         arr[i+1] := arr[i];
         brr[i+1] := brr[i]
      END;
      i := 0;
10:      arr[i+1] := a;
      brr[i+1] := b
   END
END;
PROCEDURE piksrt(n: integer; VAR arr: glsarray);
(* Programs using routine PIKSRT must define the type
TYPE
   glsarray = ARRAY [1..np] OF real;
in the main routine, with np >= n.   *)
LABEL 10;
VAR
   j,i: integer;
   a: real;
BEGIN
   FOR j := 2 to n DO BEGIN
      a := arr[j];
      FOR i := j-1 DOWNTO 1 DO BEGIN
         IF (arr[i] <= a) THEN GOTO 10;
         arr[i+1] := arr[i]
      END;
      i := 0;
10:      arr[i+1] := a
   END
END;
PROCEDURE pinvs(ie1,ie2,je1,jsf,jc1,k: integer; VAR c: glcarray;
      nci,ncj,nck: integer; VAR s: glsarray; nsi,nsj: integer);
(* Programs using routine PINVS must define the types
TYPE
   glcarray = ARRAY [1..nci,1..ncj,1..nck] OF real;
   glsarray = ARRAY [1..nsi,1..nsj] OF real;
in the main routine. *)
CONST
   zero=0.0;
   one=1.0;
   nmax=10;
VAR
   js1,jpiv,jp,je2,jcoff,j,irow,ipiv,id,icoff,i: integer;
   pivinv,piv,dum,big: real;
   pscl: ARRAY [1..nmax] OF real;
   indxr: ARRAY [1..nmax] OF integer;
BEGIN
   je2 := je1+ie2-ie1;
   js1 := je2+1;
   FOR i := ie1 to ie2 DO BEGIN
      big := zero;
      FOR j := je1 to je2 DO IF (abs(s[i,j]) > big) THEN  big := abs(s[i,j]);
      IF (big = zero) THEN BEGIN
         writeln('pause in routine PINVS');
         writeln('singular matrix - row all 0'); readln
      END;
      pscl[i] := one/big;
      indxr[i] := 0
   END;
   FOR id := ie1 to ie2 DO BEGIN
      piv := zero;
      FOR i := ie1 to ie2 DO BEGIN
         IF (indxr[i] = 0)  THEN BEGIN
            big := zero;
            FOR j := je1 to je2 DO BEGIN
               IF (abs(s[i,j]) > big)  THEN BEGIN
                  jp := j;
                  big := abs(s[i,j])
               END
            END;
            IF (big*pscl[i] > piv)  THEN BEGIN
               ipiv := i;
               jpiv := jp;
               piv := big*pscl[i]
            END
         END
      END;
      IF (s[ipiv,jpiv] = zero) THEN BEGIN
         writeln('pause in routine PINVS');
         writeln('singular matrix'); readln
      END;
      indxr[ipiv] := jpiv;
      pivinv := one/s[ipiv,jpiv];
      FOR j := je1 to jsf DO s[ipiv,j] := s[ipiv,j]*pivinv;
      s[ipiv,jpiv] := one;
      FOR i := ie1 to ie2 DO BEGIN
         IF (indxr[i] <> jpiv)  THEN BEGIN
            IF (s[i,jpiv] <> zero)  THEN BEGIN
               dum := s[i,jpiv];
               FOR j := je1 to jsf DO s[i,j] := s[i,j]-dum*s[ipiv,j];
               s[i,jpiv] := zero
            END
         END
      END
   END;
   jcoff := jc1-js1;
   icoff := ie1-je1;
   FOR i := ie1 to ie2 DO BEGIN
      irow := indxr[i]+icoff;
      FOR j := js1 to jsf DO c[irow,j+jcoff,k] := s[i,j]
   END
END;
FUNCTION plgndr(l,m: integer; x: real): real;
VAR
   fact,pll,pmm,pmmp1,somx2: real;
   i,ll: integer;
BEGIN
   IF ((m < 0) OR (m > l) OR (abs(x) > 1.0)) THEN BEGIN
      writeln('Pause in routine PLGNDR');
      writeln('bad arguments'); readln
   END;
   pmm := 1.0;
   IF (m > 0)  THEN BEGIN
      somx2 := sqrt((1.0-x)*(1.0+x));
      fact := 1.0;
      FOR i := 1 to m DO BEGIN
         pmm := -pmm*fact*somx2;
         fact := fact+2.0
      END
   END;
   IF (l = m) THEN BEGIN
      plgndr := pmm
   END ELSE BEGIN
      pmmp1 := x*(2*m+1)*pmm;
      IF (l = m+1)  THEN BEGIN
         plgndr := pmmp1
      END ELSE BEGIN
         FOR ll := m+2 to l DO BEGIN
            pll := (x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m);
            pmm := pmmp1;
            pmmp1 := pll
         END;
         plgndr := pll
      END
   END
END;
FUNCTION poidev(xm: real; VAR idum: integer): real;
(* Programs using POIDEV must declare the variables
VAR
   gloldm,glsq,glalxm,glg: real;
in the main program and should intialize gloldm to
   gloldm := -1.0;   *)
CONST
   pi=3.141592654;
VAR
   em,t,y: real;
BEGIN
   IF (xm < 12.0) THEN BEGIN
      IF (xm <> gloldm) THEN BEGIN
         gloldm := xm;
         glg := exp(-xm)
      END;
      em := -1;
      t := 1.0;
      REPEAT
         em := em+1.0;
         t := t*ran3(idum);
      UNTIL (t <= glg)
   END ELSE BEGIN
      IF (xm <> gloldm) THEN BEGIN
         gloldm := xm;
         glsq := sqrt(2.0*xm);
         glalxm := ln(xm);
         glg := xm*glalxm-gammln(xm+1.0)
      END;
      REPEAT
         REPEAT
            y := pi*ran3(idum);
            y := sin(y)/cos(y);
            em := glsq*y+xm;
         UNTIL (em >= 0.0);
         em := trunc(em);
         t := 0.9*(1.0+sqr(y))*exp(em*glalxm-gammln(em+1.0)-glg);
      UNTIL (ran3(idum) <= t)
   END;
   poidev := em
END;
PROCEDURE polcoe(x,y: glnarray; n: integer; VAR cof: glnarray);
(* Programs using routine POLCOE must define the type
TYPE
   glnarray = ARRAY [1..narray] OF real;
in the main routine. *)
VAR
   k,j,i: integer;
   phi,ff,b: real;
   s: glnarray;
BEGIN
   FOR i := 1 to n DO BEGIN
      s[i] := 0.0;
      cof[i] := 0.0
   END;
   s[n] := -x[1];
   FOR i := 2 to n DO BEGIN
      FOR j := n+1-i to n-1 DO BEGIN
         s[j] := s[j]-x[i]*s[j+1]
      END;
      s[n] := s[n]-x[i]
   END;
   FOR j := 1 to n DO BEGIN
      phi := n;
      FOR k := n-1 DOWNTO 1 DO BEGIN
         phi := k*s[k+1]+x[j]*phi
      END;
      ff := y[j]/phi;
      b := 1.0;
      FOR k := n DOWNTO 1 DO BEGIN
         cof[k] := cof[k]+b*ff;
         b := s[k]+x[j]*b
      END
   END
END;
PROCEDURE polcof(xa,ya: glnarray; n: integer; VAR cof: glnarray);
(* Programs using routine POLCOF must define the type
TYPE
   glnarray = ARRAY [1..n] OF real;
in the main routine. *)
VAR
   k,j,i: integer;
   xmin,dy: real;
   x,y: glnarray;
BEGIN
   FOR j := 1 to n DO BEGIN
      x[j] := xa[j];
      y[j] := ya[j]
   END;
   FOR j := 1 to n DO BEGIN
      polint(x,y,n+1-j,0.0,cof[j],dy);
      xmin := 1.0E38;
      k := 0;
      FOR i := 1 to n+1-j DO BEGIN
         IF (abs(x[i]) < xmin) THEN BEGIN
            xmin := abs(x[i]);
            k := i
         END;
         IF (x[i] <> 0.0) THEN y[i] := (y[i]-cof[j])/x[i]
      END;
      IF (k < (n+1-j)) THEN BEGIN
         FOR i := k+1 to n+1-j DO BEGIN
            y[i-1] := y[i];
            x[i-1] := x[i]
         END
      END
   END
END;
PROCEDURE poldiv(u: glnarray; n: integer; v: glnvarray; nv: integer;
       VAR q,r: glnarray);
(* Programs using routine POLDIV must define the types
TYPE
   glnarray = ARRAY [1..n] OF real;
   glnvarray = ARRAY [1..nv] OF real;
in the main routine. *)
VAR
   k,j: integer;
BEGIN
   FOR j := 1 to n DO BEGIN
      r[j] := u[j];
      q[j] := 0.0
   END;
   FOR k := n-nv DOWNTO 0 DO BEGIN
      q[k+1] := r[nv+k]/v[nv];
      FOR j := nv+k-1 DOWNTO k+1 DO BEGIN
         r[j] := r[j]-q[k+1]*v[j-k]
      END
   END;
   r[nv] := 0.0
END;
PROCEDURE polin2(x1a: glmarray; x2a: glnarray; ya: glmbyn;
       m,n: integer; x1,x2: real; VAR y,dy: real);
(* Programs using POLIN2 must define the types
TYPE
   glmarray = ARRAY [1..m] OF real;
   glnarray = ARRAY [1..n] OF real;
   glmbyn = ARRAY [1..m,1..n] OF real;
in the main routine. *)
VAR
   k,j: integer;
   ymtmp: glmarray;
   yntmp: glnarray;
BEGIN
   FOR j := 1 to m DO BEGIN
      FOR k := 1 to n DO BEGIN
         yntmp[k] := ya[j,k]
      END;
      polint(x2a,yntmp,n,x2,ymtmp[j],dy)
   END;
   polint(x1a,ymtmp,m,x1,y,dy)
END;
PROCEDURE polint(xa,ya: glnarray; n: integer;
       x: real; VAR y,dy: real);
(* Programs using routine POLINT must define the type
TYPE
   glnarray = ARRAY [1..n] OF real;
in the main routine. *)
VAR
   ns,m,i: integer;
   w,hp,ho,dift,dif,den: real;
   c,d: glnarray;
BEGIN
   ns := 1;
   dif := abs(x-xa[1]);
   FOR i := 1 to n DO BEGIN
      dift := abs(x-xa[i]);
      IF (dift < dif) THEN BEGIN
         ns := i;
         dif := dift
      END;
      c[i] := ya[i];
      d[i] := ya[i]
   END;
   y := ya[ns];
   ns := ns-1;
   FOR m := 1 to n-1 DO BEGIN
      FOR i := 1 to n-m DO BEGIN
         ho := xa[i]-x;
         hp := xa[i+m]-x;
         w := c[i+1]-d[i];
         den := ho-hp;
         IF (den = 0.0) THEN BEGIN
            writeln ('pause in routine POLINT'); readln END;
         den := w/den;
         d[i] := hp*den;
         c[i] := ho*den
      END;
      IF ((2*ns) < (n-m)) THEN BEGIN
         dy := c[ns+1]
      END ELSE BEGIN
         dy := d[ns];
         ns := ns-1
      END;
      y := y+dy
   END
END;
PROCEDURE powell(VAR p: glnarray; VAR xi: glnpbynp; n,np: integer;
       ftol: real; VAR iter: integer; VAR fret: real);
(* Programs using routine POWELL must define the types
TYPE
   glnarray = ARRAY [1..n] OF real;
   glnpbynp = ARRAY [1..np,1..np] OF real;
in the main routine. *)
LABEL 1,99;
CONST
   itmax=200;
VAR
   j,ibig,i: integer;
   t,fptt,fp,del: real;
   pt,ptt,xit: glnarray;
BEGIN
   fret := fnc(p);
   FOR j := 1 to n DO BEGIN
      pt[j] := p[j]
   END;
   iter := 0;
1:   iter := iter+1;
   fp := fret;
   ibig := 0;
   del := 0.0;
   FOR i := 1 to n DO BEGIN
      FOR j := 1 to n DO BEGIN
         xit[j] := xi[j,i]
      END;
      linmin(p,xit,n,fret);
      IF (abs(fp-fret) > del) THEN BEGIN
         del := abs(fp-fret);
         ibig := i
      END
   END;
   IF (2.0*abs(fp-fret) <= ftol*(abs(fp)+abs(fret))) THEN GOTO 99;
   IF (iter = itmax) THEN BEGIN
      writeln('pause in routine POWELL');
      writeln('too many interations'); readln
   END;
   FOR j := 1 to n DO BEGIN
      ptt[j] := 2.0*p[j]-pt[j];
      xit[j] := p[j]-pt[j];
      pt[j] := p[j]
   END;
   fptt := fnc(ptt);
   IF (fptt >= fp) THEN GOTO 1;
   t := 2.0*(fp-2.0*fret+fptt)*sqr(fp-fret-del)-del*sqr(fp-fptt);
   IF (t >= 0.0) THEN GOTO 1;
   linmin(p,xit,n,fret);
   FOR j := 1 to n DO BEGIN
      xi[j,ibig] := xit[j]
   END;
   GOTO 1;
99:   END;
PROCEDURE predic(data: gldarray; ndata: integer; d: glnparray;
      npoles: integer; VAR future: glnfarray; nfut: integer);
(* Programs using routine PREDIC must define the types
TYPE
   gldarray = ARRAY [1..ndata] OF real;
   glnparray = ARRAY [1..npoles] OF real;
   glnfarray = ARRAY [1..nfut] OF real;
in the main routine. *)
VAR
   k,j: integer;
   sum,discrp: real;
   reg: glnparray;
BEGIN
   FOR j := 1 to npoles DO BEGIN
      reg[j] := data[ndata+1-j]
   END;
   FOR j := 1 to nfut DO BEGIN
      discrp := 0.0;
      sum := discrp;
      FOR k := 1 to npoles DO BEGIN
         sum := sum+d[k]*reg[k]
      END;
      FOR k := npoles DOWNTO 2 DO BEGIN
         reg[k] := reg[k-1]
      END;
      reg[1] := sum;
      future[j] := sum
   END
END;
FUNCTION probks(alam: real): real;
LABEL 1;
VAR
   a2,fac,sum,term,termbf: real;
   j: integer;
BEGIN
   a2 := -2.0*alam*alam;
   fac := 2.0;
   sum := 0.0;
   termbf := 0.0;
   FOR j := 1 TO 100 DO BEGIN
      term := fac*exp(a2*sqr(j));
      sum := sum+term;
      IF (abs(term) < (0.001*termbf)) THEN BEGIN
         probks := sum;
         GOTO 1 END
      ELSE BEGIN
         fac := -fac;
         termbf := abs(term)
      END
   END;
   probks := 1.0;
1:   END;
PROCEDURE pzextr(iest: integer; xest: real; yest: glyarray;
       VAR yz,dy: glyarray; nv,nuse: integer);
(* Programs using routine PZEXTR must declare
TYPE
   glyarray = ARRAY [1..nv] OF real;
CONST
   glimax=11;
   glnmax=10;
   glncol=7;
VAR
   glx: ARRAY [1..glimax] OF real;
   glqcol: ARRAY [1..glnmax,1..glncol] OF real;
in the main routine. *)
CONST
   nmax=10;
VAR
   m1,k1,j: integer;
   q,f2,f1,delta: real;
   d: ARRAY [1..nmax] OF real;
BEGIN
   glx[iest] := xest;
   FOR j := 1 to nv DO BEGIN
      dy[j] := yest[j];
      yz[j] := yest[j]
   END;
   IF (iest = 1) THEN BEGIN
      FOR j := 1 to nv DO BEGIN
         glqcol[j,1] := yest[j]
      END
   END ELSE BEGIN
      IF (iest < nuse) THEN m1 := iest ELSE m1 := nuse;
      FOR j := 1 to nv DO BEGIN
         d[j] := yest[j]
      END;
      FOR k1 := 1 to m1-1 DO BEGIN
         delta := 1.0/(glx[iest-k1]-xest);
         f1 := xest*delta;
         f2 := glx[iest-k1]*delta;
         FOR j := 1 to nv DO BEGIN
            q := glqcol[j,k1];
            glqcol[j,k1] := dy[j];
            delta := d[j]-q;
            dy[j] := f1*delta;
            d[j] := f2*delta;
            yz[j] := yz[j]+dy[j]
         END
      END;
      FOR j := 1 to nv DO BEGIN
         glqcol[j,m1] := dy[j]
      END
   END
END;
