C
      PROGRAM tran1x
C
      PARAMETER(MWV=68000,MMWV=40)
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER INFILE*80,head_in(40)*10,head_in1*9
      CHARACTER commentFin*10
      dimension time(MWV)
      dimension XR(MWV,MMWV),X(MWV,MMWV)
      COMPLEX*16 FX(MWV),FX_amp(MWV,MMWV)
      real MEANxr(MMWV),ISUMxr(MMWV),fin_com,fin_com2
C
C  Prompt for input RECON'ed file
C
      WRITE(*,1010) 
 1010 FORMAT(1X,/' Enter input file of waveform data')
      READ(*,1020) INFILE
 1020 FORMAT(A80)
C      
C *** Read input data file  ***
C
      OPEN(11,FILE=INFILE)
        read(11,*)head_in1,(head_in(j),j=1,40)
        i=1
   10   read(11,*,end=11)time(i),(xr(i,j),j=1,40)
        i=i+1
        goto 10
   11 CLOSE(11)
c
      NT=i-1
      dt=time(2)-time(1)
      write(*,*)
      write(*,*)'NT and DT=',NT,DT
c
C *** change of + or -  ***
C
      WRITE(*,*)
      WRITE(*,*)' Change +/-? (YES=1)?? :'
      READ(*,*)MYES
C 
      IF(MYES.EQ.1)THEN
C
C        WRITE(*,*)
C        WRITE(*,*)' Change from + to - of FILE 1 ? (YES=1)?? :'
C        READ(*,*)MY1
C
C         IF(MY1.EQ.1)THEN
           DO 100 I=1,NT
		   do 101 j=1,40
             XR(I,j)=-XR(I,j)
  101      CONTINUE
  100      CONTINUE
      END IF
C
C        WRITE(*,*)' Change from + to - of FILE 2 ? (YES=1)?? :'
C        READ(*,*)MY2
C
C         IF(MY2.EQ.1)THEN
C           DO 110 I=1,NT
C             XR2(I)=-XR2(I)
C  110      CONTINUE
C         END IF
C
C        WRITE(*,*)' Change from + to - of FILE 3 ? (YES=1)?? :'
C        READ(*,*)MY3
C
C         IF(MY3.EQ.1)THEN
C           DO 120 I=1,NT
C             XR3(I)=-XR3(I)
C  120      CONTINUE
C         END IF
C
C      END IF

       Do 211 j=1,40 
	 ISUMxr(j) = 0
        DO 201 i=1,nt
		ISUMxr(j) = ISUMxr(j) + xr(i,j)

  201	  CONTINUE

       MEANxr(j) =real(ISUMxr(j)) / real(nt)
  211  CONTINUE
C
C       ISUMxr2 = 0
C
C       DO 201,i=1,nt
C
C       ISUMxr2 = ISUMxr2 + xr2(i)
C
C  201  CONTINUE
C
C       MEANxr2 = real(ISUMxr2)/real(nt)
C
C       ISUMxr3 = 0
C
C       DO 301,i=1,nt
C
C       ISUMxr3 = ISUMxr3 + xr3(i)
C
C  301  CONTINUE
C
C       MEANxr3 = real(ISUMxr3) /real( nt)
C      
      write(*,*)
      write(*,*)'入力波形の軸の補正を行いますか(YES=1,NO=0)?? :'
       READ(*,*)tr
       IF(tr.EQ.1)THEN
	  Do 56 j=1,40
           DO 55 I=1,NT
             XR(I,j)=XR(I,j)- MEANxr(j)
C             XR2(I)=XR2(I)- MEANxr2
C             XR3(I)=XR3(I)- MEANxr3
  55      CONTINUE
  56    CONTINUE
       END IF
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(*,*)' IF1,NF=',IF1,NF
C*****************************************
C
      DUR=dFLOAT(NF)*DT
      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
      WRITE(*,*)
      WRITE(*,*)' Band Pass Filter? (YES=1)?? :'
      READ(*,*)NBD
C
      IF(NBD.EQ.1)THEN
