      SUBROUTINE MDIAN1(X,N,XMED)
      DIMENSION X(N)
      CALL SORT(N,X)
      N2=N/2
      IF(2*N2.EQ.N)THEN
        XMED=0.5*(X(N2)+X(N2+1))
      ELSE
        XMED=X(N2+1)
      ENDIF
      RETURN
      END
      SUBROUTINE MDIAN2(X,N,XMED)
      DIMENSION X(N)
      PARAMETER (BIG=1.E30,AFAC=1.5,AMP=1.5)
      A=0.5*(X(1)+X(N))
      EPS=ABS(X(N)-X(1))
      AP=BIG
      AM=-BIG
1     SUM=0.
      SUMX=0.
      NP=0
      NM=0
      XP=BIG
      XM=-BIG
      DO 11 J=1,N
        XX=X(J)
        IF(XX.NE.A)THEN
          IF(XX.GT.A)THEN
            NP=NP+1
            IF(XX.LT.XP)XP=XX
          ELSE IF(XX.LT.A)THEN
            NM=NM+1
            IF(XX.GT.XM)XM=XX
          ENDIF
          DUM=1./(EPS+ABS(XX-A))
          SUM=SUM+DUM
          SUMX=SUMX+XX*DUM
        ENDIF
11    CONTINUE
      IF(NP-NM.GE.2)THEN
        AM=A
        AA=XP+MAX(0.,SUMX/SUM-A)*AMP
        IF(AA.GT.AP)AA=0.5*(A+AP)
        EPS=AFAC*ABS(AA-A)
        A=AA
        GO TO 1
      ELSE IF(NM-NP.GE.2)THEN
        AP=A
        AA=XM+MIN(0.,SUMX/SUM-A)*AMP
        IF(AA.LT.AM)AA=0.5*(A+AM)
        EPS=AFAC*ABS(AA-A)
        A=AA
        GO TO 1
      ELSE
        IF(MOD(N,2).EQ.0)THEN
          IF(NP.EQ.NM)THEN
            XMED=0.5*(XP+XM)
          ELSE IF(NP.GT.NM)THEN
            XMED=0.5*(A+XP)
          ELSE
            XMED=0.5*(XM+A)
          ENDIF
        ELSE
          IF(NP.EQ.NM)THEN
            XMED=A
          ELSE IF(NP.GT.NM)THEN
            XMED=XP
          ELSE
            XMED=XM
          ENDIF
        ENDIF
      ENDIF
      RETURN
      END
      SUBROUTINE MEDFIT(X,Y,NDATA,A,B,ABDEV)
      PARAMETER (NMAX=1000)
      EXTERNAL ROFUNC
      COMMON /ARRAYS/ NDATAT,XT(NMAX),YT(NMAX),ARR(NMAX),AA,ABDEVT
      DIMENSION X(NDATA),Y(NDATA)
      SX=0.
      SY=0.
      SXY=0.
      SXX=0.
      DO 11 J=1,NDATA
        XT(J)=X(J)
        YT(J)=Y(J)
        SX=SX+X(J)
        SY=SY+Y(J)
        SXY=SXY+X(J)*Y(J)
        SXX=SXX+X(J)**2
11    CONTINUE
      NDATAT=NDATA
      DEL=NDATA*SXX-SX**2
      AA=(SXX*SY-SX*SXY)/DEL
      BB=(NDATA*SXY-SX*SY)/DEL
      CHISQ=0.
      DO 12 J=1,NDATA
        CHISQ=CHISQ+(Y(J)-(AA+BB*X(J)))**2
12    CONTINUE
      SIGB=SQRT(CHISQ/DEL)
      B1=BB
      F1=ROFUNC(B1)
      B2=BB+SIGN(3.*SIGB,F1)
      F2=ROFUNC(B2)
1     IF(F1*F2.GT.0.)THEN
        BB=2.*B2-B1
        B1=B2
        F1=F2
        B2=BB
        F2=ROFUNC(B2)
        GOTO 1
      ENDIF
      SIGB=0.01*SIGB
2     IF(ABS(B2-B1).GT.SIGB)THEN
        BB=0.5*(B1+B2)
        IF(BB.EQ.B1.OR.BB.EQ.B2)GOTO 3
        F=ROFUNC(BB)
        IF(F*F1.GE.0.)THEN
          F1=F
          B1=BB
        ELSE
          F2=F
          B2=BB
        ENDIF
        GOTO 2
      ENDIF
