FUNCTION gamdev(VAR ia,idum: integer): real;
VAR
   am,e,s,v1,v2,x,y: real;
   j: integer;
BEGIN
   IF (ia < 1) THEN BEGIN
      writeln('pause in routine GAMDEV'); readln END;
   IF (ia < 6) THEN BEGIN
      x := 1.0;
      FOR j := 1 to ia DO x := x*ran3(idum);
      x := -ln(x);
   END ELSE BEGIN
      REPEAT
         REPEAT
            REPEAT
               v1 := 2.0*ran3(idum)-1.0;
               v2 := 2.0*ran3(idum)-1.0;
            UNTIL ((sqr(v1)+sqr(v2)) <= 1.0);
            y := v2/v1;
            am := ia-1;
            s := sqrt(2.0*am+1.0);
            x := s*y+am;
         UNTIL (x > 0.0);
         e := (1.0+sqr(y))*exp(am*ln(x/am)-s*y);
      UNTIL (ran3(idum) <= e)
   END;
   gamdev := x
END;
FUNCTION gammln(xx: real): real;
CONST
   stp = 2.50662827465;
   half = 0.5;
   one = 1.0;
   fpf = 5.5;
VAR
   x,tmp,ser: double;
   j: integer;
   cof: ARRAY [1..6] OF double;
BEGIN
   cof[1] := 76.18009173;
   cof[2] := -86.50532033;
   cof[3] := 24.01409822;
   cof[4] := -1.231739516;
   cof[5] := 0.120858003e-2;
   cof[6] := -0.536382e-5;
   x := xx-one;
   tmp := x+fpf;
   tmp := (x+half)*ln(tmp)-tmp;
   ser := one;
   FOR j := 1 to 6 DO BEGIN
      x := x+one;
      ser := ser+cof[j]/x
   END;
   gammln := sngl(tmp+ln(stp*ser))
END;
FUNCTION gammp(a,x: real): real;
VAR
   gammcf,gln: real;
BEGIN
   IF ((x < 0.0) OR (a <= 0.0)) THEN BEGIN
      writeln('pause in GAMMP - invalid arguments'); readln
   END;
   IF (x < (a+1.0)) THEN BEGIN
      gser(a,x,gammcf,gln);
      gammp := gammcf
   END ELSE BEGIN
      gcf(a,x,gammcf,gln);
      gammp := 1.0-gammcf
   END
END;
FUNCTION gammq(a,x: real): real;
VAR
   gamser,gln: real;
BEGIN
   IF ((x < 0.0) OR (a <= 0.0)) THEN BEGIN
      writeln('pause in GAMMQ - invalid arguments'); readln
   END;
   IF (x < a+1.0) THEN BEGIN
      gser(a,x,gamser,gln);
      gammq := 1.0-gamser
   END ELSE BEGIN
      gcf(a,x,gamser,gln);
      gammq := gamser
   END
END;
FUNCTION gasdev(VAR idum: integer): real;
(* Programs using GASDEV must declare the variables
VAR
   gliset: integer;
   glgset: real;
in the main routine and must intialize gliset to
   gliset := 0;   *)
VAR
   fac,r,v1,v2: real;
BEGIN
   IF  (gliset = 0)  THEN BEGIN
      REPEAT
         v1 := 2.0*ran3(idum)-1.0;
         v2 := 2.0*ran3(idum)-1.0;
         r := sqr(v1)+sqr(v2);
      UNTIL (r < 1.0);
      fac := sqrt(-2.0*ln(r)/r);
      glgset := v1*fac;
      gasdev := v2*fac;
      gliset := 1
   END ELSE BEGIN
      gasdev := glgset;
      gliset := 0
   END
END;
PROCEDURE gauleg(x1,x2: double; VAR x,w: darray; n: integer);
(* Programs using routine GAULEG must define the type
TYPE
   darray=ARRAY [1..n] OF double;
in the calling program *)
CONST
   eps = 3.0e-11; (* adjust to your floating precision *)
VAR
       m,j,i: integer;
       z1,z,xm,xl,pp,p3,p2,p1: double;
