      SUBROUTINE CALDAT(JULIAN,MM,ID,IYYY)
      PARAMETER (IGREG=2299161)
      IF(JULIAN.GE.IGREG)THEN
        JALPHA=INT(((JULIAN-1867216)-0.25)/36524.25)
        JA=JULIAN+1+JALPHA-INT(0.25*JALPHA)
      ELSE
        JA=JULIAN
      ENDIF
      JB=JA+1524
      JC=INT(6680.+((JB-2439870)-122.1)/365.25)
      JD=365*JC+INT(0.25*JC)
      JE=INT((JB-JD)/30.6001)
      ID=JB-JD-INT(30.6001*JE)
      MM=JE-1
      IF(MM.GT.12)MM=MM-12
      IYYY=JC-4715
      IF(MM.GT.2)IYYY=IYYY-1
      IF(IYYY.LE.0)IYYY=IYYY-1
      RETURN
      END
      FUNCTION CEL(QQC,PP,AA,BB)
      PARAMETER (CA=.0003, PIO2=1.5707963268)
      IF(QQC.EQ.0.)PAUSE 'failure in CEL'
      QC=ABS(QQC)
      A=AA
      B=BB
      P=PP
      E=QC
      EM=1.
      IF(P.GT.0.)THEN
        P=SQRT(P)
        B=B/P
      ELSE
        F=QC*QC
        Q=1.-F
        G=1.-P
        F=F-P
        Q=Q*(B-A*P)
        P=SQRT(F/G)
        A=(A-B)/G
        B=-Q/(G*G*P)+A*P
      ENDIF
1     F=A
      A=A+B/P
      G=E/P
      B=B+F*G
      B=B+B
      P=G+P
      G=EM
      EM=QC+EM
      IF(ABS(G-QC).GT.G*CA)THEN
        QC=SQRT(E)
        QC=QC+QC
        E=QC*EM
        GO TO 1
      ENDIF
      CEL=PIO2*(B+A*EM)/(EM*(EM+P))
      RETURN
      END
      SUBROUTINE CHDER(A,B,C,CDER,N)
      DIMENSION C(N),CDER(N)
      CDER(N)=0.
      CDER(N-1)=2*(N-1)*C(N)
      IF(N.GE.3)THEN
        DO 11 J=N-2,1,-1
          CDER(J)=CDER(J+2)+2*J*C(J+1)
11      CONTINUE
      ENDIF
      CON=2./(B-A)
      DO 12 J=1,N
        CDER(J)=CDER(J)*CON
12    CONTINUE
      RETURN
      END
      FUNCTION CHEBEV(A,B,C,M,X)
      DIMENSION C(M)
      IF ((X-A)*(X-B).GT.0.) PAUSE 'X not in range.'
      D=0.
      DD=0.
      Y=(2.*X-A-B)/(B-A)
      Y2=2.*Y
      DO 11 J=M,2,-1
        SV=D
        D=Y2*D-DD+C(J)
        DD=SV
11    CONTINUE
      CHEBEV=Y*D-DD+0.5*C(1)
      RETURN
      END
      SUBROUTINE CHEBFT(A,B,C,N,FUNC)
      REAL*8 PI
      PARAMETER (NMAX=50, PI=3.141592653589793D0)
      REAL*8 SUM
      DIMENSION C(N),F(NMAX)
      BMA=0.5*(B-A)
      BPA=0.5*(B+A)
      DO 11 K=1,N
        Y=COS(PI*(K-0.5)/N)
        F(K)=FUNC(Y*BMA+BPA)
11    CONTINUE
      FAC=2./N
      DO 13 J=1,N
        SUM=0.D0
        DO 12 K=1,N
          SUM=SUM+F(K)*COS((PI*(J-1))*((K-0.5D0)/N))
12      CONTINUE
        C(J)=FAC*SUM
13    CONTINUE
      RETURN
      END
      SUBROUTINE CHEBPC(C,D,N)
      PARAMETER (NMAX=50)
      DIMENSION C(N),D(N),DD(NMAX)
      DO 11 J=1,N
        D(J)=0.
        DD(J)=0.