C
        WRITE(*,*)
        WRITE(*,*)'  TRAPEZOIDAL BAND PATH FILTER'
        WRITE(*,*)'  ENTER fL1<fL2<fH1<fH2 (Hz)'
        READ(*,*)fL1,fL2,fH1,fH2
        OML1=fL1*6.2831853070d0
        OML2=fL2*6.2831853070d0
        OMH1=fH1*6.2831853070d0
        OMH2=fH2*6.2831853070d0
      END IF
C*****************************************
      WRITE(*,*)
      WRITE(*,*)' Integration/Differentiation:'
      WRITE(*,*)' NO Change         :                 (Enter=0)'
      WRITE(*,*)' Integration       : Ex.) acc -> vel (Enter=1)'
      WRITE(*,*)' Double Integration: Ex.) acc -> dis (Enter=2)'
      WRITE(*,*)' Differentiation:  : Ex.) Vel -> acc (Enter=3)'
      READ(*,*)NIT
C*****************************************
	DO 305 j=1,40
C
      DO 300 I=1,NT
        FX(I)=dCMPLX(XR(I,j),0.d0)
C        FX2(I)=dCMPLX(XR2(I),0.d0)
C        FX3(I)=dCMPLX(XR3(I),0.d0)
  300 CONTINUE
      DO 310 I=NT+1,NF
        FX(I)=(0.d0,0.d0)
C        FX2(I)=(0.d0,0.d0)
C        FX3(I)=(0.d0,0.d0)
  310 CONTINUE
C
      CALL WVFAST(NF,FX,MWV,-1)
C      CALL WVFAST(NF,FX2,MWV,-1)
C      CALL WVFAST(NF,FX3,MWV,-1)
C      DUR=dFLOAT(NF)*DT
      DO 400 I=1,NF
        FX(I)=FX(I)*DT
C        FX2(I)=FX2(I)*DT
C        FX3(I)=FX3(I)*DT
  400 CONTINUE
C
c      FOM=0.d0
C     
      IF(NBD.EQ.1)THEN
        DO 500 I = 1,NOM
           OMEGA=dFLOAT(I-1)*DOM
           IF(OMEGA.LE.OML1)then
             FX(I)=(0.d0,0.d0)
C             FX2(I)=(0.d0,0.d0)
C             FX3(I)=(0.d0,0.d0)
           end if
           IF(OMEGA.GE.OMH2)then
             FX(I)=(0.d0,0.d0)
C             FX2(I)=(0.d0,0.d0)
C             FX3(I)=(0.d0,0.d0)
           end if
           AA=1.d0
           IF(OMEGA.LT.OML2)AA=(OMEGA-OML1)/(OML2-OML1)
           IF(OMEGA.GT.OMH1)AA=(OMH2-OMEGA)/(OMH2-OMH1)
	   FX(I) = AA*FX(I)
C	   FX2(I) = AA*FX2(I)
C	   FX3(I) = AA*FX3(I)
  500   CONTINUE
C
      END IF
C
c
      IF(NIT.EQ.1)THEN
        DO 550 I = 2,NOM
           OMEGA=dFLOAT(I-1)*DOM
	   FX(I) = FX(I)/dcmplx(0.d0,omega)
C	   FX2(I) = FX2(I)/dcmplx(0.d0,omega)
C	   FX3(I) = FX3(I)/dcmplx(0.d0,omega)
  550   CONTINUE
      END IF
c
      IF(NIT.EQ.2)THEN
        DO 560 I = 2,NOM
           OMEGA=dFLOAT(I-1)*DOM
	   FX(I) = FX(I)/(-omega*omega)
C	   FX2(I) = FX2(I)/(-omega*omega)
C	   FX3(I) = FX3(I)/(-omega*omega)
  560   CONTINUE
      END IF
c
      IF(NIT.EQ.3)THEN
        DO 570 I = 2,NOM
           OMEGA=dFLOAT(I-1)*DOM
	   FX(I) = FX(I)*dcmplx(0.d0,omega)
C	   FX2(I) = FX2(I)*dcmplx(0.d0,omega)
C	   FX3(I) = FX3(I)*dcmplx(0.d0,omega)
  570   CONTINUE
      END IF
C
c * Output Fourier Spectrum data *
c
C      write(*,*)
C      write(*,*)'Fourier Spectrum data are stored in "amp3.csv"'
C      open(10,file='amp3.csv',status='unknown')
C      write(10,*)'frequency(Hz),',head2,',',head3,',',head4
      DO 680 I = 2,NOM
