      FUNCTION RAN(ISEED)
      PARAMETER(IA=7141,IC=54773,IM=259200)
      ISEED=MOD(ISEED*IA+IC,IM)
      RAN=FLOAT(ISEED)/FLOAT(IM)
      RETURN
      END
      FUNCTION RAN0(IDUM)
      DIMENSION V(97)
      DATA IFF /0/
      IF(IDUM.LT.0.OR.IFF.EQ.0)THEN
        IFF=1
        ISEED=ABS(IDUM)
        IDUM=1
        DO 11 J=1,97
          DUM=RAN(ISEED)
11      CONTINUE
        DO 12 J=1,97
          V(J)=RAN(ISEED)
12      CONTINUE
        Y=RAN(ISEED)
      ENDIF
      J=1+INT(97.*Y)
      IF(J.GT.97.OR.J.LT.1)PAUSE
      Y=V(J)
      RAN0=Y
      V(J)=RAN(ISEED)
      RETURN
      END
      FUNCTION RAN1(IDUM)
      DIMENSION R(97)
      PARAMETER (M1=259200,IA1=7141,IC1=54773,RM1=3.8580247E-6)
      PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=7.4373773E-6)
      PARAMETER (M3=243000,IA3=4561,IC3=51349)
      DATA IFF /0/
      IF (IDUM.LT.0.OR.IFF.EQ.0) THEN
        IFF=1
        IX1=MOD(IC1-IDUM,M1)
        IX1=MOD(IA1*IX1+IC1,M1)
        IX2=MOD(IX1,M2)
        IX1=MOD(IA1*IX1+IC1,M1)
        IX3=MOD(IX1,M3)
        DO 11 J=1,97
          IX1=MOD(IA1*IX1+IC1,M1)
          IX2=MOD(IA2*IX2+IC2,M2)
          R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
11      CONTINUE
        IDUM=1
      ENDIF
      IX1=MOD(IA1*IX1+IC1,M1)
      IX2=MOD(IA2*IX2+IC2,M2)
      IX3=MOD(IA3*IX3+IC3,M3)
      J=1+(97*IX3)/M3
      IF(J.GT.97.OR.J.LT.1)PAUSE
      RAN1=R(J)
      R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
      RETURN
      END
      FUNCTION RAN2(IDUM)
      PARAMETER (M=714025,IA=1366,IC=150889,RM=1.4005112E-6)
      DIMENSION IR(97)
      DATA IFF /0/
      IF(IDUM.LT.0.OR.IFF.EQ.0)THEN
        IFF=1
        IDUM=MOD(IC-IDUM,M)
        DO 11 J=1,97
          IDUM=MOD(IA*IDUM+IC,M)
          IR(J)=IDUM
11      CONTINUE
        IDUM=MOD(IA*IDUM+IC,M)
        IY=IDUM
      ENDIF
      J=1+(97*IY)/M
      IF(J.GT.97.OR.J.LT.1)PAUSE
      IY=IR(J)
      RAN2=IY*RM
      IDUM=MOD(IA*IDUM+IC,M)
      IR(J)=IDUM
      RETURN
      END
      FUNCTION RAN3(IDUM)
C         IMPLICIT REAL*4(M)
C         PARAMETER (MBIG=4000000.,MSEED=1618033.,MZ=0.,FAC=2.5E-7)
      PARAMETER (MBIG=1000000000,MSEED=161803398,MZ=0,FAC=1.E-9)
      DIMENSION MA(55)
      DATA IFF /0/
      IF(IDUM.LT.0.OR.IFF.EQ.0)THEN
        IFF=1
        MJ=MSEED-IABS(IDUM)
        MJ=MOD(MJ,MBIG)
        MA(55)=MJ
        MK=1
        DO 11 I=1,54
          II=MOD(21*I,55)
          MA(II)=MK
          MK=MJ-MK
          IF(MK.LT.MZ)MK=MK+MBIG
          MJ=MA(II)
11      CONTINUE
        DO 13 K=1,4
          DO 12 I=1,55
            MA(I)=MA(I)-MA(1+MOD(I+30,55))
            IF(MA(I).LT.MZ)MA(I)=MA(I)+MBIG
