      SUBROUTINE TOEPLZ(R,X,Y,N)
      PARAMETER (NMAX=100)
      DIMENSION R(*),X(N),Y(N),G(NMAX),H(NMAX)
      IF(R(N).EQ.0.) GO TO 99
      X(1)=Y(1)/R(N)
      IF(N.EQ.1)RETURN
      G(1)=R(N-1)/R(N)
      H(1)=R(N+1)/R(N)
      DO 15 M=1,N
        M1=M+1
        SXN=-Y(M1)
        SD=-R(N)
        DO 11 J=1,M
          SXN=SXN+R(N+M1-J)*X(J)
          SD=SD+R(N+M1-J)*G(M-J+1)
11      CONTINUE
        IF(SD.EQ.0.)GO TO 99
        X(M1)=SXN/SD
        DO 12 J=1,M
          X(J)=X(J)-X(M1)*G(M-J+1)
12      CONTINUE
        IF(M1.EQ.N)RETURN
        SGN=-R(N-M1)
        SHN=-R(N+M1)
        SGD=-R(N)
        DO 13 J=1,M
          SGN=SGN+R(N+J-M1)*G(J)
          SHN=SHN+R(N+M1-J)*H(J)
          SGD=SGD+R(N+J-M1)*H(M-J+1)
13      CONTINUE
        IF(SD.EQ.0..OR.SGD.EQ.0.)GO TO 99
        G(M1)=SGN/SGD
        H(M1)=SHN/SD
        K=M
        M2=(M+1)/2
        PP=G(M1)
        QQ=H(M1)
        DO 14 J=1,M2
          PT1=G(J)
          PT2=G(K)
          QT1=H(J)
          QT2=H(K)
          G(J)=PT1-PP*QT2
          G(K)=PT2-PP*QT1
          H(J)=QT1-QQ*PT2
          H(K)=QT2-QQ*PT1
          K=K-1
14      CONTINUE
15    CONTINUE
      PAUSE 'never get here'
99    PAUSE 'Levinson method fails: singular principal minor'
      END
      SUBROUTINE TPTEST(DATA1,DATA2,N,T,PROB)
      DIMENSION DATA1(N),DATA2(N)
      CALL AVEVAR(DATA1,N,AVE1,VAR1)
      CALL AVEVAR(DATA2,N,AVE2,VAR2)
      COV=0.
      DO 11 J=1,N
        COV=COV+(DATA1(J)-AVE1)*(DATA2(J)-AVE2)
11    CONTINUE
      DF=N-1
      COV=COV/DF
      SD=SQRT((VAR1+VAR2-2.*COV)/N)
      T=(AVE1-AVE2)/SD
      PROB=BETAI(0.5*DF,0.5,DF/(DF+T**2))
      RETURN
      END
      SUBROUTINE TQLI(D,E,N,NP,Z)
      DIMENSION D(NP),E(NP),Z(NP,NP)
      IF (N.GT.1) THEN
        DO 11 I=2,N
          E(I-1)=E(I)
11      CONTINUE
        E(N)=0.
        DO 15 L=1,N
          ITER=0
1         DO 12 M=L,N-1
            DD=ABS(D(M))+ABS(D(M+1))
            IF (ABS(E(M))+DD.EQ.DD) GO TO 2
12        CONTINUE
          M=N
2         IF(M.NE.L)THEN
            IF(ITER.EQ.30)PAUSE 'too many iterations'
            ITER=ITER+1
            G=(D(L+1)-D(L))/(2.*E(L))
            R=SQRT(G**2+1.)
            G=D(M)-D(L)+E(L)/(G+SIGN(R,G))
            S=1.
            C=1.
            P=0.
            DO 14 I=M-1,L,-1
              F=S*E(I)
              B=C*E(I)
              IF(ABS(F).GE.ABS(G))THEN
                C=G/F
                R=SQRT(C**2+1.)
                E(I+1)=F*R
                S=1./R
                C=C*S
              ELSE
                S=F/G
                R=SQRT(S**2+1.)
                E(I+1)=G*R
                C=1./R  
                S=S*C
              ENDIF
              G=D(I+1)-P
              R=(D(I)-G)*S+2.*C*B
              P=S*R
              D(I+1)=G+P
              G=C*R-B
              DO 13 K=1,N
                F=Z(K,I+1)
                Z(K,I+1)=S*Z(K,I)+C*F
                Z(K,I)=C*Z(K,I)-S*F
13            CONTINUE
14          CONTINUE
            D(L)=D(L)-P
            E(L)=G
            E(M)=0.
            GO TO 1
          ENDIF
15      CONTINUE
      ENDIF
      RETURN
      END
      SUBROUTINE TRAPZD(FUNC,A,B,S,N)
      IF (N.EQ.1) THEN
        S=0.5*(B-A)*(FUNC(A)+FUNC(B))
        IT=1
      ELSE
        TNM=IT
        DEL=(B-A)/TNM
        X=A+0.5*DEL
        SUM=0.
        DO 11 J=1,IT
          SUM=SUM+FUNC(X)
          X=X+DEL
11      CONTINUE
        S=0.5*(S+(B-A)*SUM/TNM)
        IT=2*IT
      ENDIF
      RETURN
      END
      SUBROUTINE TRED2(A,N,NP,D,E)
      DIMENSION A(NP,NP),D(NP),E(NP)
      IF(N.GT.1)THEN
        DO 18 I=N,2,-1  
          L=I-1
          H=0.
          SCALE=0.
          IF(L.GT.1)THEN
            DO 11 K=1,L
              SCALE=SCALE+ABS(A(I,K))
