PROCEDURE eclass(VAR nf: glnarray; n: integer;
       lista,listb: glmarray; m: integer);
(* Programs using routine ECLASS must define the types
TYPE
   glnarray = ARRAY [1..n] OF integer;
   glmarray = ARRAY [1..m] OF integer;
in the main routine. *)
VAR
   l,k,j: integer;
BEGIN
   FOR k := 1 to n DO BEGIN
      nf[k] := k
   END;
   FOR l := 1 to m DO BEGIN
      j := lista[l];
      WHILE (nf[j] <> j) DO j := nf[j];
      k := listb[l];
      WHILE (nf[k] <> k) DO k := nf[k];
      IF (j <> k) THEN nf[j] := k
   END;
   FOR j := 1 to n DO WHILE (nf[j] <> nf[nf[j]]) DO nf[j] := nf[nf[j]];
END;

PROCEDURE eclazz(VAR nf: glnarray; n: integer);
(* Programs using routine ECLAZZ must supply a boolean function
equiv(i,j:integer):boolean which indicates (TRUE or FALSE) whether
i and j belong to the same equivalence class. They must also
define the type
TYPE
   glnarray = ARRAY [1..n] OF integer;
in the main routine. *)
VAR
   kk,jj: integer;
BEGIN
   nf[1] := 1;
   FOR jj := 2 to n DO BEGIN
      nf[jj] := jj;
      FOR kk := 1 to jj-1 DO BEGIN
          nf[kk] := nf[nf[kk]];
          IF (equiv(jj,kk)) THEN nf[nf[nf[kk]]] := jj
      END
   END;
   FOR jj := 1 to n DO nf[jj] := nf[nf[jj]]
END;
PROCEDURE eigsrt(VAR d: glnp; VAR v: glnpnp; n: integer);
(* Programs using routine EIGSRT must define the types
TYPE
   glnp = ARRAY [1..np] OF real;
   glnpnp = ARRAY [1..np,1..np] OF real;
where np is the physical dimension of the arrays to be used (v and d) *)
VAR
   k,j,i: integer;
   p: real;
BEGIN
   FOR i := 1 to n-1 DO BEGIN
      k := i;
      p := d[i];
      FOR j := i+1 to n DO BEGIN
         IF (d[j] >= p) THEN BEGIN
            k := j;
            p := d[j]
         END
      END;
      IF (k <> i) THEN BEGIN
         d[k] := d[i];
         d[i] := p;
         FOR j := 1 to n DO BEGIN
            p := v[j,i];
            v[j,i] := v[j,k];
            v[j,k] := p
         END
      END
   END
END;

FUNCTION el2(x,qqc,aa,bb: real): real;
LABEL 1;
CONST
   pi=3.14159265;
   ca=0.0003;
   cb=1.0e-9;
VAR
   a,b,c,d,e,f,g: real;
   em,eye,p,qc,y,z: real;
   l: integer;
BEGIN
   IF (x = 0.0) THEN el2 := 0.0
   ELSE IF (qqc <> 0.0) THEN BEGIN
      qc := qqc;
      a := aa;
      b := bb;
      c := sqr(x);
      d := 1.0+c;
      p := sqrt((1.0+c*sqr(qc))/d);
      d := x/d;
      c := d/(2.0*p);
      z := a-b;
      eye := a;
      a := 0.5*(b+a);
      y := abs(1.0/x);
      f := 0.0;
      l := 0;
      em := 1.0;
      qc := abs(qc);
1:      b := eye*qc+b;
      e := em*qc;
      g := e/p;
      d := f*g+d;
      f := c;
      eye := a;
      p := g+p;
      c := 0.5*(d/p+c);
      g := em;
      em := qc+em;
      a := 0.5*(b/em+a);
      y := -e/y+y;
      IF (y = 0.0) THEN y := sqrt(e)*cb;
      IF (abs(g-qc) > ca*g) THEN BEGIN
         qc := sqrt(e)*2.0;
         l := l+l;
         IF (y < 0.0) THEN l := l+1;
         GOTO 1
      END;
      IF (y < 0.0) THEN l := l+1;
      e := (arctan(em/y)+pi*l)*a/em;
      IF (x < 0.0) THEN e := -e;
      el2 := e+c*z
   END ELSE BEGIN
      writeln('pause in routine EL2'); readln
   END
END;
PROCEDURE elmhes(VAR a: glnpnp; n: integer);
(* Programs using routine ELMHES must define the type
TYPE
   glnpnp = ARRAY [1..np,1..np]
where 'np by np' is the physical dimension of the matrix to be reduced. *)
VAR
   m,j,i: integer;
   y,x: real;