12        CONTINUE
13      CONTINUE
        INEXT=0
        INEXTP=31
        IDUM=1
      ENDIF
      INEXT=INEXT+1
      IF(INEXT.EQ.56)INEXT=1
      INEXTP=INEXTP+1
      IF(INEXTP.EQ.56)INEXTP=1
      MJ=MA(INEXT)-MA(INEXTP)
      IF(MJ.LT.MZ)MJ=MJ+MBIG
      MA(INEXT)=MJ
      RAN3=MJ*FAC
      RETURN
      END
      FUNCTION RAN4(IDUM)
      PARAMETER (IM=11979,IA=430,IC=2531,NACC=24)
      DIMENSION INP(64),JOT(64),KEY(64),POW(65)
      DATA IFF/0/
      IF(IDUM.LT.0.OR.IFF.EQ.0)THEN
        IFF=1
        IDUM=MOD(IDUM,IM)
        POW(1)=0.5
        DO 11 J=1,64
          IDUM=MOD(IDUM*IA+IC,IM)
          KEY(J)=(2*IDUM)/IM
          INP(J)=MOD((4*IDUM)/IM,2)
          POW(J+1)=0.5*POW(J)
11      CONTINUE
        NEWKEY=1
      ENDIF
      ISAV=INP(64)
      IF(ISAV.NE.0)THEN
        INP(4)=1-INP(4)
        INP(3)=1-INP(3)
        INP(1)=1-INP(1)
      ENDIF
      DO 12 J=64,2,-1
        INP(J)=INP(J-1)
12    CONTINUE
      INP(1)=ISAV
      CALL DES(INP,KEY,NEWKEY,0,JOT)
      RAN4=0.0
      DO 13 J=1,NACC
        IF(JOT(J).NE.0)RAN4=RAN4+POW(J)
13    CONTINUE
      RETURN
      END
      SUBROUTINE RANK(N,INDX,IRANK)
      DIMENSION INDX(N),IRANK(N)
      DO 11 J=1,N
        IRANK(INDX(J))=J
11    CONTINUE
      RETURN
      END
      SUBROUTINE RATINT(XA,YA,N,X,Y,DY)
      PARAMETER (NMAX=10,TINY=1.E-25)
      DIMENSION XA(N),YA(N),C(NMAX),D(NMAX)
      NS=1
      HH=ABS(X-XA(1))
      DO 11 I=1,N
        H=ABS(X-XA(I))
        IF (H.EQ.0.)THEN
          Y=YA(I)
          DY=0.0
          RETURN
        ELSE IF (H.LT.HH) THEN
          NS=I
          HH=H
        ENDIF
        C(I)=YA(I)
        D(I)=YA(I)+TINY
11    CONTINUE
      Y=YA(NS)
      NS=NS-1
      DO 13 M=1,N-1
        DO 12 I=1,N-M
          W=C(I+1)-D(I)
          H=XA(I+M)-X
          T=(XA(I)-X)*D(I)/H
          DD=T-C(I+1)
          IF(DD.EQ.0.)PAUSE
          DD=W/DD
          D(I)=C(I+1)*DD
          C(I)=T*DD
12      CONTINUE
        IF (2*NS.LT.N-M)THEN
          DY=C(NS+1)
        ELSE
          DY=D(NS)
          NS=NS-1
        ENDIF
        Y=Y+DY
13    CONTINUE
      RETURN
      END
      SUBROUTINE REALFT(DATA,N,ISIGN)
      REAL*8 WR,WI,WPR,WPI,WTEMP,THETA
      DIMENSION DATA(*)
      THETA=6.28318530717959D0/2.0D0/DBLE(N)
      C1=0.5
      IF (ISIGN.EQ.1) THEN
        C2=-0.5
        CALL FOUR1(DATA,N,+1)
      ELSE
        C2=0.5
        THETA=-THETA
      ENDIF
      WPR=-2.0D0*DSIN(0.5D0*THETA)**2
      WPI=DSIN(THETA)
      WR=1.0D0+WPR
      WI=WPI
      N2P3=2*N+3
      DO 11 I=2,N/2+1
        I1=2*I-1
        I2=I1+1
        I3=N2P3-I2
        I4=I3+1
        WRS=SNGL(WR)
        WIS=SNGL(WI)
        H1R=C1*(DATA(I1)+DATA(I3))
        H1I=C1*(DATA(I2)-DATA(I4))
        H2R=-C2*(DATA(I2)+DATA(I4))
        H2I=C2*(DATA(I1)-DATA(I3))
        DATA(I1)=H1R+WRS*H2R-WIS*H2I
        DATA(I2)=H1I+WRS*H2I+WIS*H2R
        DATA(I3)=H1R-WRS*H2R+WIS*H2I
        DATA(I4)=-H1I+WRS*H2I+WIS*H2R
        WTEMP=WR
        WR=WR*WPR-WI*WPI+WR
        WI=WI*WPR+WTEMP*WPI+WI