11          CONTINUE
            IF(SCALE.EQ.0.)THEN
              E(I)=A(I,L)
            ELSE
              DO 12 K=1,L
                A(I,K)=A(I,K)/SCALE
                H=H+A(I,K)**2
12            CONTINUE
              F=A(I,L)
              G=-SIGN(SQRT(H),F)
              E(I)=SCALE*G
              H=H-F*G
              A(I,L)=F-G
              F=0.
              DO 15 J=1,L
                A(J,I)=A(I,J)/H
                G=0.
                DO 13 K=1,J
                  G=G+A(J,K)*A(I,K)
13              CONTINUE
                IF(L.GT.J)THEN
                  DO 14 K=J+1,L
                    G=G+A(K,J)*A(I,K)
14                CONTINUE
                ENDIF
                E(J)=G/H
                F=F+E(J)*A(I,J)
15            CONTINUE
              HH=F/(H+H)
              DO 17 J=1,L
                F=A(I,J)
                G=E(J)-HH*F
                E(J)=G
                DO 16 K=1,J
                  A(J,K)=A(J,K)-F*E(K)-G*A(I,K)
16              CONTINUE
17            CONTINUE
            ENDIF
          ELSE
            E(I)=A(I,L)
          ENDIF
          D(I)=H
18      CONTINUE
      ENDIF
      D(1)=0.
      E(1)=0.
      DO 23 I=1,N
        L=I-1
        IF(D(I).NE.0.)THEN
          DO 21 J=1,L
            G=0.
            DO 19 K=1,L
              G=G+A(I,K)*A(K,J)
19          CONTINUE
            DO 20 K=1,L
              A(K,J)=A(K,J)-G*A(K,I)
20          CONTINUE
21        CONTINUE
        ENDIF
        D(I)=A(I,I)
        A(I,I)=1.
        IF(L.GE.1)THEN
          DO 22 J=1,L
            A(I,J)=0.
            A(J,I)=0.
22        CONTINUE
        ENDIF
23    CONTINUE
      RETURN
      END
      SUBROUTINE TRIDAG(A,B,C,R,U,N)
      PARAMETER (NMAX=100)
      DIMENSION GAM(NMAX),A(N),B(N),C(N),R(N),U(N)
      IF(B(1).EQ.0.)PAUSE
      BET=B(1)
      U(1)=R(1)/BET
      DO 11 J=2,N
        GAM(J)=C(J-1)/BET
        BET=B(J)-A(J)*GAM(J)
        IF(BET.EQ.0.)PAUSE
        U(J)=(R(J)-A(J)*U(J-1))/BET
11    CONTINUE
      DO 12 J=N-1,1,-1
        U(J)=U(J)-GAM(J+1)*U(J+1)
12    CONTINUE
      RETURN
      END
      SUBROUTINE TTEST(DATA1,N1,DATA2,N2,T,PROB)
      DIMENSION DATA1(N1),DATA2(N2)
      CALL AVEVAR(DATA1,N1,AVE1,VAR1)
      CALL AVEVAR(DATA2,N2,AVE2,VAR2)
      DF=N1+N2-2
      VAR=((N1-1)*VAR1+(N2-1)*VAR2)/DF
      T=(AVE1-AVE2)/SQRT(VAR*(1./N1+1./N2))
      PROB=BETAI(0.5*DF,0.5,DF/(DF+T**2))
      RETURN
      END
      SUBROUTINE TUTEST(DATA1,N1,DATA2,N2,T,PROB)
      DIMENSION DATA1(N1),DATA2(N2)
      CALL AVEVAR(DATA1,N1,AVE1,VAR1)
      CALL AVEVAR(DATA2,N2,AVE2,VAR2)
      T=(AVE1-AVE2)/SQRT(VAR1/N1+VAR2/N2)
      DF=(VAR1/N1+VAR2/N2)**2/((VAR1/N1)**2/(N1-1)+(VAR2/N2)**2/(N2-1))
      PROB=BETAI(0.5*DF,0.5,DF/(DF+T**2))
      RETURN
      END
      SUBROUTINE TWOFFT(DATA1,DATA2,FFT1,FFT2,N)
      DIMENSION DATA1(N),DATA2(N)
      COMPLEX FFT1(N),FFT2(N),H1,H2,C1,C2
      C1=CMPLX(0.5,0.0)
      C2=CMPLX(0.0,-0.5)
      DO 11 J=1,N
        FFT1(J)=CMPLX(DATA1(J),DATA2(J))
11    CONTINUE
      CALL FOUR1(FFT1,N,1)
      FFT2(1)=CMPLX(AIMAG(FFT1(1)),0.0)
      FFT1(1)=CMPLX(REAL(FFT1(1)),0.0)
      N2=N+2
      DO 12 J=2,N/2+1
        H1=C1*(FFT1(J)+CONJG(FFT1(N2-J)))
        H2=C2*(FFT1(J)-CONJG(FFT1(N2-J)))
        FFT1(J)=H1
        FFT1(N2-J)=CONJG(H1)
        FFT2(J)=H2
        FFT2(N2-J)=CONJG(H2)
12    CONTINUE
      RETURN
      END
