PROCEDURE qcksrt(n: integer; VAR arr: glarray);
(* Programs using routine QCKSRT must define the type
TYPE
   glarray = ARRAY [1..np] OF real;
in the main routine, with np >= n.   *)
LABEL 11,21,22,30,99;
CONST
   m=7;
   nstack=50;
   fm=7875;
   fa=211.0;
   fc=1663.0;
VAR
   l,jstack,j,ir,iq,i: integer;
   fx,fmi,a: real;
   istack: ARRAY[1..nstack] OF integer;
BEGIN
   fmi := 1.0/fm;
   jstack := 0;
   l := 1;
   ir := n;
   fx := 0.0;
   WHILE true DO BEGIN
      IF ((ir-l) < m) THEN BEGIN
         FOR j := l+1 to ir DO BEGIN
            a := arr[j];
            FOR i := j-1 DOWNTO 1 DO BEGIN
               IF (arr[i] <= a) THEN GOTO 11;
               arr[i+1] := arr[i]
            END;
            i := 0;
11:            arr[i+1] := a
         END;
         IF (jstack = 0) THEN GOTO 99;
         ir := istack[jstack];
         l := istack[jstack-1];
         jstack := jstack-2
      END ELSE BEGIN
         i := l;
         j := ir;
         fx := (fx*fa+fc)/fm;
         fx := fx-trunc(fx);
         iq := l+(ir-l+1)*trunc(fx*fmi);
         a := arr[iq];
         arr[iq] := arr[l];
21:         IF (j > 0) THEN BEGIN
            IF (a < arr[j]) THEN BEGIN
               j := j-1;
               GOTO 21
            END
         END;
         IF (j <= i) THEN BEGIN
            arr[i] := a;
            GOTO 30
         END;
         arr[i] := arr[j];
         i := i+1;
22:         IF (i <= n) THEN IF (a > arr[i]) THEN BEGIN
            i := i+1;
            GOTO 22
         END;
         IF (j <= i) THEN BEGIN
            arr[j] := a;
            i := j;
            GOTO 30
         END;
         arr[j] := arr[i];
         j := j-1;
         GOTO 21;
30:         jstack := jstack+2;
         IF (jstack > nstack) THEN BEGIN
            writeln('pause in QCKSRT - NSTACK must be made larger'); readln
         END;
         IF ((ir-i) >= (i-l)) THEN BEGIN
            istack[jstack] := ir;
            istack[jstack-1] := i+1;
            ir := i-1
         END ELSE BEGIN
            istack[jstack] := i-1;
            istack[jstack-1] := l;
            l := i+1
         END
      END
   END;
99:   END;
PROCEDURE qgaus(a,b: real; VAR ss: real);
(* Programs using routine QGAUS must externally define
a function func(x:real):real which is to be integrated *)
VAR
   j: integer;
   xr,xm,dx: real;
   w,x: ARRAY[1..5] OF real;
BEGIN
   x[1] := 0.1488743389;
   x[2] := 0.4333953941;
   x[3] := 0.6794095682;
   x[4] := 0.8650633666;
   x[5] := 0.97390652;
   w[1] := 0.2955242247;
   w[2] := 0.2692667193;
   w[3] := 0.2190863625;
   w[4] := 0.1494513491;
   w[5] := 0.06667134;
   xm := 0.5*(b+a);
   xr := 0.5*(b-a);
   ss := 0;
   FOR j := 1 to 5 DO BEGIN
      dx := xr*x[j];
      ss := ss+w[j]*(func(xm+dx)+func(xm-dx))
   END;
   ss := xr*ss
END;
PROCEDURE qromb(a,b: real; VAR ss: real);
(* Programs using routine QROMB must define type
TYPE
   glnarray = ARRAY [1..n] OF real;
just as for routine POLINT which it calls. In this case
n should have a value no smaller than the constant k below. *)
LABEL 99;
CONST
   eps=1.0e-6;
   jmax=20;
   jmaxp=21;   (* jmax+1 *)
   k=5;
VAR
   i,j: integer;
   dss: real;
   h,s: ARRAY[1..jmaxp] OF real;
   c,d: glnarray;
BEGIN
   h[1] := 1.0;
   FOR j := 1 to jmax DO BEGIN
      trapzd(a,b,s[j],j);
      IF  (j >= k)  THEN BEGIN
         FOR i := 1 to k DO BEGIN
            c[i] := h[j-k+i];
            d[i] := s[j-k+i]
         END;
         polint(c,d,k,0.0,ss,dss);
         IF (abs(dss) < eps*abs(ss)) THEN GOTO 99
      END;
      s[j+1] := s[j];
      h[j+1] := 0.25*h[j]
   END;
   writeln('pause in QROMB - too many steps'); readln;
99:   END;
PROCEDURE qromo(a,b: real; VAR ss: real);
(* Programs using routine QROMO must define the type
TYPE
   glnarray = ARRAY [1..n] OF real;
in the main routine, where n is equal to the constant k below. The routine
func(x:real):real in the calling routine must return the value of the function
to be integrated.  You must choose MIDPNT, MIDSQL, MIDSQU or MIDINF at the
indicated point below *)
LABEL 99;
CONST
   eps=1.0e-6;
   jmax=14;
   jmaxp=15;   (* jmaxp=jmax+1 *)
   k=5;
   km=4;      (* km=k-1 *)
VAR
   i,j: integer;
   dss: real;
   h,s: ARRAY [1..jmaxp] OF real;
   c,d: glnarray;