11    CONTINUE
      D(1)=C(N)
      DO 13 J=N-1,2,-1
        DO 12 K=N-J+1,2,-1
          SV=D(K)
          D(K)=2.*D(K-1)-DD(K)
          DD(K)=SV
12      CONTINUE
        SV=D(1)
        D(1)=-DD(1)+C(J)
        DD(1)=SV
13    CONTINUE
      DO 14 J=N,2,-1
        D(J)=D(J-1)-DD(J)
14    CONTINUE
      D(1)=-DD(1)+0.5*C(1)
      RETURN
      END
      SUBROUTINE CHINT(A,B,C,CINT,N)
      DIMENSION C(N),CINT(N)
      CON=0.25*(B-A)
      SUM=0.
      FAC=1.
      DO 11 J=2,N-1
        CINT(J)=CON*(C(J-1)-C(J+1))/(J-1)
        SUM=SUM+FAC*CINT(J)
        FAC=-FAC
11    CONTINUE
      CINT(N)=CON*C(N-1)/(N-1)
      SUM=SUM+FAC*CINT(N)
      CINT(1)=2.*SUM
      RETURN
      END
      SUBROUTINE CHSONE(BINS,EBINS,NBINS,KNSTRN,DF,CHSQ,PROB)
      DIMENSION BINS(NBINS),EBINS(NBINS)
      DF=NBINS-1-KNSTRN
      CHSQ=0.
      DO 11 J=1,NBINS
        IF(EBINS(J).LE.0.)PAUSE 'bad expected number'
        CHSQ=CHSQ+(BINS(J)-EBINS(J))**2/EBINS(J)
11    CONTINUE
      PROB=GAMMQ(0.5*DF,0.5*CHSQ)
      RETURN
      END
      SUBROUTINE CHSTWO(BINS1,BINS2,NBINS,KNSTRN,DF,CHSQ,PROB)
      DIMENSION BINS1(NBINS),BINS2(NBINS)
      DF=NBINS-1-KNSTRN
      CHSQ=0.
      DO 11 J=1,NBINS
        IF(BINS1(J).EQ.0..AND.BINS2(J).EQ.0.)THEN
          DF=DF-1.
        ELSE
          CHSQ=CHSQ+(BINS1(J)-BINS2(J))**2/(BINS1(J)+BINS2(J))
        ENDIF
11    CONTINUE
      PROB=GAMMQ(0.5*DF,0.5*CHSQ)
      RETURN
      END
      SUBROUTINE CNTAB1(NN,NI,NJ,CHISQ,DF,PROB,CRAMRV,CCC)
      PARAMETER (MAXI=100,MAXJ=100,TINY=1.E-30)
      DIMENSION NN(NI,NJ),SUMI(MAXI),SUMJ(MAXJ)
      SUM=0
      NNI=NI
      NNJ=NJ
      DO 12 I=1,NI
        SUMI(I)=0.
        DO 11 J=1,NJ
          SUMI(I)=SUMI(I)+NN(I,J)
          SUM=SUM+NN(I,J)
11      CONTINUE
        IF(SUMI(I).EQ.0.)NNI=NNI-1
12    CONTINUE
      DO 14 J=1,NJ
        SUMJ(J)=0.
        DO 13 I=1,NI
          SUMJ(J)=SUMJ(J)+NN(I,J)
13      CONTINUE
        IF(SUMJ(J).EQ.0.)NNJ=NNJ-1
14    CONTINUE
      DF=NNI*NNJ-NNI-NNJ+1
      CHISQ=0.
      DO 16 I=1,NI
        DO 15 J=1,NJ
          EXPCTD=SUMJ(J)*SUMI(I)/SUM
          CHISQ=CHISQ+(NN(I,J)-EXPCTD)**2/(EXPCTD+TINY)
