C
      PROGRAM correct
C
      PARAMETER(MWV=33000)
      CHARACTER INFILE*20
      REAL time(MWV)
      REAL XR1(MWV),XR2(MWV),XR3(MWV),X1(MWV),X2(MWV),X3(MWV)
      COMPLEX*8 FX1(MWV),FX2(MWV),FX3(MWV)
C
      WRITE(*,*)
      WRITE(*,*) 'FFT and Band-Pass Filter for Taiwan Data'
      WRITE(*,*) 'Read from "***_dis" '
      WRITE(*,*) '------------------------------------------------------
     +-----------'
      WRITE(*,*)
C
C  Prompt for input RECON'ed file
C
      WRITE(*,1010) 
 1010 FORMAT(1X,'Enter name "***_dis')
      READ(*,1020) INFILE
 1020 FORMAT(A20)
C      
C *** Read input data file  ***
C
      OPEN(11,FILE=INFILE)
        i=1
   10   read(11,*,end=11)time(i),xr1(i),xr2(i),xr3(i)
        i=i+1
        goto 10
   11 CLOSE(11)
c
      NT=i-1
      dt=time(2)-time(1)
      write(*,*)'NT and DT=',NT,DT
      dur=(nt-1)*dt
      write(*,*)'Duration (sec)=',dur
c
c *** baseline correction *** 
c
      write(*,*)
      write(*,*)'Enter Initial and Final Time for the first half'
      read(*,*)tinit,tfinal
      rat1=0.0
      rat2=0.0
      rat3=0.0
      nrat=0
      do 100 i=1,nt
        time1=dt*(i-1)
        if((time1.le.tinit).or.(time1.gt.tfinal))goto 100
        amp11=xr1(i-1)
        amp21=xr2(i-1)
        amp31=xr3(i-1)
        amp12=xr1(i)
        amp22=xr2(i)
        amp32=xr3(i)
        rat1=rat1+(amp12-amp11)/dt
        rat2=rat2+(amp22-amp21)/dt
        rat3=rat3+(amp32-amp31)/dt
        nrat=nrat+1
  100 continue
      rat11=rat1/nrat
      rat12=rat2/nrat
      rat13=rat3/nrat
      write(*,*)'First rat1,rat2,rat3: ',rat11,rat12,rat13
      do 110 i=1,nt
        time1=dt*(i-1)
        if(time1.lt.tinit)goto 110
        init1=i
        x11=xr1(i)
        x12=xr2(i)
        x13=xr3(i)
        goto 120
  110 continue
  120 write(*,*)'init,x11,x12,x13 ',init1,x11,x12,x13
      b11=x11-rat11*time(init1)
      b12=x12-rat12*time(init1)
      b13=x13-rat13*time(init1)
      write(*,*)'b11,b12,b13 ',b11,b12,b13