BEGIN
   m := (n+1) DIV 2;
   xm := 0.5*(x2+x1);
   xl := 0.5*(x2-x1);
   FOR i := 1 to m DO BEGIN
      z := cos(3.141592654*(i-0.25)/(n+0.5));
      REPEAT
             p1 := 1.0;
             p2 := 0.0;
             FOR j := 1 to n DO BEGIN
               p3 := p2;
               p2 := p1;
               p1 := ((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j
             END;
             pp := n*(z*p1-p2)/(z*z-1.0);
             z1 := z;
             z := z1-p1/pp;
      UNTIL (abs(z-z1) <= eps);
      x[i] := xm-xl*z;
      x[n+1-i] := xm+xl*z;
      w[i] := 2.0*xl/((1.0-z*z)*pp*pp);
      w[n+1-i] := w[i]
   END
END;
PROCEDURE gaussj(VAR a: glnpbynp; n,np: integer;
       VAR b: glnpbymp; m,mp: integer);
(* Programs using GAUSSJ must define the types
TYPE
   glnpbynp = ARRAY [1..np,1..np] OF real;
   glnpbymp = ARRAY [1..np,1..mp] OF real;
   glnp = ARRAY [1..np] OF integer;
in the main routine. *)
VAR
   big,dum,pivinv: real;
   i,icol,irow,j,k,l,ll: integer;
   indxc,indxr,ipiv: glnp;
BEGIN
   FOR j := 1 to n DO BEGIN
      ipiv[j] := 0
   END;
   FOR i := 1 to n DO BEGIN
      big := 0.0;
      FOR j := 1 to n DO BEGIN
         IF (ipiv[j] <> 1) THEN BEGIN
            FOR k := 1 to n DO BEGIN
               IF (ipiv[k] = 0) THEN BEGIN
                  IF (abs(a[j,k]) >= big) THEN BEGIN
                     big := abs(a[j,k]);
                     irow := j;
                     icol := k
                  END
               END ELSE IF (ipiv[k] > 1) THEN BEGIN
                  writeln('pause 1 in GAUSSJ - singular matrix'); readln
               END
            END
         END
      END;
      ipiv[icol] := ipiv[icol]+1;
      IF (irow <> icol) THEN BEGIN
         FOR l := 1 to n DO BEGIN
            dum := a[irow,l];
            a[irow,l] := a[icol,l];
            a[icol,l] := dum
         END;
         FOR l := 1 to m DO BEGIN
            dum := b[irow,l];
            b[irow,l] := b[icol,l];
            b[icol,l] := dum
         END
      END;
      indxr[i] := irow;
      indxc[i] := icol;
      IF (a[icol,icol] = 0.0) THEN BEGIN
         writeln('pause 2 in GAUSSJ - singular matrix'); readln
      END;
      pivinv := 1.0/a[icol,icol];
      a[icol,icol] := 1.0;
      FOR l := 1 to n DO BEGIN
         a[icol,l] := a[icol,l]*pivinv
      END;
      FOR l := 1 to m DO BEGIN
         b[icol,l] := b[icol,l]*pivinv
      END;
      FOR ll := 1 to n DO BEGIN
         IF (ll <> icol) THEN BEGIN
            dum := a[ll,icol];
            a[ll,icol] := 0.0;
            FOR l := 1 to n DO BEGIN
               a[ll,l] := a[ll,l]-a[icol,l]*dum
            END;
            FOR l := 1 to m DO BEGIN
               b[ll,l] := b[ll,l]-b[icol,l]*dum
            END
         END
      END
   END;
   FOR l := n DOWNTO 1 DO BEGIN
      IF (indxr[l] <> indxc[l]) THEN BEGIN
         FOR k := 1 to n DO BEGIN
            dum := a[k,indxr[l]];
            a[k,indxr[l]] := a[k,indxc[l]];
            a[k,indxc[l]] := dum
         END
      END
   END
END;
PROCEDURE gcf(a,x: real; VAR gammcf,gln: real);
LABEL 1;
CONST
   itmax=100;
   eps=3.0e-7;
VAR
   n: integer;
   gold,g,fac,b1,b0,anf,ana,an,a1,a0: real;
BEGIN
   gln := gammln(a);
   gold := 0.0;
   a0 := 1.0;
   a1 := x;
   b0 := 0.0;
   b1 := 1.0;
   fac := 1.0;
   FOR n := 1 to itmax DO BEGIN
      an := 1.0*n;
      ana := an-a;
      a0 := (a1+a0*ana)*fac;
      b0 := (b1+b0*ana)*fac;
      anf := an*fac;
      a1 := x*a0+anf*a1;
      b1 := x*b0+anf*b1;
      IF (a1 <> 0.0) THEN BEGIN
         fac := 1.0/a1;
         g := b1*fac;
         IF (abs((g-gold)/g) < eps) THEN GOTO 1;
         gold := g
      END
   END;
   writeln('pause in GCF - a too large, itmax too small'); readln;
1:   gammcf := exp(-x+a*ln(x)-gln)*g
END;
FUNCTION golden(ax,bx,cx,tol: real; VAR xmin: real): real;
(* Programs using routine GOLDEN must supply an external
function func(x:real):real whose minimum is to be found. *)
CONST
   r=0.61803399;
VAR
   f0,f1,f2,f3,c: real;
   x0,x1,x2,x3: real;
BEGIN
   c := 1.0-r;
   x0 := ax;
   x3 := cx;
   IF (abs(cx-bx) > abs(bx-ax)) THEN BEGIN
      x1 := bx;
      x2 := bx+c*(cx-bx)
   END ELSE BEGIN
      x2 := bx;
      x1 := bx-c*(bx-ax)
   END;
   f1 := func(x1);
   f2 := func(x2);
   WHILE (abs(x3-x0) > tol*(abs(x1)+abs(x2))) DO BEGIN
      IF (f2 < f1) THEN BEGIN
         x0 := x1;
         x1 := x2;
         x2 := r*x1+c*x3;
         f0 := f1;
         f1 := f2;
         f2 := func(x2)
      END ELSE BEGIN
         x3 := x2;
         x2 := x1;
         x1 := r*x2+c*x0;
         f3 := f2;
         f2 := f1;
         f1 := func(x1)
      END
   END;
   IF (f1 < f2) THEN BEGIN
      golden := f1;
      xmin := x1
   END ELSE BEGIN
      golden := f2;
      xmin := x2
   END
END;
PROCEDURE gser(a,x: real; VAR gamser,gln: real);
LABEL 1;
CONST
   itmax=100;
   eps=3.0e-7;
VAR
   n: integer;
   sum,del,ap: real;
BEGIN
   gln := gammln(a);
   IF (x <= 0.0) THEN BEGIN
      IF (x < 0.0) THEN BEGIN
         writeln('pause in GSER - x less than 0'); readln
      END;
      gamser := 0.0
   END ELSE BEGIN
      ap := a;
      sum := 1.0/a;
      del := sum;
      FOR n := 1 to itmax DO BEGIN
         ap := ap+1.0;
         del := del*x/ap;
         sum := sum+del;
         IF (abs(del) < abs(sum)*eps) THEN GOTO 1
      END;
      writeln('pause in GSER - a too large, itmax too small'); readln;
1:      gamser := sum*exp(-x+a*ln(x)-gln)
   END
END;
