C
c *** Calculation of Transfer Function ***
c
      PROGRAM trans
C
      PARAMETER(MWV=33000)
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER INFILE*20
      dimension time(MWV)
      dimension XR1(MWV),XR2(MWV),XR3(MWV),X1(MWV),X2(MWV)
      COMPLEX*16 FX1(MWV),FX2(MWV),FX3(MWV)
      COMPLEX*16 S11(MWV),S12(MWV),S21(MWV),S22(MWV)
C
C  Prompt for input RECON'ed file
C
      WRITE(*,1010) 
 1010 FORMAT(1X,'Enter input file name (*.csv file)')
      READ(*,1020) INFILE
 1020 FORMAT(A20)
C
C *** Read input data file  ***
C
      OPEN(11,FILE=INFILE)
c      OPEN(11,FILE='047-test.csv')
        read(11,*)
        i=1
   10   read(11,*,end=11)time(i),x1(i),x2(i)
        i=i+1
        goto 10
   11 CLOSE(11)
c
      NT=i-1
      dt=time(2)-time(1)
      write(*,*)'Number of Time =',NT
      write(*,*)'Delta Time (s) =',DT
c
c      twin=10.0d0
c      tint=1.0d0
      write(*,*)
      write(*,*)'Enter Duration of Time Window (s) ='
      read(*,*)twin
      write(*,*)
      write(*,*)'Enter Interval Time among Time Windows (s) ='
      read(*,*)tint
      ntwin=int(twin/DT)
      ntint=int(tint/DT)
      write(*,*)
      write(*,*)'Number of Time for Time Window   =',ntwin
      write(*,*)'Number of Time for Time Interval =',ntint
      write(*,*)
      write(*,*)'Enter Starting Time of Time Windows (s) ='
      read(*,*)stime
c
      NF=2
      DO I=2,15
        NF=2*NF
        IF(NF.GE.ntwin)GOTO 20
      end do
      WRITE(*,*)' NT>2**15 !!!!!'
      STOP
   20 WRITE(*,*)'Number of time for FFT=',NF
      DUR=dFLOAT(NF)*DT
c      FOM=0.d0
      DOM=6.2831853070d0/DUR
      NOM=NF/2+1
      EOM=DOM*NOM
      WRITE(*,*)
      WRITE(*,*)' DUR,DOM,NOM,EOM=',DUR,DOM,NOM,EOM
      WRITE(*,*)' Highest Frequency (hz)=',EOM/6.2831853070d0
C
C *** BAND PATH ***
C
c      WRITE(*,*)
c      WRITE(*,*)' Band Pass Filter? (YES=1)?? :'
c      READ(*,*)NBD
C
c      IF(NBD.EQ.1)THEN
c        WRITE(*,*)
c        WRITE(*,*)'  TRAPEZOIDAL BAND PATH FILTER'
c        WRITE(*,*)'  ENTER fL1<fL2<fH1<fH2 (Hz)'
c        READ(*,*)fL1,fL2,fH1,fH2
c        OML1=fL1*6.2831853070d0
c        OML2=fL2*6.2831853070d0
c        OMH1=fH1*6.2831853070d0
c        OMH2=fH2*6.2831853070d0
c      END IF
c
      DO I=1,NOM
        S11(I)=(0.d0,0.d0)
        S12(I)=(0.d0,0.d0)
        S21(I)=(0.d0,0.d0)
        S22(I)=(0.d0,0.d0)
      end do
c
c *** Loop of Time Window ***
c
      write(*,*)
      write(*,*)'Time Window Wave Data are stored in "twin.csv"'
      open(11,file='twin.csv',status='unknown')
c
      nwind=0
      duration=dt*dfloat(ntwin)
      timeL1=0.d0
      timeL2=duration/5.d0
      timeH1=duration*4.d0/5.d0
      timeH2=duration
      write(*,*)
      write(*,*)'timeL1,L2,H1,H4 =',timeL1,timeL2,timeH1,timeH2
   30 write(*,*)'Calculating Time Window =',nwind+1
      swtime=dt*dfloat(ntint*nwind)
      if(swtime.lt.stime)goto35
c
      write(*,*)'Starting Time (s) =',swtime
      write(11,*)nwind+1,',input,output'
      do 100 i=1,ntwin
        if(i+ntint*nwind.gt.NT)goto 100
        t=dt*dfloat(i-1)
        AA=1.d0
        IF(t.LT.timeL2)AA=(t-timeL1)/timeL2
        IF(t.GT.timeH1)AA=(timeH2-t)/timeL2
        xr1(i)=AA*x1(i+ntint*nwind)
        xr2(i)=AA*x2(i+ntint*nwind)
        write(11,*)dt*dfloat(i+ntint*nwind-1),',',xr1(i),',',xr2(i)
  100 continue