c
      write(*,*)
      write(*,*)'Enter Initial and Final Time for the second half'
      read(*,*)tinit,tfinal
      rat1=0.0
      rat2=0.0
      rat3=0.0
      nrat=0
      do 140 i=1,nt
        time1=dt*(i-1)
        if((time1.le.tinit).or.(time1.gt.tfinal))goto 140
        amp11=xr1(i-1)
        amp21=xr2(i-1)
        amp31=xr3(i-1)
        amp12=xr1(i)
        amp22=xr2(i)
        amp32=xr3(i)
        rat1=rat1+(amp12-amp11)/dt
        rat2=rat2+(amp22-amp21)/dt
        rat3=rat3+(amp32-amp31)/dt
        nrat=nrat+1
  140 continue
      rat21=rat1/nrat
      rat22=rat2/nrat
      rat23=rat3/nrat
      write(*,*)'Second rat1,rat2,rat3: ',rat21,rat22,rat23
      do 150 i=1,nt
        time1=dt*(i-1)
        if(time1.lt.tinit)goto 150
        init2=i
        x21=xr1(i)
        x22=xr2(i)
        x23=xr3(i)
        goto 155
  150 continue
  155 write(*,*)'init,x21,x22,x23 ',init2,x21,x22,x23
      b21=x21-rat21*time(init2)
      b22=x22-rat22*time(init2)
      b23=x23-rat23*time(init2)
      write(*,*)'b21,b22,b23 ',b21,b22,b23
      cross1=(b21-b11)/(RAT11-RAT21)
      cross2=(b22-b12)/(RAT12-RAT22)
      cross3=(b23-b13)/(RAT13-RAT23)
      WRITE(*,*)'crossing time (s)= ',cross1,cross2,cross3 
      write(*,*)'Is those ok? (yes=0,no=1)'
      read(*,*)ncross
      if(ncross.eq.1)then
         write(*,*)'Enter cross1,cross2,cross3'
         read(*,*)cross1,cross2,cross3
         do 160 i=1,nt
           time1=dt*(i-1)
           if(time1.lt.cross1)goto 160
           init=i
           x21=xr1(i)
           goto 161
  160    continue
  161    b21=x21-rat21*time(init)
         do 162 i=1,nt
           time1=dt*(i-1)
           if(time1.lt.cross2)goto 162
           init=i
           x22=xr2(i)
           goto 163
  162    continue
  163    b22=x22-rat22*time(init)
         do 164 i=1,nt
           time1=dt*(i-1)
           if(time1.lt.cross3)goto 164
           init=i
           x23=xr3(i)
           goto 165
  164    continue
  165    b23=x23-rat23*time(init)
         write(*,*)'b21,b22,b23 ',b21,b22,b23
      end if
c
      do 170 i=1,nt
        time1=dt*(i-1)
        if(time1.lt.cross1)then
          xr1(i)=xr1(i)-(rat11*time1+b11)
          x211=xr1(i)
        else
          xr1(i)=xr1(i)-(rat21*time1+b21)
          if(ncross.eq.1)xr1(i)=xr1(i)+x211
        end if
        if(time1.lt.cross2)then
          xr2(i)=xr2(i)-(rat12*time1+b12)
          x221=xr2(i)
        else
          xr2(i)=xr2(i)-(rat22*time1+b22)
          if(ncross.eq.1)xr2(i)=xr2(i)+x221+x22
        end if
        if(time1.lt.cross3)then
          xr3(i)=xr3(i)-(rat13*time1+b13)
          x231=xr3(i)
        else
          xr3(i)=xr3(i)-(rat23*time1+b23)
          if(ncross.eq.1)xr2(i)=xr2(i)+x231
        end if
  170 continue
c
C *** FOURIER TRANSFORM ***
C
      NF=2
      DO 200 I=2,16
        NF=2*NF
        IF(NF.GE.NT)GOTO 30
  200 CONTINUE
      WRITE(*,*)' NT>2**16 !!!!!'
      STOP
   30 IF1=I
      WRITE(*,*)
      WRITE(*,*)' IF,NF=',IF1,NF
C
      DO 300 I=1,NT
        FX1(I)=CMPLX(XR1(I),0.0)
        FX2(I)=CMPLX(XR2(I),0.0)
        FX3(I)=CMPLX(XR3(I),0.0)
  300 CONTINUE
      DO 310 I=NT+1,NF
        FX1(I)=(0.0,0.0)
        FX2(I)=(0.0,0.0)
        FX3(I)=(0.0,0.0)
  310 CONTINUE
c
      CALL WVFAST(NF,FX1,MWV,-1)
      CALL WVFAST(NF,FX2,MWV,-1)
      CALL WVFAST(NF,FX3,MWV,-1)
      DUR=FLOAT(NF)*DT
      DO 400 I=1,NF
        FX1(I)=FX1(I)*DT
        FX2(I)=FX2(I)*DT
        FX3(I)=FX3(I)*DT
  400 CONTINUE