15      CONTINUE
16    CONTINUE
      PROB=GAMMQ(0.5*DF,0.5*CHISQ)
      CRAMRV=SQRT(CHISQ/(SUM*MIN(NNI-1,NNJ-1)))
      CCC=SQRT(CHISQ/(CHISQ+SUM))
      RETURN
      END
      SUBROUTINE CNTAB2(NN,NI,NJ,H,HX,HY,HYGX,HXGY,UYGX,UXGY,UXY)
      PARAMETER (MAXI=100,MAXJ=100,TINY=1.E-30)
      DIMENSION NN(NI,NJ),SUMI(MAXI),SUMJ(MAXJ)
      SUM=0
      DO 12 I=1,NI
        SUMI(I)=0.0
        DO 11 J=1,NJ
          SUMI(I)=SUMI(I)+NN(I,J)
          SUM=SUM+NN(I,J)
11      CONTINUE
12    CONTINUE
      DO 14 J=1,NJ
        SUMJ(J)=0.
        DO 13 I=1,NI
          SUMJ(J)=SUMJ(J)+NN(I,J)
13      CONTINUE
14    CONTINUE
      HX=0.
      DO 15 I=1,NI
        IF(SUMI(I).NE.0.)THEN
          P=SUMI(I)/SUM
          HX=HX-P*ALOG(P)
        ENDIF
15    CONTINUE
      HY=0.
      DO 16 J=1,NJ
        IF(SUMJ(J).NE.0.)THEN
          P=SUMJ(J)/SUM
          HY=HY-P*ALOG(P)
        ENDIF
16    CONTINUE
      H=0.
      DO 18 I=1,NI
        DO 17 J=1,NJ
          IF(NN(I,J).NE.0)THEN
            P=NN(I,J)/SUM
            H=H-P*ALOG(P)
          ENDIF
17      CONTINUE
18    CONTINUE
      HYGX=H-HX
      HXGY=H-HY
      UYGX=(HY-HYGX)/(HY+TINY)
      UXGY=(HX-HXGY)/(HX+TINY)
      UXY=2.*(HX+HY-H)/(HX+HY+TINY)
      RETURN
      END
      SUBROUTINE CONVLV(DATA,N,RESPNS,M,ISIGN,ANS)
      PARAMETER(NMAX=8192)
      DIMENSION DATA(N),RESPNS(N)
      COMPLEX FFT(NMAX),ANS(N)
      DO 11 I=1,(M-1)/2
        RESPNS(N+1-I)=RESPNS(M+1-I)
11    CONTINUE
      DO 12 I=(M+3)/2,N-(M-1)/2
        RESPNS(I)=0.0
12    CONTINUE
      CALL TWOFFT(DATA,RESPNS,FFT,ANS,N)
      NO2=N/2
      DO 13 I=1,NO2+1
        IF (ISIGN.EQ.1) THEN
          ANS(I)=FFT(I)*ANS(I)/NO2
        ELSE IF (ISIGN.EQ.-1) THEN
          IF (CABS(ANS(I)).EQ.0.0) PAUSE 'deconvolving at a response zer
     *o'
          ANS(I)=FFT(I)/ANS(I)/NO2
        ELSE
          PAUSE 'no meaning for ISIGN'
        ENDIF
13    CONTINUE
      ANS(1)=CMPLX(REAL(ANS(1)),REAL(ANS(NO2+1)))
      CALL REALFT(ANS,NO2,-1)
      RETURN
      END
      SUBROUTINE CORREL(DATA1,DATA2,N,ANS)
      PARAMETER(NMAX=8192)
      DIMENSION DATA1(N),DATA2(N)
      COMPLEX FFT(NMAX),ANS(N)
      CALL TWOFFT(DATA1,DATA2,FFT,ANS,N)
      NO2=FLOAT(N)/2.0
      DO 11 I=1,N/2+1
        ANS(I)=FFT(I)*CONJG(ANS(I))/NO2