BEGIN
   IF (n > 2) THEN BEGIN
      FOR m := 2 to n-1 DO BEGIN
         x := 0.0;
         i := m;
         FOR j := m to n DO BEGIN
            IF (abs(a[j,m-1]) > abs(x)) THEN BEGIN
               x := a[j,m-1];
               i := j
            END
         END;
         IF (i <> m) THEN BEGIN
            FOR j := m-1 to n DO BEGIN
               y := a[i,j];
               a[i,j] := a[m,j];
               a[m,j] := y
            END;
            FOR j := 1 to n DO BEGIN
               y := a[j,i];
               a[j,i] := a[j,m];
               a[j,m] := y
            END
         END;
         IF (x <> 0.0) THEN BEGIN
            FOR i := m+1 to n DO BEGIN
               y := a[i,m-1];
               IF (y <> 0.0) THEN BEGIN
                  y := y/x;
                  a[i,m-1] := y;
                  FOR j := m to n DO a[i,j] := a[i,j]-y*a[m,j];
                  FOR j := 1 to n DO a[j,m] := a[j,m]+y*a[j,i]
               END
            END
         END
      END
   END
END;
FUNCTION erf(x: real): real;
BEGIN
   IF (x < 0.0) THEN BEGIN
      erf := -gammp(0.5,sqr(x))
   END ELSE BEGIN
      erf := gammp(0.5,sqr(x))
   END
END;
FUNCTION erfc(VAR x: real): real;
BEGIN
   IF (x < 0.0) THEN BEGIN
      erfc := 1.0+gammp(0.5,sqr(x))
   END ELSE BEGIN
      erfc := gammq(0.5,sqr(x))
   END
END;
FUNCTION erfcc(x: real): real;
VAR
   t,z,ans: real;
BEGIN
   z := abs(x);
   t := 1.0/(1.0+0.5*z);
   ans := t*exp(-z*z-1.26551223+t*(1.00002368+
      t*(0.37409196+t*(0.09678418+t*(-0.18628806+
      t*(0.27886807+t*(-1.13520398+t*(1.48851587+
      t*(-0.82215223+t*0.17087277)))))))));
   IF (x >= 0.0) THEN erfcc := ans
   ELSE erfcc := 2.0-ans
END;
PROCEDURE eulsum(VAR sum: real; term: real; jterm: integer);
(* Programs using routine EULSUM must declare the variable array
   glwksp: ARRAY [1..np] OF real;
where np is a physical dimension larger than any value of jterm
to be used. Also declare
VAR
   glnterm: integer; *)
VAR
   j: integer;
   tmp,dum: real;
BEGIN
   IF (jterm = 1) THEN BEGIN
      glnterm := 1;
      glwksp[1] := term;
      sum := 0.5*term
   END ELSE BEGIN
      tmp := glwksp[1];
      glwksp[1] := term;
      FOR j := 1 to glnterm DO BEGIN
         dum := glwksp[j+1];
         glwksp[j+1] := 0.5*(glwksp[j]+tmp);
         tmp := dum
      END;
      IF (abs(glwksp[glnterm+1]) <= abs(glwksp[glnterm])) THEN BEGIN
         sum := sum+0.5*glwksp[glnterm+1];
         glnterm := glnterm+1
      END ELSE BEGIN
         sum := sum+glwksp[glnterm+1]
      END
   END
END;
FUNCTION evlmem(fdt: real; cof: glmarray; m: integer; pm: real): real;
(* Programs using routine EVLMEM must define the types
TYPE
   glmarray = ARRAY [1..m] OF real;
where m is the dimension of the array of coefficients. *)
VAR
   wr,wi,wpr,wpi,wtemp,theta: double;
   sumi,sumr: real;
   i: integer;
BEGIN
   theta := 6.28318530717959*fdt;
   wpr := cos(theta);
   wpi := sin(theta);
   wr := 1.0;
   wi := 0.0;
   sumr := 1.0;
   sumi := 0.0;
   FOR i := 1 to m DO BEGIN
      wtemp := wr;
      wr := wr*wpr-wi*wpi;
      wi := wi*wpr+wtemp*wpi;
      sumr := sumr-cof[i]*sngl(wr);
      sumi := sumi-cof[i]*sngl(wi)
   END;
   evlmem := pm/(sqr(sumr)+sqr(sumi))
END;
FUNCTION expdev(VAR idum: integer): real;
BEGIN
   expdev := -ln(ran3(idum))
END;