BEGIN
   h[1] := 1.0;
   FOR j := 1 to jmax DO BEGIN
(* Here you must choose the appropriate integration method *)
      midsql(a,b,s[j],j);
      IF  (j >= k)  THEN BEGIN
         FOR i := 1 to k DO BEGIN
            c[i] := h[j-k+i];
            d[i] := s[j-k+i]
         END;
         polint(c,d,k,0.0,ss,dss);
         IF (abs(dss) < (eps*abs(ss))) THEN GOTO 99
      END;
      s[j+1] := s[j];
      h[j+1] := h[j]/9.0
   END;
   writeln('pause in QROMO - too many steps');
   readln;
99: END;
PROCEDURE qroot(p: glnarray; n: integer; VAR b,c: real; eps: real);
(* Programs using procedure QROOT must define the types
TYPE
   glnarray = ARRAY [1..n] OF real;
   glnvarray = ARRAY [1..3] OF real;
in the main routine *)
LABEL 99;
CONST
   itmax=20;
   tiny=1.0e-6;
VAR
   iter,i: integer;
   sc,sb,s,rc,rb,r,epsf,dv,delc,delb: real;
   q,qq,rem: glnarray;
   d: glnvarray;
BEGIN
   d[3] := 1.0;
   FOR iter := 1 to itmax DO BEGIN
      d[2] := b;
      d[1] := c;
      poldiv(p,n,d,3,q,rem);
      s := rem[1];
      r := rem[2];
      poldiv(q,n-1,d,3,qq,rem);
      sc := -rem[1];
      rc := -rem[2];
      FOR i := n-1 DOWNTO 1 DO BEGIN
         q[i+1] := q[i]
      END;
      q[1] := 0.0;
      poldiv(q,n,d,3,qq,rem);
      sb := -rem[1];
      rb := -rem[2];
      dv := 1.0/(sb*rc-sc*rb);
      delb := (r*sc-s*rc)*dv;
      delc := (-r*sb+s*rb)*dv;
      b := b+delb;
      c := c+delc;
      IF(((abs(delb) <= eps*abs(b)) OR (abs(b) < tiny)) AND
      ((abs(delc) <= eps*abs(c)) OR (abs(c) < tiny))) THEN GOTO 99
   END;
   writeln('pause in routine QROOT - too many iterations');
99:   END;
PROCEDURE qsimp(a,b: real; VAR s: real);
LABEL 99;
CONST
   eps=1.0e-6;
   jmax=20;
VAR
   j: integer;
   st,ost,os: real;
BEGIN
   ost := -1.0e30;
   os := -1.0e30;
   FOR j := 1 to jmax DO BEGIN
      trapzd(a,b,st,j);
      s := (4.0*st-ost)/3.0;
      IF (abs(s-os) < eps*abs(os)) THEN GOTO 99;
      os := s;
      ost := st
   END;
   writeln ('pause in QSIMP - too many steps'); readln;
99:   END;
PROCEDURE qtrap(a,b: real;VAR s: real);
LABEL 99;
CONST
   eps=1.0e-6;
   jmax=20;
VAR
   j: integer;
   olds: real;
BEGIN
   olds := -1.0e30;
   FOR j := 1 to jmax DO BEGIN
      trapzd(a,b,s,j);
      IF  (abs(s-olds) < eps*abs(olds)) THEN GOTO 99;
      olds := s
   END;
   writeln ('pause in QTRAP - too many steps'); readln;
99:   END;
PROCEDURE quad3d(x1,x2: real; VAR ss: real);
(* Evaluates 3-dimensional integral with z integration innermost, then
y integration, and finally x integration. Unlike FORTRAN version, calls QGAUS
(here called QGAUS3) recursively. Programs using routine QUAD3D must define
the integrand by
FUNCTION func(x,y,z: real): real;
and functions for the limits of integration by
FUNCTION y1(x: real): real;
FUNCTION y2(x: real): real;
FUNCTION z1(x,y: real): real;
FUNCTION z2(x,y: real): real;
Also global variables
VAR
   glx,gly: real;
are required. *)
   PROCEDURE qgaus3(a,b: real; VAR ss: real; n: integer);
   VAR
      j: integer;
      xr,xm,dx: real;
      w,x: ARRAY [1..5] OF real;
   FUNCTION f(x: real; n: integer): real;
      VAR
         ss: real;
      BEGIN
         IF n=1 THEN BEGIN
            glx := x;
            qgaus3(y1(glx),y2(glx),ss,2);
            f := ss
         END ELSE IF n=2 THEN BEGIN
            gly := x;
            qgaus3(z1(glx,gly),z2(glx,gly),ss,3);
            f := ss
         END ELSE f := func(glx,gly,x)
      END;
   BEGIN
      x[1] := 0.1488743389;
      x[2] := 0.4333953941;
      x[3] := 0.6794095682;
      x[4] := 0.8650633666;
      x[5] := 0.97390652;
      w[1] := 0.2955242247;
      w[2] := 0.2692667193;
      w[3] := 0.2190863625;
      w[4] := 0.1494513491;
      w[5] := 0.06667134;
      xm := 0.5*(b+a);
      xr := 0.5*(b-a);
      ss := 0;
      FOR j := 1 TO 5 DO BEGIN
         dx := xr*x[j];
         ss := ss+w[j]*(f(xm+dx,n)+f(xm-dx,n))
      END;
      ss := xr*ss
   END;
BEGIN
   qgaus3(x1,x2,ss,1)
END;
