      SUBROUTINE ECLASS(NF,N,LISTA,LISTB,M)
      DIMENSION NF(N),LISTA(M),LISTB(M)
      DO 11 K=1,N
        NF(K)=K
11    CONTINUE
      DO 12 L=1,M
        J=LISTA(L)
1       IF(NF(J).NE.J)THEN
          J=NF(J)
        GOTO 1
        ENDIF
        K=LISTB(L)
2       IF(NF(K).NE.K)THEN
          K=NF(K)
        GOTO 2
        ENDIF
        IF(J.NE.K)NF(J)=K
12    CONTINUE
      DO 13 J=1,N
3       IF(NF(J).NE.NF(NF(J)))THEN
          NF(J)=NF(NF(J))
        GOTO 3
        ENDIF
13    CONTINUE
      RETURN
      END
      SUBROUTINE ECLAZZ(NF,N,EQUIV)
      LOGICAL EQUIV
      DIMENSION NF(N)
      NF(1)=1
      DO 12 JJ=2,N
        NF(JJ)=JJ
        DO 11 KK=1,JJ-1
          NF(KK)=NF(NF(KK))
          IF (EQUIV(JJ,KK)) NF(NF(NF(KK)))=JJ
11      CONTINUE
12    CONTINUE
      DO 13 JJ=1,N
        NF(JJ)=NF(NF(JJ))
13    CONTINUE
      RETURN
      END
      SUBROUTINE EIGSRT(D,V,N,NP)
      DIMENSION D(NP),V(NP,NP)
      DO 13 I=1,N-1
        K=I
        P=D(I)
        DO 11 J=I+1,N
          IF(D(J).GE.P)THEN
            K=J
            P=D(J)
          ENDIF
11      CONTINUE
        IF(K.NE.I)THEN
          D(K)=D(I)
          D(I)=P
          DO 12 J=1,N
            P=V(J,I)
            V(J,I)=V(J,K)
            V(J,K)=P
12        CONTINUE
        ENDIF
13    CONTINUE
      RETURN
      END
      FUNCTION EL2(X,QQC,AA,BB)
      PARAMETER(PI=3.14159265, CA=.0003, CB=1.E-9)
      IF(X.EQ.0.)THEN
        EL2=0.
      ELSE IF(QQC.NE.0.)THEN
        QC=QQC
        A=AA
        B=BB
        C=X**2
        D=1.+C
        P=SQRT((1.+QC**2*C)/D)
        D=X/D
        C=D/(2.*P)
        Z=A-B
        EYE=A
        A=0.5*(B+A)
        Y=ABS(1./X)
        F=0.
        L=0
        EM=1.
        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.EQ.0.)Y=SQRT(E)*CB
        IF(ABS(G-QC).GT.CA*G)THEN
          QC=SQRT(E)*2.
          L=L+L
          IF(Y.LT.0.)L=L+1
          GO TO 1
        ENDIF
        IF(Y.LT.0.)L=L+1
        E=(ATAN(EM/Y)+PI*L)*A/EM
        IF(X.LT.0.)E=-E
        EL2=E+C*Z
      ELSE
        PAUSE 'failure in EL2'
      ENDIF
      RETURN
      END
      SUBROUTINE ELMHES(A,N,NP)
      DIMENSION A(NP,NP)
      IF(N.GT.2)THEN
        DO 17 M=2,N-1
          X=0.
          I=M
          DO 11 J=M,N
            IF(ABS(A(J,M-1)).GT.ABS(X))THEN
              X=A(J,M-1)
              I=J
            ENDIF
11        CONTINUE
          IF(I.NE.M)THEN
            DO 12 J=M-1,N
              Y=A(I,J)
              A(I,J)=A(M,J)
              A(M,J)=Y
12          CONTINUE
            DO 13 J=1,N
              Y=A(J,I)
              A(J,I)=A(J,M)
              A(J,M)=Y
13          CONTINUE
          ENDIF
          IF(X.NE.0.)THEN
            DO 16 I=M+1,N
              Y=A(I,M-1)
              IF(Y.NE.0.)THEN
                Y=Y/X
                A(I,M-1)=Y
                DO 14 J=M,N
                  A(I,J)=A(I,J)-Y*A(M,J)
14              CONTINUE
                DO 15 J=1,N
                  A(J,M)=A(J,M)+Y*A(J,I)
15              CONTINUE
              ENDIF
16          CONTINUE
          ENDIF
17      CONTINUE
      ENDIF
      RETURN
      END
      FUNCTION ERF(X)
      IF(X.LT.0.)THEN
        ERF=-GAMMP(.5,X**2)
      ELSE
        ERF=GAMMP(.5,X**2)
      ENDIF
      RETURN
      END
      FUNCTION ERFC(X)
      IF(X.LT.0.)THEN
        ERFC=1.+GAMMP(.5,X**2)
      ELSE
        ERFC=GAMMQ(.5,X**2)
      ENDIF
      RETURN
      END
      FUNCTION ERFCC(X)
      Z=ABS(X)      
      T=1./(1.+0.5*Z)
      ERFCC=T*EXP(-Z*Z-1.26551223+T*(1.00002368+T*(.37409196+
     *    T*(.09678418+T*(-.18628806+T*(.27886807+T*(-1.13520398+
     *    T*(1.48851587+T*(-.82215223+T*.17087277)))))))))
      IF (X.LT.0.) ERFCC=2.-ERFCC
      RETURN
      END
      SUBROUTINE EULSUM(SUM,TERM,JTERM,WKSP)
      DIMENSION WKSP(JTERM)
      IF(JTERM.EQ.1)THEN
        NTERM=1
        WKSP(1)=TERM
        SUM=0.5*TERM
      ELSE
        TMP=WKSP(1)
        WKSP(1)=TERM
        DO 11 J=1,NTERM
          DUM=WKSP(J+1)
          WKSP(J+1)=0.5*(WKSP(J)+TMP)
          TMP=DUM
11      CONTINUE
        IF(ABS(WKSP(NTERM+1)).LE.ABS(WKSP(NTERM)))THEN
          SUM=SUM+0.5*WKSP(NTERM+1)
          NTERM=NTERM+1
        ELSE
          SUM=SUM+WKSP(NTERM+1)
        ENDIF
      ENDIF
      RETURN
      END
      FUNCTION EVLMEM(FDT,COF,M,PM)
      DIMENSION COF(M)
      REAL*8 WR,WI,WPR,WPI,WTEMP,THETA
      THETA=6.28318530717959D0*FDT
      WPR=DCOS(THETA)
      WPI=DSIN(THETA)
      WR=1.D0
      WI=0.D0
      SUMR=1.
      SUMI=0.
      DO 11 I=1,M
        WTEMP=WR
        WR=WR*WPR-WI*WPI
        WI=WI*WPR+WTEMP*WPI
        SUMR=SUMR-COF(I)*SNGL(WR)
        SUMI=SUMI-COF(I)*SNGL(WI)
11    CONTINUE
      EVLMEM=PM/(SUMR**2+SUMI**2)
      RETURN
      END
      FUNCTION EXPDEV(IDUM)
      EXPDEV=-LOG(RAN1(IDUM))
      RETURN
      END