C        f=dFLOAT(I-1)*DOM/6.2831853070d0
        FX_amp(I,j)=FX(I)
  680 CONTINUE
C 1900 format(f15.8,6(',',e15.8))
C      close(10)
C
C *** INVERSE FOURIER TRANSFORM  ***
C
      DO 600 I=2,NF/2
           FX(NF+2-I)=dCMPLX(dREAL(FX(I)),-dIMAG(FX(I)))
C           FX2(NF+2-I)=dCMPLX(dREAL(FX2(I)),-dIMAG(FX2(I)))
C           FX3(NF+2-I)=dCMPLX(dREAL(FX3(I)),-dIMAG(FX3(I)))
  600 CONTINUE
C
      CALL WVFAST(NF,FX,MWV,1)
C      CALL WVFAST(NF,FX2,MWV,1)
C      CALL WVFAST(NF,FX3,MWV,1)
C
C *** WAVE DATA ***
C
      DO 700 I = 1,NF
         XR(I,j)=dREAL(FX(I))/DUR
C         XR2(I)=dREAL(FX2(I))/DUR
C         XR3(I)=dREAL(FX3(I))/DUR
  700 CONTINUE
  305 CONTINUE

C
      write(*,*)
      write(*,*)'Fourier Spectrum data are stored in "amp3.csv"'
      open(10,file='amp3.csv',status='unknown')
      write(10,1901)'frequency(Hz)',(head_in(j),j=1,40)
 1901 format(a13,40(',',a10))
      DO 681 I = 2,NOM
        f=dFLOAT(I-1)*DOM/6.2831853070d0
        write(10,1900)f,(cdabs(FX_amp(I,j)),j=1,40)
  681 CONTINUE
 1900 format(f15.8,40(',',e15.8))
      close(10)
c
c *** Cut ***
c
      do 304 i=1,nt
	 do 306 j=1,40

        x(i,j)=xr(i,j)
C        x2(i)=xr2(i)
C        x3(i)=xr3(i)
  306 end do
  304 end do

      write(*,*)
      write(*,*)'Cutting data? (No=0,Yes=1)'
      read(*,*)icut
      if(icut.eq.1)then
        write(*,*)'Enter Starting Time (sec)'
        read(*,*)stime
        nstime=int(stime/dt)
        write(*,*)'Number of Cutting Data =',nstime
        write(*,*)'Enter Duration (sec)'
        read(*,*)dur
        ndur=int(dur/dt)+1
        write(*,*)'Number of data =',ndur
	 do 307 j=1,40
        do i=1,ndur
          if((i+nstime).le.0)then
            x(i,j)=0.d0
C            x2(i)=0.d0
C            x3(i)=0.d0
	  else
            x(i,j)=xr(i+nstime,j)
C            x2(i)=xr2(i+nstime)
C            x3(i)=xr3(i+nstime)
          end if
        end do
  307	 Continue
	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(*,2010) 
 2010 FORMAT(/1X,'Enter output file name of waveform (csv format file)')
      READ(*,2020) INFILE
 2020 FORMAT(A30)
      OPEN(11,FILE=INFILE,STATUS='unknown')
      write(11,2001)head_in1,(head_in(j),j=1,40)
 2001 format(a7,40(',',a10))
      do i=1,nt,iskip
      write(11,2000)dt*(i-1),(X(I,j),j=1,40)
c      write(11,2000)dt*(i-1),X1(I)-X1(1),x2(i)-X2(1),x3(i)-x3(1)
      end do
 2000 format(f10.5,40(',',e15.8))
      CLOSE(11)
	 fin_com=rand()*1000000/3
C	 fin_com= mod(int(fin_com2),3)
	 commentFin="Done"
	if(fin_com.GE.0.and.fin_com.LT.3)then
	  commentFin="完了"
      end if
	if (fin_com.GE.3. And .fin_com.LT.6)then
	  commentFin="官僚"
      end if
	if (fin_com.GE.6. and .fin_com.LT.9) then
	  commentFin="かんりょう"
      end if
	if (fin_com.GE.9) then
	  commentFin="感無量"
      end if

      write(*,*)commentFin,fin_com
C      read(*,*)
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