11    CONTINUE
      IF (ISIGN.EQ.1) THEN
        H1R=DATA(1)
        DATA(1)=H1R+DATA(2)
        DATA(2)=H1R-DATA(2)
      ELSE
        H1R=DATA(1)
        DATA(1)=C1*(H1R+DATA(2))
        DATA(2)=C1*(H1R-DATA(2))
        CALL FOUR1(DATA,N,-1)
      ENDIF
      RETURN
      END
      SUBROUTINE RED(IZ1,IZ2,JZ1,JZ2,JM1,JM2,JMF,IC1,JC1,JCF,KC,
     *    C,NCI,NCJ,NCK,S,NSI,NSJ)
      DIMENSION C(NCI,NCJ,NCK),S(NSI,NSJ)
      LOFF=JC1-JM1
      IC=IC1
      DO 14 J=JZ1,JZ2
        DO 12 L=JM1,JM2
          VX=C(IC,L+LOFF,KC)
          DO 11 I=IZ1,IZ2
            S(I,L)=S(I,L)-S(I,J)*VX
11        CONTINUE
12      CONTINUE
        VX=C(IC,JCF,KC)
        DO 13 I=IZ1,IZ2
          S(I,JMF)=S(I,JMF)-S(I,J)*VX
13      CONTINUE
        IC=IC+1
14    CONTINUE
      RETURN
      END
      SUBROUTINE RK4(Y,DYDX,N,X,H,YOUT,DERIVS)
      PARAMETER (NMAX=10)
      DIMENSION Y(N),DYDX(N),YOUT(N),YT(NMAX),DYT(NMAX),DYM(NMAX)
      HH=H*0.5
      H6=H/6.
      XH=X+HH
      DO 11 I=1,N
        YT(I)=Y(I)+HH*DYDX(I)
11    CONTINUE
      CALL DERIVS(XH,YT,DYT)
      DO 12 I=1,N
        YT(I)=Y(I)+HH*DYT(I)
12    CONTINUE
      CALL DERIVS(XH,YT,DYM)
      DO 13 I=1,N
        YT(I)=Y(I)+H*DYM(I)
        DYM(I)=DYT(I)+DYM(I)
13    CONTINUE
      CALL DERIVS(X+H,YT,DYT)
      DO 14 I=1,N
        YOUT(I)=Y(I)+H6*(DYDX(I)+DYT(I)+2.*DYM(I))
14    CONTINUE
      RETURN
      END
      SUBROUTINE RKDUMB(VSTART,NVAR,X1,X2,NSTEP,DERIVS)
      PARAMETER (NMAX=10)
      COMMON /PATH/ XX(200),Y(10,200)
      DIMENSION VSTART(NVAR),V(NMAX),DV(NMAX)
      DO 11 I=1,NVAR
        V(I)=VSTART(I)
        Y(I,1)=V(I)
11    CONTINUE
      XX(1)=X1
      X=X1
      H=(X2-X1)/NSTEP
      DO 13 K=1,NSTEP
        CALL DERIVS(X,V,DV)
        CALL RK4(V,DV,NVAR,X,H,V,DERIVS)
        IF(X+H.EQ.X)PAUSE 'Stepsize not significant in RKDUMB.'
        X=X+H
        XX(K+1)=X
        DO 12 I=1,NVAR
          Y(I,K+1)=V(I)
12      CONTINUE
13    CONTINUE
      RETURN
      END
      SUBROUTINE RKQC(Y,DYDX,N,X,HTRY,EPS,YSCAL,HDID,HNEXT,DERIVS)
      PARAMETER (NMAX=10,FCOR=.0666666667,
     *    ONE=1.,SAFETY=0.9,ERRCON=6.E-4)
      EXTERNAL DERIVS
      DIMENSION Y(N),DYDX(N),YSCAL(N),YTEMP(NMAX),YSAV(NMAX),DYSAV(NMAX)
      PGROW=-0.20
      PSHRNK=-0.25
      XSAV=X
      DO 11 I=1,N
        YSAV(I)=Y(I)
        DYSAV(I)=DYDX(I)