3     A=AA
      B=BB
      ABDEV=ABDEVT/NDATA
      RETURN
      END
      SUBROUTINE MEMCOF(DATA,N,M,PM,COF,WK1,WK2,WKM)
      DIMENSION DATA(N),COF(M),WK1(N),WK2(N),WKM(M)
      P=0.
      DO 11 J=1,N
        P=P+DATA(J)**2
11    CONTINUE
      PM=P/N
      WK1(1)=DATA(1)
      WK2(N-1)=DATA(N)
      DO 12 J=2,N-1
        WK1(J)=DATA(J)
        WK2(J-1)=DATA(J)
12    CONTINUE
      DO 17 K=1,M
        PNEUM=0.
        DENOM=0.
        DO 13 J=1,N-K
          PNEUM=PNEUM+WK1(J)*WK2(J)
          DENOM=DENOM+WK1(J)**2+WK2(J)**2
13      CONTINUE
        COF(K)=2.*PNEUM/DENOM
        PM=PM*(1.-COF(K)**2)
        IF(K.NE.1)THEN
          DO 14 I=1,K-1
            COF(I)=WKM(I)-COF(K)*WKM(K-I)
14        CONTINUE
        ENDIF
        IF(K.EQ.M)RETURN
        DO 15 I=1,K
          WKM(I)=COF(I)
15      CONTINUE
        DO 16 J=1,N-K-1
          WK1(J)=WK1(J)-WKM(K)*WK2(J)
          WK2(J)=WK2(J+1)-WKM(K)*WK1(J+1)
16      CONTINUE
17    CONTINUE
      PAUSE 'never get here'
      END
      SUBROUTINE MIDINF(FUNK,AA,BB,S,N)
      FUNC(X)=FUNK(1./X)/X**2
      B=1./AA
      A=1./BB
      IF (N.EQ.1) THEN
        S=(B-A)*FUNC(0.5*(A+B))
        IT=1
      ELSE
        TNM=IT
        DEL=(B-A)/(3.*TNM)
        DDEL=DEL+DEL
        X=A+0.5*DEL
        SUM=0.
        DO 11 J=1,IT
          SUM=SUM+FUNC(X)
          X=X+DDEL
          SUM=SUM+FUNC(X)
          X=X+DEL
11      CONTINUE
        S=(S+(B-A)*SUM/TNM)/3.
        IT=3*IT
      ENDIF
      RETURN
      END
      SUBROUTINE MIDPNT(FUNC,A,B,S,N)
      IF (N.EQ.1) THEN
        S=(B-A)*FUNC(0.5*(A+B))
        IT=1
      ELSE
        TNM=IT
        DEL=(B-A)/(3.*TNM)
        DDEL=DEL+DEL
        X=A+0.5*DEL
        SUM=0.
        DO 11 J=1,IT
          SUM=SUM+FUNC(X)
          X=X+DDEL
          SUM=SUM+FUNC(X)
          X=X+DEL
11      CONTINUE
        S=(S+(B-A)*SUM/TNM)/3.
        IT=3*IT
      ENDIF
      RETURN
      END
      SUBROUTINE MMID(Y,DYDX,NVAR,XS,HTOT,NSTEP,YOUT,DERIVS)
      PARAMETER (NMAX=10)
      DIMENSION Y(NVAR),DYDX(NVAR),YOUT(NVAR),YM(NMAX),YN(NMAX)
      H=HTOT/NSTEP
      DO 11 I=1,NVAR
        YM(I)=Y(I)
        YN(I)=Y(I)+H*DYDX(I)
11    CONTINUE
      X=XS+H
      CALL DERIVS(X,YN,YOUT)
      H2=2.*H
      DO 13 N=2,NSTEP
        DO 12 I=1,NVAR
          SWAP=YM(I)+H2*YOUT(I)
          YM(I)=YN(I)
          YN(I)=SWAP
12      CONTINUE
        X=X+H
        CALL DERIVS(X,YN,YOUT)
13    CONTINUE
      DO 14 I=1,NVAR
        YOUT(I)=0.5*(YM(I)+YN(I)+H*YOUT(I))