c
c      FOM=0.0
      DOM=6.2831853070e0/DUR
      NOM=NF/2+1
      EOM=DOM*NOM
      WRITE(*,*)
      WRITE(*,*)' DUR,DOM,NOM,EOM=',DUR,DOM,NOM,EOM
C
C *** BAND PATH ***
C
      WRITE(*,*)
      WRITE(*,*)' BAND PASS? (YES=1)?? :'
      READ(*,*)NBD
C
      IF(NBD.EQ.1)THEN
C
        WRITE(*,*)
        WRITE(*,*)'  TRAPEZOIDAL BAND PATH FILTER'
c        WRITE(*,*)'  INITIAL AND FINAL OM=',DOM*IOM,DOM*(NDATA-1)
        WRITE(*,*)'  ENTER OML1<OML2<OMH1<OMH2'
        READ(*,*)OML1,OML2,OMH1,OMH2
C
        DO 500 I = 1,NOM
           OMEGA=FLOAT(I-1)*DOM
           IF(OMEGA.LE.OML1)then
              FX1(I)=(0.0,0.0)
              FX2(I)=(0.0,0.0)
              FX3(I)=(0.0,0.0)
           end if
           IF(OMEGA.GE.OMH2)then
              FX1(I)=(0.0,0.0)
              FX2(I)=(0.0,0.0)
              FX3(I)=(0.0,0.0)
           end if
           AA=1.0
           IF(OMEGA.LT.OML2)AA=(OMEGA-OML1)/(OML2-OML1)
           IF(OMEGA.GT.OMH1)AA=(OMH2-OMEGA)/(OMH2-OMH1)
	   FX1(I) = AA*FX1(I)
	   FX2(I) = AA*FX2(I)
	   FX3(I) = AA*FX3(I)
  500   CONTINUE
C
      END IF
C
      WRITE(*,*)
      WRITE(*,*)' Integ.(=1), Double Integ.(=2) or Diff.(=3) ?:'
      READ(*,*)NIT
c
      IF(NIT.EQ.1)THEN
        DO 550 I = 2,NOM
           OMEGA=FLOAT(I-1)*DOM
	   FX1(I) = FX1(I)/cmplx(0.0,omega)
	   FX2(I) = FX2(I)/cmplx(0.0,omega)
	   FX3(I) = FX3(I)/cmplx(0.0,omega)
  550   CONTINUE
      END IF
c
      IF(NIT.EQ.2)THEN
        DO 560 I = 2,NOM
           OMEGA=FLOAT(I-1)*DOM
	   FX1(I) = FX1(I)/(-omega*omega)
	   FX2(I) = FX2(I)/(-omega*omega)
	   FX3(I) = FX3(I)/(-omega*omega)
  560   CONTINUE
      END IF
c
      IF(NIT.EQ.3)THEN
        DO 570 I = 2,NOM
           OMEGA=FLOAT(I-1)*DOM
	   FX1(I) = FX1(I)*cmplx(0.0,omega)
	   FX2(I) = FX2(I)*cmplx(0.0,omega)
 	   FX3(I) = FX3(I)*cmplx(0.0,omega)
 570   CONTINUE
      END IF
c
c * Output Fourier Spectrum data *
c
      write(*,*)
      write(*,*)'The Fourier Spectrum data are stored in "amp_dat"'
      open(10,file='amp_dat',status='unknown')
      DO 680 I = 2,NOM
        om=FLOAT(I-1)*DOM
        fr=om/6.283185307
        write(10,1900)om,fr,cabs(FX1(I)),cabs(FX2(I)),cabs(FX3(I))
  680 CONTINUE
 1900 format(f8.2,f8.3,3e12.4)
      close(10)
C
C *** INVERSE FOURIER TRANSFORM  ***
C
      DO 600 I=2,NF/2
           FX1(NF+2-I)=CMPLX(REAL(FX1(I)),-AIMAG(FX1(I)))
           FX2(NF+2-I)=CMPLX(REAL(FX2(I)),-AIMAG(FX2(I)))
           FX3(NF+2-I)=CMPLX(REAL(FX3(I)),-AIMAG(FX3(I)))
  600 CONTINUE