11    CONTINUE
      H=HTRY
1     HH=0.5*H
      CALL RK4(YSAV,DYSAV,N,XSAV,HH,YTEMP,DERIVS)
      X=XSAV+HH
      CALL DERIVS(X,YTEMP,DYDX)
      CALL RK4(YTEMP,DYDX,N,X,HH,Y,DERIVS)
      X=XSAV+H
      IF(X.EQ.XSAV)PAUSE 'Stepsize not significant in RKQC.'
      CALL RK4(YSAV,DYSAV,N,XSAV,H,YTEMP,DERIVS)
      ERRMAX=0.
      DO 12 I=1,N
        YTEMP(I)=Y(I)-YTEMP(I)
        ERRMAX=MAX(ERRMAX,ABS(YTEMP(I)/YSCAL(I)))
12    CONTINUE
      ERRMAX=ERRMAX/EPS
      IF(ERRMAX.GT.ONE) THEN
        H=SAFETY*H*(ERRMAX**PSHRNK)
        GOTO 1
      ELSE
        HDID=H
        IF(ERRMAX.GT.ERRCON)THEN
          HNEXT=SAFETY*H*(ERRMAX**PGROW)
        ELSE
          HNEXT=4.*H
        ENDIF
      ENDIF
      DO 13 I=1,N
        Y(I)=Y(I)+YTEMP(I)*FCOR
13    CONTINUE
      RETURN
      END
      FUNCTION ROFUNC(B)
      PARAMETER (NMAX=1000)
      COMMON /ARRAYS/ NDATA,X(NMAX),Y(NMAX),ARR(NMAX),AA,ABDEV
      N1=NDATA+1
      NML=N1/2
      NMH=N1-NML
      DO 11 J=1,NDATA
        ARR(J)=Y(J)-B*X(J)
11    CONTINUE
      CALL SORT(NDATA,ARR)
      AA=0.5*(ARR(NML)+ARR(NMH))
      SUM=0.
      ABDEV=0.
      DO 12 J=1,NDATA
        D=Y(J)-(B*X(J)+AA)
        ABDEV=ABDEV+ABS(D)
        SUM=SUM+X(J)*SIGN(1.0,D)
12    CONTINUE
      ROFUNC=SUM
      RETURN
      END
      FUNCTION RTBIS(FUNC,X1,X2,XACC)
      PARAMETER (JMAX=40)
      FMID=FUNC(X2)
      F=FUNC(X1)
      IF(F*FMID.GE.0.) PAUSE 'Root must be bracketed for bisection.'
      IF(F.LT.0.)THEN
        RTBIS=X1
        DX=X2-X1
      ELSE
        RTBIS=X2
        DX=X1-X2
      ENDIF
      DO 11 J=1,JMAX
        DX=DX*.5
        XMID=RTBIS+DX
        FMID=FUNC(XMID)
        IF(FMID.LE.0.)RTBIS=XMID
        IF(ABS(DX).LT.XACC .OR. FMID.EQ.0.) RETURN
11    CONTINUE
      PAUSE 'too many bisections'
      END
      FUNCTION RTFLSP(FUNC,X1,X2,XACC)
      PARAMETER (MAXIT=30)
      FL=FUNC(X1)
      FH=FUNC(X2)
      IF(FL*FH.GT.0.) PAUSE 'Root must be bracketed for false position.'
      IF(FL.LT.0.)THEN
        XL=X1
        XH=X2
      ELSE
        XL=X2
        XH=X1
        SWAP=FL
        FL=FH
        FH=SWAP
      ENDIF
      DX=XH-XL
      DO 11 J=1,MAXIT
        RTFLSP=XL+DX*FL/(FL-FH)
        F=FUNC(RTFLSP)
        IF(F.LT.0.) THEN
          DEL=XL-RTFLSP
          XL=RTFLSP
          FL=F
        ELSE
          DEL=XH-RTFLSP
          XH=RTFLSP
          FH=F
        ENDIF
        DX=XH-XL
        IF(ABS(DEL).LT.XACC.OR.F.EQ.0.)RETURN