14    CONTINUE
      RETURN
      END
      SUBROUTINE MNBRAK(AX,BX,CX,FA,FB,FC,FUNC)
      PARAMETER (GOLD=1.618034, GLIMIT=100., TINY=1.E-20)
      FA=FUNC(AX)
      FB=FUNC(BX)
      IF(FB.GT.FA)THEN
        DUM=AX
        AX=BX
        BX=DUM
        DUM=FB
        FB=FA
        FA=DUM
      ENDIF
      CX=BX+GOLD*(BX-AX)
      FC=FUNC(CX)
1     IF(FB.GE.FC)THEN
        R=(BX-AX)*(FB-FC)
        Q=(BX-CX)*(FB-FA)
        U=BX-((BX-CX)*Q-(BX-AX)*R)/(2.*SIGN(MAX(ABS(Q-R),TINY),Q-R))
        ULIM=BX+GLIMIT*(CX-BX)
        IF((BX-U)*(U-CX).GT.0.)THEN
          FU=FUNC(U)
          IF(FU.LT.FC)THEN
            AX=BX
            FA=FB
            BX=U
            FB=FU
            GO TO 1
          ELSE IF(FU.GT.FB)THEN
            CX=U
            FC=FU
            GO TO 1
          ENDIF
          U=CX+GOLD*(CX-BX)
          FU=FUNC(U)
        ELSE IF((CX-U)*(U-ULIM).GT.0.)THEN
          FU=FUNC(U)
          IF(FU.LT.FC)THEN
            BX=CX
            CX=U
            U=CX+GOLD*(CX-BX)
            FB=FC
            FC=FU
            FU=FUNC(U)
          ENDIF
        ELSE IF((U-ULIM)*(ULIM-CX).GE.0.)THEN
          U=ULIM
          FU=FUNC(U)
        ELSE
          U=CX+GOLD*(CX-BX)
          FU=FUNC(U)
        ENDIF
        AX=BX
        BX=CX
        CX=U
        FA=FB
        FB=FC
        FC=FU
        GO TO 1
      ENDIF
      RETURN
      END
      SUBROUTINE MNEWT(NTRIAL,X,N,TOLX,TOLF)
      PARAMETER (NP=15)
      DIMENSION X(NP),ALPHA(NP,NP),BETA(NP),INDX(NP)
      DO 13  K=1,NTRIAL
        CALL USRFUN(X,ALPHA,BETA)
        ERRF=0.
        DO 11 I=1,N
          ERRF=ERRF+ABS(BETA(I))
11      CONTINUE
        IF(ERRF.LE.TOLF)RETURN
        CALL LUDCMP(ALPHA,N,NP,INDX,D)
        CALL LUBKSB(ALPHA,N,NP,INDX,BETA)
        ERRX=0.
        DO 12 I=1,N
          ERRX=ERRX+ABS(BETA(I))
          X(I)=X(I)+BETA(I)
12      CONTINUE
        IF(ERRX.LE.TOLX)RETURN
13    CONTINUE
      RETURN
      END
      SUBROUTINE MOMENT(DATA,N,AVE,ADEV,SDEV,VAR,SKEW,CURT)
      DIMENSION DATA(N)
      IF(N.LE.1)PAUSE 'N must be at least 2'
      S=0.
      DO 11 J=1,N
        S=S+DATA(J)
11    CONTINUE
      AVE=S/N
      ADEV=0.
      VAR=0.
      SKEW=0.
      CURT=0.
      DO 12 J=1,N
        S=DATA(J)-AVE
        ADEV=ADEV+ABS(S)
        P=S*S
        VAR=VAR+P
        P=P*S
        SKEW=SKEW+P
        P=P*S
        CURT=CURT+P
12    CONTINUE
      ADEV=ADEV/N
      VAR=VAR/(N-1)
      SDEV=SQRT(VAR)
      IF(VAR.NE.0.)THEN
        SKEW=SKEW/(N*SDEV**3)
        CURT=CURT/(N*VAR**2)-3.  
      ELSE
        PAUSE 'no skew or kurtosis when zero variance'
      ENDIF
      RETURN
      END
      SUBROUTINE MPROVE(A,ALUD,N,NP,INDX,B,X)
      PARAMETER (NMAX=100)
      DIMENSION A(NP,NP),ALUD(NP,NP),INDX(N),B(N),X(N),R(NMAX)
      REAL*8 SDP
      DO 12 I=1,N
        SDP=-B(I)
        DO 11 J=1,N
          SDP=SDP+DBLE(A(I,J))*DBLE(X(J))
