      SUBROUTINE LAGUER(A,M,X,EPS,POLISH)
      COMPLEX A(*),X,DX,X1,B,D,F,G,H,SQ,GP,GM,G2,ZERO
      LOGICAL POLISH
      PARAMETER (ZERO=(0.,0.),EPSS=6.E-8,MAXIT=100)
      DXOLD=CABS(X)
      DO 12 ITER=1,MAXIT
        B=A(M+1)
        ERR=CABS(B)
        D=ZERO
        F=ZERO
        ABX=CABS(X)
        DO 11 J=M,1,-1
          F=X*F+D
          D=X*D+B
          B=X*B+A(J)
          ERR=CABS(B)+ABX*ERR
11      CONTINUE
        ERR=EPSS*ERR
        IF(CABS(B).LE.ERR) THEN
          DX=ZERO
          RETURN
        ELSE
          G=D/B
          G2=G*G
          H=G2-2.*F/B
          SQ=CSQRT((M-1)*(M*H-G2))
          GP=G+SQ
          GM=G-SQ
          IF(CABS(GP).LT.CABS(GM)) GP=GM
          DX=M/GP
        ENDIF
        X1=X-DX
        IF(X.EQ.X1)RETURN
        X=X1
        CDX=CABS(DX)
        IF(ITER.GT.6.AND.CDX.GE.DXOLD)RETURN
        DXOLD=CDX
        IF(.NOT.POLISH)THEN
          IF(CABS(DX).LE.EPS*CABS(X))RETURN
        ENDIF
12    CONTINUE
      PAUSE 'too many iterations'
      RETURN
      END
      SUBROUTINE LFIT(X,Y,SIG,NDATA,A,MA,LISTA,MFIT,COVAR,NCVM,CHISQ,
     *FUNCS)
      PARAMETER (MMAX=50)
      DIMENSION X(NDATA),Y(NDATA),SIG(NDATA),A(MA),LISTA(MA),
     *    COVAR(NCVM,NCVM),BETA(MMAX),AFUNC(MMAX)
      KK=MFIT+1
      DO 12 J=1,MA
        IHIT=0
        DO 11 K=1,MFIT
          IF (LISTA(K).EQ.J) IHIT=IHIT+1
11      CONTINUE
        IF (IHIT.EQ.0) THEN
          LISTA(KK)=J
          KK=KK+1
        ELSE IF (IHIT.GT.1) THEN
          PAUSE 'Improper set in LISTA'
        ENDIF
12    CONTINUE
      IF (KK.NE.(MA+1)) PAUSE 'Improper set in LISTA'
      DO 14 J=1,MFIT
        DO 13 K=1,MFIT
          COVAR(J,K)=0.
13      CONTINUE
        BETA(J)=0.
14    CONTINUE
      DO 18 I=1,NDATA
        CALL FUNCS(X(I),AFUNC,MA)
        YM=Y(I)
        IF(MFIT.LT.MA) THEN
          DO 15 J=MFIT+1,MA
            YM=YM-A(LISTA(J))*AFUNC(LISTA(J))
15        CONTINUE
        ENDIF
        SIG2I=1./SIG(I)**2
        DO 17 J=1,MFIT
          WT=AFUNC(LISTA(J))*SIG2I
          DO 16 K=1,J
            COVAR(J,K)=COVAR(J,K)+WT*AFUNC(LISTA(K))
16        CONTINUE
          BETA(J)=BETA(J)+YM*WT
17      CONTINUE
18    CONTINUE
      IF (MFIT.GT.1) THEN
        DO 21 J=2,MFIT
          DO 19 K=1,J-1
            COVAR(K,J)=COVAR(J,K)
19        CONTINUE
21      CONTINUE
      ENDIF
      CALL GAUSSJ(COVAR,MFIT,NCVM,BETA,1,1)
      DO 22 J=1,MFIT
        A(LISTA(J))=BETA(J)
22    CONTINUE
      CHISQ=0.
      DO 24 I=1,NDATA
        CALL FUNCS(X(I),AFUNC,MA)
        SUM=0.
        DO 23 J=1,MA
          SUM=SUM+A(J)*AFUNC(J)
23      CONTINUE
        CHISQ=CHISQ+((Y(I)-SUM)/SIG(I))**2