11    CONTINUE
      PAUSE 'RTFLSP exceed maximum iterations'
      END
      FUNCTION RTNEWT(FUNCD,X1,X2,XACC)
      PARAMETER (JMAX=20)
      RTNEWT=.5*(X1+X2)
      DO 11 J=1,JMAX
        CALL FUNCD(RTNEWT,F,DF)
        DX=F/DF
        RTNEWT=RTNEWT-DX
        IF((X1-RTNEWT)*(RTNEWT-X2).LT.0.)PAUSE 'jumped out of brackets'
        IF(ABS(DX).LT.XACC) RETURN
11    CONTINUE
      PAUSE 'RTNEWT exceeding maximum iterations'
      END
      FUNCTION RTSAFE(FUNCD,X1,X2,XACC)
      PARAMETER (MAXIT=100)
      CALL FUNCD(X1,FL,DF)
      CALL FUNCD(X2,FH,DF)
      IF(FL*FH.GE.0.) PAUSE 'root must be bracketed'
      IF(FL.LT.0.)THEN
        XL=X1
        XH=X2
      ELSE
        XH=X1
        XL=X2
        SWAP=FL
        FL=FH
        FH=SWAP
      ENDIF
      RTSAFE=.5*(X1+X2)
      DXOLD=ABS(X2-X1)
      DX=DXOLD
      CALL FUNCD(RTSAFE,F,DF)
      DO 11 J=1,MAXIT
        IF(((RTSAFE-XH)*DF-F)*((RTSAFE-XL)*DF-F).GE.0.
     *      .OR. ABS(2.*F).GT.ABS(DXOLD*DF) ) THEN
          DXOLD=DX
          DX=0.5*(XH-XL)
          RTSAFE=XL+DX
          IF(XL.EQ.RTSAFE)RETURN
        ELSE
          DXOLD=DX
          DX=F/DF
          TEMP=RTSAFE
          RTSAFE=RTSAFE-DX
          IF(TEMP.EQ.RTSAFE)RETURN
        ENDIF
        IF(ABS(DX).LT.XACC) RETURN
        CALL FUNCD(RTSAFE,F,DF)
        IF(F.LT.0.) THEN
          XL=RTSAFE
          FL=F
        ELSE
          XH=RTSAFE
          FH=F
        ENDIF
11    CONTINUE
      PAUSE 'RTSAFE exceeding maximum iterations'
      RETURN
      END
      FUNCTION RTSEC(FUNC,X1,X2,XACC)
      PARAMETER (MAXIT=30)
      FL=FUNC(X1)
      F=FUNC(X2)
      IF(ABS(FL).LT.ABS(F))THEN
        RTSEC=X1
        XL=X2
        SWAP=FL
        FL=F
        F=SWAP
      ELSE
        XL=X1
        RTSEC=X2
      ENDIF
      DO 11 J=1,MAXIT
        DX=(XL-RTSEC)*F/(F-FL)
        XL=RTSEC
        FL=F
        RTSEC=RTSEC+DX
        F=FUNC(RTSEC)
        IF(ABS(DX).LT.XACC.OR.F.EQ.0.)RETURN
11    CONTINUE
      PAUSE 'RTSEC exceed maximum iterations'
      END
      SUBROUTINE RZEXTR(IEST,XEST,YEST,YZ,DY,NV,NUSE)
      PARAMETER (IMAX=11,NMAX=10,NCOL=7)
      DIMENSION X(IMAX),YEST(NV),YZ(NV),DY(NV),D(NMAX,NCOL),FX(NCOL)
      X(IEST)=XEST
      IF(IEST.EQ.1) THEN
        DO 11 J=1,NV
          YZ(J)=YEST(J)
          D(J,1)=YEST(J)
          DY(J)=YEST(J)
11      CONTINUE
      ELSE
        M1=MIN(IEST,NUSE)
        DO 12 K=1,M1-1
          FX(K+1)=X(IEST-K)/XEST
12      CONTINUE
        DO 14 J=1,NV
          YY=YEST(J)
          V=D(J,1)
          C=YY
          D(J,1)=YY
          DO 13 K=2,M1
            B1=FX(K)*V
            B=B1-C
            IF(B.NE.0.) THEN
              B=(C-V)/B
              DDY=C*B
              C=B1*B
            ELSE
              DDY=V
            ENDIF
            V=D(J,K)
            D(J,K)=DDY
            YY=YY+DDY
13        CONTINUE
          DY(J)=DDY
          YZ(J)=YY
14      CONTINUE
      ENDIF
      RETURN
      END