11      CONTINUE
        R(I)=SDP
12    CONTINUE
      CALL LUBKSB(ALUD,N,NP,INDX,R)
      DO 13 I=1,N
        X(I)=X(I)-R(I)
13    CONTINUE
      RETURN
      END
      SUBROUTINE MRQCOF(X,Y,SIG,NDATA,A,MA,LISTA,MFIT,ALPHA,BETA,NALP,
     *CHISQ,FUNCS)
      PARAMETER (MMAX=20)
      DIMENSION X(NDATA),Y(NDATA),SIG(NDATA),ALPHA(NALP,NALP),BETA(MA),
     *    DYDA(MMAX),LISTA(MFIT),A(MA)
      DO 12 J=1,MFIT
        DO 11 K=1,J
          ALPHA(J,K)=0.
11      CONTINUE
        BETA(J)=0.
12    CONTINUE
      CHISQ=0.
      DO 15 I=1,NDATA
        CALL FUNCS(X(I),A,YMOD,DYDA,MA)
        SIG2I=1./(SIG(I)*SIG(I))
        DY=Y(I)-YMOD
        DO 14 J=1,MFIT
          WT=DYDA(LISTA(J))*SIG2I
          DO 13 K=1,J
            ALPHA(J,K)=ALPHA(J,K)+WT*DYDA(LISTA(K))
13        CONTINUE
          BETA(J)=BETA(J)+DY*WT
14      CONTINUE
        CHISQ=CHISQ+DY*DY*SIG2I
15    CONTINUE
      DO 17 J=2,MFIT
        DO 16 K=1,J-1
          ALPHA(K,J)=ALPHA(J,K)
16      CONTINUE
17    CONTINUE
      RETURN
      END
      SUBROUTINE MRQMIN(X,Y,SIG,NDATA,A,MA,LISTA,MFIT,
     *    COVAR,ALPHA,NCA,CHISQ,FUNCS,ALAMDA)
      PARAMETER (MMAX=20)
      DIMENSION X(NDATA),Y(NDATA),SIG(NDATA),A(MA),LISTA(MA),
     *  COVAR(NCA,NCA),ALPHA(NCA,NCA),ATRY(MMAX),BETA(MMAX),DA(MMAX)
      IF(ALAMDA.LT.0.)THEN
        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 permutation in LISTA'
          ENDIF
12      CONTINUE
        IF (KK.NE.(MA+1)) PAUSE 'Improper permutation in LISTA'
        ALAMDA=0.001
        CALL MRQCOF(X,Y,SIG,NDATA,A,MA,LISTA,MFIT,ALPHA,BETA,NCA,CHISQ,F
     *UNCS)
        OCHISQ=CHISQ
        DO 13 J=1,MA
          ATRY(J)=A(J)
13      CONTINUE
      ENDIF
      DO 15 J=1,MFIT
        DO 14 K=1,MFIT
          COVAR(J,K)=ALPHA(J,K)
14      CONTINUE
        COVAR(J,J)=ALPHA(J,J)*(1.+ALAMDA)
        DA(J)=BETA(J)
15    CONTINUE
      CALL GAUSSJ(COVAR,MFIT,NCA,DA,1,1)
      IF(ALAMDA.EQ.0.)THEN
        CALL COVSRT(COVAR,NCA,MA,LISTA,MFIT)
        RETURN
      ENDIF
      DO 16 J=1,MFIT
        ATRY(LISTA(J))=A(LISTA(J))+DA(J)
16    CONTINUE
      CALL MRQCOF(X,Y,SIG,NDATA,ATRY,MA,LISTA,MFIT,COVAR,DA,NCA,CHISQ,FU
     *NCS)
      IF(CHISQ.LT.OCHISQ)THEN
        ALAMDA=0.1*ALAMDA
        OCHISQ=CHISQ
        DO 18 J=1,MFIT
          DO 17 K=1,MFIT
            ALPHA(J,K)=COVAR(J,K)
17        CONTINUE
          BETA(J)=DA(J)
          A(LISTA(J))=ATRY(LISTA(J))
18      CONTINUE
      ELSE
        ALAMDA=10.*ALAMDA
        CHISQ=OCHISQ
      ENDIF
      RETURN
      END