24    CONTINUE
      CALL COVSRT(COVAR,NCVM,MA,LISTA,MFIT)
      RETURN
      END
      SUBROUTINE REVCST(X,Y,IORDER,NCITY,N,DE)
C Rev. 12/13/85
      DIMENSION X(NCITY),Y(NCITY),IORDER(NCITY),N(6),XX(4),YY(4)
      ALEN(X1,X2,Y1,Y2)=SQRT((X2-X1)**2+(Y2-Y1)**2)
      N(3)=1+MOD((N(1)+NCITY-2),NCITY)
      N(4)=1+MOD(N(2),NCITY)
      DO 11 J=1,4
          II=IORDER(N(J))
          XX(J)=X(II)
          YY(J)=Y(II)
11    CONTINUE
      DE=-ALEN(XX(1),XX(3),YY(1),YY(3))-ALEN(XX(2),XX(4),YY(2),YY(4))
     *     +ALEN(XX(1),XX(4),YY(1),YY(4))+ALEN(XX(2),XX(3),YY(2),YY(3))
      RETURN
      END

      SUBROUTINE REVERS(IORDER,NCITY,N)
      DIMENSION IORDER(NCITY),N(6)
      NN=(1+MOD(N(2)-N(1)+NCITY,NCITY))/2
      DO 11 J=1,NN
          K=1+MOD((N(1)+J-2),NCITY)
          L=1+MOD((N(2)-J+NCITY),NCITY)    
          ITMP=IORDER(K)
          IORDER(K)=IORDER(L)
          IORDER(L)=ITMP
11    CONTINUE
      RETURN
      END

      SUBROUTINE TRNCST(X,Y,IORDER,NCITY,N,DE)
      DIMENSION X(NCITY),Y(NCITY),IORDER(NCITY),N(6),XX(6),YY(6)
      ALEN(X1,X2,Y1,Y2)=SQRT((X2-X1)**2+(Y2-Y1)**2)
      N(4)=1+MOD(N(3),NCITY)
      N(5)=1+MOD((N(1)+NCITY-2),NCITY)
      N(6)=1+MOD(N(2),NCITY)
      DO 11 J=1,6
          II=IORDER(N(J))
          XX(J)=X(II)
          YY(J)=Y(II)
11    CONTINUE
      DE=-ALEN(XX(2),XX(6),YY(2),YY(6))-ALEN(XX(1),XX(5),YY(1),YY(5))
     *     -ALEN(XX(3),XX(4),YY(3),YY(4))+ALEN(XX(1),XX(3),YY(1),YY(3))
     *     +ALEN(XX(2),XX(4),YY(2),YY(4))+ALEN(XX(5),XX(6),YY(5),YY(6))
      RETURN
      END    

      SUBROUTINE TRNSPT(IORDER,NCITY,N)
      PARAMETER(MAXCIT=1000)
      DIMENSION IORDER(NCITY),JORDER(MAXCIT),N(6)
      M1=1+MOD((N(2)-N(1)+NCITY),NCITY)
      M2=1+MOD((N(5)-N(4)+NCITY),NCITY)
      M3=1+MOD((N(3)-N(6)+NCITY),NCITY)
      NN=1  
      DO 11 J=1,M1
          JJ=1+MOD((J+N(1)-2),NCITY)
          JORDER(NN)=IORDER(JJ)
          NN=NN+1
11    CONTINUE
      IF (M2.GT.0) THEN
          DO 12 J=1,M2
              JJ=1+MOD((J+N(4)-2),NCITY)
              JORDER(NN)=IORDER(JJ)
              NN=NN+1
12        CONTINUE
      ENDIF
      IF (M3.GT.0) THEN
          DO 13 J=1,M3
              JJ=1+MOD((J+N(6)-2),NCITY)
              JORDER(NN)=IORDER(JJ)
              NN=NN+1
13        CONTINUE
      ENDIF
      DO 14 J=1,NCITY
          IORDER(J)=JORDER(J)
