C
      PROGRAM FRatio
C
      PARAMETER(MWV=33000)
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER INFILE1*90,head1*90,head2*90,head3*90,head4*90
      CHARACTER INFILE2*90,head5*90,head6*90,head7*90,head8*90
      dimension time1(MWV),time2(MWV)
      dimension XR1(MWV),XR2(MWV),XR3(MWV),X1(MWV),X2(MWV),X3(MWV)
      dimension YR1(MWV),YR2(MWV),YR3(MWV),Y1(MWV),Y2(MWV),Y3(MWV)
      COMPLEX*16 FX1(MWV),FX2(MWV),FX3(MWV)
      COMPLEX*16 FY1(MWV),FY2(MWV),FY3(MWV)
      COMPLEX*16 YX1(MWV),YX2(MWV),YX3(MWV)
C
C  Prompt for input RECON'ed file
C
      write(*,*)
      write(*,*)' Fourier Amplitude Ratio (file No.2/file No.1)'
      WRITE(*,1010) 
      READ(*,1020) INFILE1
      WRITE(*,1015)
      READ(*,1020) INFILE2
 1010 FORMAT(1X,/'Enter input file No.1 of waveform data')
 1015 FORMAT(1X,/'Enter input file No.2 of waveform data')
 1020 FORMAT(A90)
C
      OPEN(11,FILE=INFILE1)
        read(11,*)head1,head2,head3,head4
        i=1
   10   read(11,*,end=11)time1(i),xr1(i),xr2(i),xr3(i)
        i=i+1
        goto 10
   11 CLOSE(11)
c
      OPEN(12,FILE=INFILE2)
        read(12,*)head5,head6,head7,head8
        j=1
   12   read(12,*,end=13)time2(j),yr1(j),yr2(j),yr3(j)
        j=j+1
        goto 12
   13 CLOSE(12)
c
      NT1=i-1
      dt1=time1(2)-time1(1)
      NT2=j-1
      dt2=time2(2)-time2(1)
      write(*,*)
      write(*,*)'NT1 and DT1=',NT1,DT1
      write(*,*)'NT2 and DT2=',NT2,DT2
      write(*,*)' if NT1.ne.NT2 or DT1.ne.DT2, then stop!!'
      if(NT1.ne.NT2)stop
      if(DT1.ne.DT2)stop
      NT=NT1
      dt=dt1
C
C *** FOURIER TRANSFORM ***
C
      NF=2
      DO 200 I=2,15
        NF=2*NF
        IF(NF.GE.NT)GOTO 30
  200 CONTINUE
      WRITE(*,*)' NT>2**15 !!!!!'
      STOP
   30 IF1=I
      WRITE(*,*)
      WRITE(*,*)' IF1,NF=',IF1,NF
C
      DO 300 I=1,NT
        FX1(I)=dCMPLX(XR1(I),0.d0)
        FX2(I)=dCMPLX(XR2(I),0.d0)
        FX3(I)=dCMPLX(XR3(I),0.d0)
        FY1(I)=dCMPLX(YR1(I),0.d0)
        FY2(I)=dCMPLX(YR2(I),0.d0)
        FY3(I)=dCMPLX(YR3(I),0.d0)
  300 CONTINUE
      DO 310 I=NT+1,NF
        FX1(I)=(0.d0,0.d0)
        FX2(I)=(0.d0,0.d0)
        FX3(I)=(0.d0,0.d0)
        FY1(I)=(0.d0,0.d0)
        FY2(I)=(0.d0,0.d0)
        FY3(I)=(0.d0,0.d0)
  310 CONTINUE
C
      CALL WVFAST(NF,FX1,MWV,-1)
      CALL WVFAST(NF,FX2,MWV,-1)
      CALL WVFAST(NF,FX3,MWV,-1)
      CALL WVFAST(NF,FY1,MWV,-1)
      CALL WVFAST(NF,FY2,MWV,-1)
      CALL WVFAST(NF,FY3,MWV,-1)
      DUR=dFLOAT(NF)*DT
      DO 400 I=1,NF
        YX1(I)=FY1(I)/FX1(I)
        YX2(I)=FY2(I)/FX2(I)
        YX3(I)=FY3(I)/FX3(I)
  400 CONTINUE
C
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 * Output Fourier Amplitude Ratio data *
c
      write(*,*)
      write(*,*)'Fourier Amplitude Ratios are stored in "ARatio.csv"'
      open(10,file='ARatio.csv',status='unknown')
c      write(10,*)'frequency(Hz),',head2,',',head3,',',head4
      write(10,*)'frequency(Hz),comp1, comp2, comp3'
c      DO 680 I = 2,NOM
      DO 680 I = 1,NOM
        f=dFLOAT(I-1)*DOM/6.2831853070d0
        write(10,1900)f,cdabs(YX1(I)),cdabs(YX2(I)),cdabs(YX3(I))
  680 CONTINUE
 1900 format(f8.4,6(',',e12.4))
      close(10)
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