C
      CALL WVFAST(NF,FX1,MWV,1)
      CALL WVFAST(NF,FX2,MWV,1)
      CALL WVFAST(NF,FX3,MWV,1)
C
C *** WAVE DATA ***
C
      DO 700 I = 1,NF
         time(I) = FLOAT(I-1)*DT
         XR1(I)=REAL(FX1(I))/DUR
         XR2(I)=REAL(FX2(I))/DUR
         XR3(I)=REAL(FX3(I))/DUR
  700 CONTINUE
c
c *** Cut ***
c
      do i=1,nt
        x1(i)=xr1(i)
        x2(i)=xr2(i)
        x3(i)=xr3(i)
      end do

      write(*,*)
      write(*,*)'Cut data? (No=0,Yes=1)'
      read(*,*)icut
      if(icut.eq.1)then
        write(*,*)'Enter Starting Time (sec)'
        read(*,*)stime
        nstime=stime/dt
        write(*,*)'Number of Cutting Data =',nstime
        write(*,*)'Enter Duration (sec)'
        read(*,*)dur
        ndur=dur/dt+1
        write(*,*)'Number of data =',ndur
        do i=1,ndur
          if((i+nstime).le.0)then
            x1(i)=0.0
            x2(i)=0.0
            x3(i)=0.0
	  else
            x1(i)=xr1(i+nstime)
            x2(i)=xr2(i+nstime)
            x3(i)=xr3(i+nstime)
          end if
        end do
	nt=ndur
        write(*,*)'New nt =',nt
      end if
C
c *** output ***
c
      write(*,*)
      write(*,*)'Skip Data (no=0; yes=skip number)?'
      read(*,*)iskip
      iskip=iskip+1
      WRITE(*,*)'wave2_dat:' 
      OPEN(11,FILE='wave2_dat',STATUS='unknown')
      do i=1,nt,iskip
      write(11,2000)dt*(i-1),X1(I)-x1(1),x2(i)-x2(1),x3(i)-x3(1)
      end do
 2000 format(4e13.5)
      CLOSE(11)
c
        STOP
	END
C
C *** FOURIER TRANSFORM OR INVERSE FOURIER TRANSFORM
C
      SUBROUTINE WVFAST(N,A,MOM,IND)
C
C *** N   : NUMBER OF FOURIER COEFF. OR WAVE ( = 2**I )
C *** A(I): FOURIER COEFFICIENTS ( N/2+1 ) OR WAVE ( REAL PART )
C    ( FOURIER SPECTRUM = FOURIER COEFFICIENTS * DURATION/2.0 )
C *** MOM : N IN MAIN DIMMENSION
C         : IF MWV > 2**13(=8192), CAN NOT USE !!!
C *** IND : = -1 : TRANS.,   = 1 : INVERSE TRANS.
C
      COMPLEX*8 A(MOM),TEMP,THETA
C
C *** CALUCURATION OF FFT
C
      PAI=3.1415926540e0
      J=1
      DO 140 I=1,N
        IF(I.GE.J) GO TO 110
        TEMP=A(J)
        A(J)=A(I)
        A(I)=TEMP
  110   M=N/2
  120   IF(J.LE.M) GO TO 130
        J=J-M
        M=M/2
        IF(M.GE.2) GO TO 120
  130   J=J+M
  140 CONTINUE
      KMAX=1
  150 IF(KMAX.GE.N) RETURN
      ISTEP=KMAX*2
      DO 170 K=1,KMAX
        THETA=CMPLX(0.0,PAI*FLOAT(IND*(K-1))/FLOAT(KMAX))
        DO 160 I=K,N,ISTEP
          J=I+KMAX
          TEMP=A(J)*CEXP(THETA)
          A(J)=A(I)-TEMP
          A(I)=A(I)+TEMP
  160   CONTINUE
  170 CONTINUE
      KMAX=ISTEP
      GO TO 150
C
      END