14    CONTINUE
      RETURN
      END

      SUBROUTINE METROP(DE,T,ANS)
      PARAMETER(JDUM=1)
      LOGICAL ANS
      ANS=(DE.LT.0.0).OR.(RAN3(JDUM).LT.EXP(-DE/T))
      RETURN
      END
      SUBROUTINE LINMIN(P,XI,N,FRET)
      PARAMETER (NMAX=50,TOL=1.E-4)
      EXTERNAL F1DIM
      DIMENSION P(N),XI(N)
      COMMON /F1COM/ NCOM,PCOM(NMAX),XICOM(NMAX)
      NCOM=N
      DO 11 J=1,N
        PCOM(J)=P(J)
        XICOM(J)=XI(J)
11    CONTINUE
      AX=0.
      XX=1.
      BX=2.
      CALL MNBRAK(AX,XX,BX,FA,FX,FB,F1DIM)
      FRET=BRENT(AX,XX,BX,F1DIM,TOL,XMIN)
      DO 12 J=1,N
        XI(J)=XMIN*XI(J)
        P(J)=P(J)+XI(J)
12    CONTINUE
      RETURN
      END
      SUBROUTINE LOCATE(XX,N,X,J)
      DIMENSION XX(N)
      JL=0
      JU=N+1
10    IF(JU-JL.GT.1)THEN
        JM=(JU+JL)/2
        IF((XX(N).GT.XX(1)).EQV.(X.GT.XX(JM)))THEN
          JL=JM
        ELSE
          JU=JM
        ENDIF
      GO TO 10
      ENDIF
      J=JL
      RETURN
      END
      SUBROUTINE LUBKSB(A,N,NP,INDX,B)
      DIMENSION A(NP,NP),INDX(N),B(N)
      II=0
      DO 12 I=1,N
        LL=INDX(I)
        SUM=B(LL)
        B(LL)=B(I)
        IF (II.NE.0)THEN
          DO 11 J=II,I-1
            SUM=SUM-A(I,J)*B(J)
11        CONTINUE
        ELSE IF (SUM.NE.0.) THEN
          II=I
        ENDIF
        B(I)=SUM
12    CONTINUE
      DO 14 I=N,1,-1
        SUM=B(I)
        IF(I.LT.N)THEN
          DO 13 J=I+1,N
            SUM=SUM-A(I,J)*B(J)
13        CONTINUE
        ENDIF
        B(I)=SUM/A(I,I)
14    CONTINUE
      RETURN
      END
      SUBROUTINE LUDCMP(A,N,NP,INDX,D)
      PARAMETER (NMAX=100,TINY=1.0E-20)
      DIMENSION A(NP,NP),INDX(N),VV(NMAX)
      D=1.
      DO 12 I=1,N
        AAMAX=0.
        DO 11 J=1,N
          IF (ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J))
11      CONTINUE
        IF (AAMAX.EQ.0.) PAUSE 'Singular matrix.'
        VV(I)=1./AAMAX
12    CONTINUE
      DO 19 J=1,N
        IF (J.GT.1) THEN
          DO 14 I=1,J-1
            SUM=A(I,J)
            IF (I.GT.1)THEN
              DO 13 K=1,I-1
                SUM=SUM-A(I,K)*A(K,J)
13            CONTINUE
              A(I,J)=SUM
            ENDIF
14        CONTINUE
        ENDIF
        AAMAX=0.
        DO 16 I=J,N
          SUM=A(I,J)
          IF (J.GT.1)THEN
            DO 15 K=1,J-1
              SUM=SUM-A(I,K)*A(K,J)
15          CONTINUE
            A(I,J)=SUM
          ENDIF
          DUM=VV(I)*ABS(SUM)
          IF (DUM.GE.AAMAX) THEN
            IMAX=I
            AAMAX=DUM
          ENDIF
16      CONTINUE
        IF (J.NE.IMAX)THEN
          DO 17 K=1,N
            DUM=A(IMAX,K)
            A(IMAX,K)=A(J,K)
            A(J,K)=DUM
17        CONTINUE
          D=-D
          VV(IMAX)=VV(J)
        ENDIF
        INDX(J)=IMAX
        IF(J.NE.N)THEN
          IF(A(J,J).EQ.0.)A(J,J)=TINY
          DUM=1./A(J,J)
          DO 18 I=J+1,N
            A(I,J)=A(I,J)*DUM
18        CONTINUE
        ENDIF
19    CONTINUE
      IF(A(N,N).EQ.0.)A(N,N)=TINY
      RETURN
      END