11    CONTINUE
      ANS(1)=CMPLX(REAL(ANS(1)),REAL(ANS(N/2+1)))
      CALL REALFT(ANS,N/2,-1)
      RETURN
      END
      SUBROUTINE COSFT(Y,N,ISIGN)
      REAL*8 WR,WI,WPR,WPI,WTEMP,THETA
      DIMENSION Y(N)
      THETA=3.14159265358979D0/DBLE(N)
      WR=1.0D0
      WI=0.0D0
      WPR=-2.0D0*DSIN(0.5D0*THETA)**2
      WPI=DSIN(THETA)
      SUM=Y(1)
      M=N/2
      DO 11 J=1,M-1
        WTEMP=WR
        WR=WR*WPR-WI*WPI+WR
        WI=WI*WPR+WTEMP*WPI+WI
        Y1=0.5*(Y(J+1)+Y(N-J+1))
        Y2=(Y(J+1)-Y(N-J+1))
        Y(J+1)=Y1-WI*Y2
        Y(N-J+1)=Y1+WI*Y2
        SUM=SUM+WR*Y2
11    CONTINUE
      CALL REALFT(Y,M,+1)
      Y(2)=SUM
      DO 12 J=4,N,2
        SUM=SUM+Y(J)
        Y(J)=SUM
12    CONTINUE
      IF (ISIGN.EQ.-1) THEN
        EVEN=Y(1)
        ODD=Y(2)
        DO 13 I=3,N-1,2
          EVEN=EVEN+Y(I)
          ODD=ODD+Y(I+1)
13      CONTINUE
        ENF0=2.0*(EVEN-ODD)
        SUMO=Y(1)-ENF0
        SUME=(2.0*ODD/FLOAT(N))-SUMO
        Y(1)=0.5*ENF0
        Y(2)=Y(2)-SUME
        DO 14 I=3,N-1,2
          Y(I)=Y(I)-SUMO
          Y(I+1)=Y(I+1)-SUME
14      CONTINUE
      ENDIF
      RETURN
      END
      SUBROUTINE COVSRT(COVAR,NCVM,MA,LISTA,MFIT)
      DIMENSION COVAR(NCVM,NCVM),LISTA(MFIT)
      DO 12 J=1,MA-1
        DO 11 I=J+1,MA
          COVAR(I,J)=0.
11      CONTINUE
12    CONTINUE
      DO 14 I=1,MFIT-1
        DO 13 J=I+1,MFIT
          IF(LISTA(J).GT.LISTA(I)) THEN
            COVAR(LISTA(J),LISTA(I))=COVAR(I,J)
          ELSE
            COVAR(LISTA(I),LISTA(J))=COVAR(I,J)
          ENDIF
13      CONTINUE
14    CONTINUE
      SWAP=COVAR(1,1)
      DO 15 J=1,MA
        COVAR(1,J)=COVAR(J,J)
        COVAR(J,J)=0.
15    CONTINUE
      COVAR(LISTA(1),LISTA(1))=SWAP
      DO 16 J=2,MFIT
        COVAR(LISTA(J),LISTA(J))=COVAR(1,J)
16    CONTINUE
      DO 18 J=2,MA
        DO 17 I=1,J-1
          COVAR(I,J)=COVAR(J,I)
17      CONTINUE
18    CONTINUE
      RETURN
      END
      SUBROUTINE CRANK(N,W,S)
      DIMENSION W(N)
      S=0.
      J=1
1     IF(J.LT.N)THEN
        IF(W(J+1).NE.W(J))THEN
          W(J)=J
          J=J+1
        ELSE
          DO 11 JT=J+1,N
            IF(W(JT).NE.W(J))GO TO 2
11        CONTINUE
          JT=N+1
2         RANK=0.5*(J+JT-1)
          DO 12 JI=J,JT-1
            W(JI)=RANK
12        CONTINUE
          T=JT-J
          S=S+T**3-T
          J=JT
        ENDIF
      GO TO 1
      ENDIF
      IF(J.EQ.N)W(N)=N
      RETURN
      END