C
C *** FOURIER TRANSFORM ***
C
      DO 300 I=1,ntwin
        FX1(I)=dCMPLX(XR1(I),0.d0)
        FX2(I)=dCMPLX(XR2(I),0.d0)
  300 CONTINUE
      DO 310 I=ntwin+1,NF
        FX1(I)=(0.d0,0.d0)
        FX2(I)=(0.d0,0.d0)
  310 CONTINUE
C
      CALL WVFAST(NF,FX1,MWV,-1)
      CALL WVFAST(NF,FX2,MWV,-1)
      DO 320 I=1,NF
        FX1(I)=FX1(I)*DT
        FX2(I)=FX2(I)*DT
  320 CONTINUE
C
c      IF(NBD.EQ.1)THEN
C
c        DO I = 1,NOM
c           OMEGA=dFLOAT(I-1)*DOM
c           IF(OMEGA.LE.OML1)then
c             FX1(I)=(0.d0,0.d0)
c             FX2(I)=(0.d0,0.d0)
c           end if
c           IF(OMEGA.GE.OMH2)then
c             FX1(I)=(0.d0,0.d0)
c             FX2(I)=(0.d0,0.d0)
c           end if
c           AA=1.d0
c           IF(OMEGA.LT.OML2)AA=(OMEGA-OML1)/(OML2-OML1)
c           IF(OMEGA.GT.OMH1)AA=(OMH2-OMEGA)/(OMH2-OMH1)
c	   FX1(I) = AA*FX1(I)
c	   FX2(I) = AA*FX2(I)
c        End Do
C
c      END IF
c
      DO 400 I=2,NOM
        S11(I)=S11(I)+FX1(I)*dconjg(FX1(I))
        S12(I)=S12(I)+FX1(I)*dconjg(FX2(I))
        S21(I)=S21(I)+FX2(I)*dconjg(FX1(I))
        S22(I)=S22(I)+FX2(I)*dconjg(FX2(I))
  400 CONTINUE
c
   35 nwind=nwind+1
      if(ntint*(nwind+1).lt.NT)goto 30
      close(11)
c
c *** Transfer Function ***
c
      DO 500 I=2,NOM
        if(cdabs(S11(I)).lt.1.d-10)then
          FX1(I)=(0.d0,0.d0)
          FX3(I)=(0.d0,0.d0)
        else
          FX1(I)=S21(I)/S11(I)
          FX3(I)=cdsqrt(S22(I)/S11(I))
        end if
        if(cdabs(S12(I)).lt.1.d-10)then
          FX2(I)=(0.d0,0.d0)
        else
          FX2(I)=S22(I)/S12(I)
        end if
  500 CONTINUE
C
c * Output Fourier Spectrum data *
c
      write(*,*)
      write(*,*)'Fourier Spectrum data are stored in "trans.csv"'
      open(10,file='trans.csv',status='unknown')
      write(10,*)'f (Hz),H1,H2,H3'
      DO 650 I = 2,NOM
        f=dFLOAT(I-1)*DOM/6.2831853070d0
        write(10,1900)f,cdabs(FX1(I)),cdabs(FX2(I)),cdabs(FX3(I))
  650 CONTINUE
 1900 format(f8.4,6(',',e12.4))
      close(10)
C
C *** INVERSE FOURIER TRANSFORM  ***
C
      DO 700 I=2,NF/2
           FX1(NF+2-I)=dconjg(FX1(I))
           FX2(NF+2-I)=dconjg(FX2(I))
           FX3(NF+2-I)=dconjg(FX3(I))
  700 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 800 I = 1,NF
         XR1(I)=dREAL(FX1(I))/DUR
         XR2(I)=dREAL(FX2(I))/DUR
         XR3(I)=dREAL(FX3(I))/DUR
  800 CONTINUE
C
c *** output ***
c
      write(*,*)
      WRITE(*,2010)
 2010 FORMAT(/1X,'Enter output file name of waveform (csv format file)')
      READ(*,2020) INFILE
 2020 FORMAT(A20)
      OPEN(11,FILE=INFILE,STATUS='unknown')
      write(11,*)'time (s),H1,H2,H3'
      do i=1,NF
      write(11,2000)dt*(i-1),Xr1(I),xr2(i),Xr3(I)
      end do
 2000 format(f6.3,4(',',f12.4))
      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
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 A(MOM),TEMP,THETA
C
C *** CALUCURATION OF FFT
C
      PAI=3.1415926540d0
      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=dCMPLX(0.d0,PAI*dFLOAT(IND*(K-1))/dFLOAT(KMAX))
        DO 160 I=K,N,ISTEP
          J=I+KMAX
          TEMP=A(J)*CdEXP(THETA)
          A(J)=A(I)-TEMP
          A(I)=A(I)+TEMP
  160   CONTINUE
  170 CONTINUE
      KMAX=ISTEP
      GO TO 150
C
      END
