C
      PROGRAM tranmex
C
      PARAMETER(MWV=30000)
      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(*,*)
      WRITE(*,*) 'Read from "wave1_dat" '
      WRITE(*,*) '---------------------'
      WRITE(*,*)
C
C *** Read input data file  ***
C
      OPEN(11,FILE='wave1_dat',status='unknown')
        i=1
   10   read(11,*,end=11)time(i),xr1(i),xr2(i),xr3(i)
        i=i+1
        goto 10
   11 CLOSE(11)
c
      DT=time(2)-time(1)
      NT=i-1
      write(*,*)'NT,DT=',NT,DT
      dt1=DT
      dt2=DT
      dt3=DT
      nt1=nt
      nt2=nt
      nt3=nt
c
C *** change of + or -  ***
C
      WRITE(*,*)
      WRITE(*,*)' Change +/-? (YES=1)?? :'
      READ(*,*)MYES
C 
      IF(MYES.EQ.1)THEN
C
        WRITE(*,*)
        WRITE(*,*)' Change from + to - of FILE 1 ? (YES=1)?? :'
        READ(*,*)MY1
C
         IF(MY1.EQ.1)THEN
           DO 100 I=1,NT1
             XR1(I)=-XR1(I)
  100      CONTINUE
         END IF
C
        WRITE(*,*)' Change from + to - of FILE 2 ? (YES=1)?? :'
        READ(*,*)MY2
C
         IF(MY2.EQ.1)THEN
           DO 110 I=1,NT2
             XR2(I)=-XR2(I)
  110      CONTINUE
         END IF
C
        WRITE(*,*)' Change from + to - of FILE 3 ? (YES=1)?? :'
        READ(*,*)MY3
C
         IF(MY3.EQ.1)THEN
           DO 120 I=1,NT3
             XR3(I)=-XR3(I)
  120      CONTINUE
         END IF
C
      END IF
C
C *** FOURIER TRANSFORM ***
C
      NF1=2
      DO 200 I=2,15
        NF1=2*NF1
        IF(NF1.GE.NT1)GOTO 30
  200 CONTINUE
      WRITE(*,*)' NT1>2**15 !!!!!'
      STOP
   30 IF1=I
      WRITE(*,*)
      WRITE(*,*)' IF1,NF1=',IF1,NF1
C
      NF2=2
      DO 210 I=2,15
        NF2=2*NF2
        IF(NF2.GE.NT2)GOTO 31
  210 CONTINUE
      WRITE(*,*)' NT2>2**15 !!!!!'
      STOP
   31 IF2=I
      WRITE(*,*)' IF2,NF2=',IF2,NF2
C
      NF3=2
      DO 220 I=2,15
        NF3=2*NF3
        IF(NF3.GE.NT3)GOTO 32
  220 CONTINUE
      WRITE(*,*)' NT3>2**15 !!!!!'
      STOP
   32 IF3=I
      WRITE(*,*)' IF3,NF3=',IF3,NF3
C
      DO 300 I=1,NT1
        FX1(I)=CMPLX(XR1(I),0.0)
  300 CONTINUE
      DO 310 I=NT1+1,NF1
        FX1(I)=(0.0,0.0)
  310 CONTINUE
C
      DO 320 I=1,NT2
        FX2(I)=CMPLX(XR2(I),0.0)
  320 CONTINUE
      DO 330 I=NT2+1,NF2
        FX2(I)=(0.0,0.0)
  330 CONTINUE
C
      DO 340 I=1,NT3
        FX3(I)=CMPLX(XR3(I),0.0)
  340 CONTINUE
      DO 350 I=NT3+1,NF3
        FX3(I)=(0.0,0.0)
  350 CONTINUE
C
      CALL WVFAST(NF1,FX1,MWV,-1)
      DUR1=FLOAT(NF1)*DT1
      DO 400 I=1,NF1
        FX1(I)=FX1(I)*DT1
  400 CONTINUE
C
c      FOM1=0.0
      DOM1=6.2831853070e0/DUR1
      NOM1=NF1/2+1
      EOM1=DOM1*NOM1
      WRITE(*,*)
      WRITE(*,*)' DUR1,DOM1,NOM1,EOM1=',DUR1,DOM1,NOM1,EOM1
C
      CALL WVFAST(NF2,FX2,MWV,-1)
      DUR2=FLOAT(NF2)*DT2
      DO 410 I=1,NF2
        FX2(I)=FX2(I)*DT2
  410 CONTINUE
C
c      FOM2=0.0
      DOM2=6.2831853070e0/DUR2
      NOM2=NF2/2+1
      EOM2=DOM2*NOM2
      WRITE(*,*)' DUR2,DOM2,NOM2,EOM2=',DUR2,DOM2,NOM2,EOM2
C
      CALL WVFAST(NF3,FX3,MWV,-1)
      DUR3=FLOAT(NF3)*DT3
      DO 420 I=1,NF3
        FX3(I)=FX3(I)*DT3
  420 CONTINUE
C
c      FOM3=0.0
      DOM3=6.2831853070e0/DUR3
      NOM3=NF3/2+1
      EOM3=DOM3*NOM3
      WRITE(*,*)' DUR3,DOM3,NOM3,EOM3=',DUR3,DOM3,NOM3,EOM3
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,NOM1
           OMEGA=FLOAT(I-1)*DOM1
           IF(OMEGA.LE.OML1)FX1(I)=(0.0,0.0)
           IF(OMEGA.GE.OMH2)FX1(I)=(0.0,0.0)
           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)
  500   CONTINUE
C
        DO 510 I = 1,NOM2
           OMEGA=FLOAT(I-1)*DOM2
           IF(OMEGA.LE.OML1)FX2(I)=(0.0,0.0)
           IF(OMEGA.GE.OMH2)FX2(I)=(0.0,0.0)
           AA=1.0
           IF(OMEGA.LT.OML2)AA=(OMEGA-OML1)/(OML2-OML1)
           IF(OMEGA.GT.OMH1)AA=(OMH2-OMEGA)/(OMH2-OMH1)
	   FX2(I) = AA*FX2(I)
  510   CONTINUE
C
        DO 520 I = 1,NOM3
           OMEGA=FLOAT(I-1)*DOM3
           IF(OMEGA.LE.OML1)FX3(I)=(0.0,0.0)
           IF(OMEGA.GE.OMH2)FX3(I)=(0.0,0.0)
           AA=1.0
           IF(OMEGA.LT.OML2)AA=(OMEGA-OML1)/(OML2-OML1)
           IF(OMEGA.GT.OMH1)AA=(OMH2-OMEGA)/(OMH2-OMH1)
	   FX3(I) = AA*FX3(I)
  520   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,NOM1
           OMEGA=FLOAT(I-1)*DOM1
	   FX1(I) = FX1(I)/cmplx(0.0,omega)
  550   CONTINUE
        DO 552 I = 2,NOM2
           OMEGA=FLOAT(I-1)*DOM2
	   FX2(I) = FX2(I)/cmplx(0.0,omega)
  552   CONTINUE
        DO 554 I = 2,NOM3
           OMEGA=FLOAT(I-1)*DOM3
	   FX3(I) = FX3(I)/cmplx(0.0,omega)
  554   CONTINUE
      END IF
c
      IF(NIT.EQ.2)THEN
        DO 560 I = 2,NOM1
           OMEGA=FLOAT(I-1)*DOM1
	   FX1(I) = FX1(I)/(-omega*omega)
  560   CONTINUE
        DO 562 I = 2,NOM2
           OMEGA=FLOAT(I-1)*DOM2
	   FX2(I) = FX2(I)/(-omega*omega)
  562   CONTINUE
        DO 564 I = 2,NOM3
           OMEGA=FLOAT(I-1)*DOM3
	   FX3(I) = FX3(I)/(-omega*omega)
  564   CONTINUE
      END IF
c
      IF(NIT.EQ.3)THEN
        DO 570 I = 2,NOM1
           OMEGA=FLOAT(I-1)*DOM1
	   FX1(I) = FX1(I)*cmplx(0.0,omega)
  570   CONTINUE
        DO 572 I = 2,NOM2
           OMEGA=FLOAT(I-1)*DOM2
	   FX2(I) = FX2(I)*cmplx(0.0,omega)
  572   CONTINUE
        DO 574 I = 2,NOM3
           OMEGA=FLOAT(I-1)*DOM3
	   FX3(I) = FX3(I)*cmplx(0.0,omega)
  574   CONTINUE
      END IF
C
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,NOM1
        OMEGA=FLOAT(I-1)*DOM1
        write(10,1900)omega,FX1(I),FX2(I),FX3(I)
  680 CONTINUE
 1900 format(f8.4,6e12.4)
      close(10)
C
C *** INVERSE FOURIER TRANSFORM  ***
C
      DO 600 I=2,NF1/2
           FX1(NF1+2-I)=CMPLX(REAL(FX1(I)),-AIMAG(FX1(I)))
  600 CONTINUE
C
      DO 610 I=2,NF2/2
           FX2(NF2+2-I)=CMPLX(REAL(FX2(I)),-AIMAG(FX2(I)))
  610 CONTINUE
C
      DO 620 I=2,NF3/2
           FX3(NF3+2-I)=CMPLX(REAL(FX3(I)),-AIMAG(FX3(I)))
  620 CONTINUE
C
      CALL WVFAST(NF1,FX1,MWV,1)
      CALL WVFAST(NF2,FX2,MWV,1)
      CALL WVFAST(NF3,FX3,MWV,1)
C
C *** WAVE DATA ***
C
      DO 700 I = 1,NF1
         time(I) = FLOAT(I-1)*DT1
         XR1(I)=REAL(FX1(I))/DUR1
  700 CONTINUE
C
      DO 710 I = 1,NF2
         time(I) = FLOAT(I-1)*DT2
         XR2(I)=REAL(FX2(I))/DUR2
  710 CONTINUE
C
      DO 720 I = 1,NF3
         time(I) = FLOAT(I-1)*DT3
         XR3(I)=REAL(FX3(I))/DUR3
  720 CONTINUE
C
C *** MAX AND MIN AMPLITUDE ***
C
      XMIN1=0.0 
      XMIN2=0.0 
      XMIN3=0.0 
      XMAX1=0.0 
      XMAX2=0.0 
      XMAX3=0.0 
C
      DO 800 I=1,NF1
        IF(XR1(I).GT.XMAX1)XMAX1=XR1(I)
        IF(XR1(I).LT.XMIN1)XMIN1=XR1(I)
  800 CONTINUE
C
      DO 810 I=1,NF2
        IF(XR2(I).GT.XMAX2)XMAX2=XR2(I)
        IF(XR2(I).LT.XMIN2)XMIN2=XR2(I)
  810 CONTINUE
C
      DO 820 I=1,NF3
        IF(XR3(I).GT.XMAX3)XMAX3=XR3(I)
        IF(XR3(I).LT.XMIN3)XMIN3=XR3(I)
  820 CONTINUE
C
      AMAX=XMAX1
      AMIN=XMIN1
      IF(XMAX2.GT.AMAX)AMAX=XMAX2
      IF(XMAX3.GT.AMAX)AMAX=XMAX3
      IF(XMIN2.LT.AMIN)AMIN=XMIN2
      IF(XMIN3.LT.AMIN)AMIN=XMIN3
C
      AAMAX=AMAX
      IF(ABS(AMIN).GT.AAMAX)AAMAX=ABS(AMIN)
C      
      TMIN1O=0.0
      TMAX1O=time(NF1)
      TMIN1 =-0.05*TMAX1O
      TMAX1 = 1.05*TMAX1O
      WRITE(*,*)
      WRITE(*,*)'ORIGINAL TMIN1,TMAX1 =',TMIN1O,TMAX1O     
      WRITE(*,*)'OUTPUT   TMIN1,TMAX1 =',TMIN1,TMAX1     
C      
      TMIN2O=0.0
      TMAX2O=time(NF2)
      TMIN2 =-0.05*TMAX2O
      TMAX2 = 1.05*TMAX2O
      WRITE(*,*)'ORIGINAL TMIN2,TMAX2 =',TMIN2O,TMAX2O     
      WRITE(*,*)'OUTPUT   TMIN2,TMAX2 =',TMIN2,TMAX2     
C      
      TMIN3O=0.0
      TMAX3O=time(NF3)
      TMIN3 =-0.05*TMAX3O
      TMAX3 = 1.05*TMAX3O
      WRITE(*,*)'ORIGINAL TMIN3,TMAX3 =',TMIN3O,TMAX3O     
      WRITE(*,*)'OUTPUT   TMIN3,TMAX3 =',TMIN3,TMAX3     
C
      WRITE(*,*)
      WRITE(*,*)'XMIN1,XMIN2,XMIN3 =',XMIN1,XMIN2,XMIN3
      WRITE(*,*)'XMAX1,XMAX2,XMAX3 =',XMAX1,XMAX2,XMAX3
      WRITE(*,*)'AMIN,AMAX,AAMAX=',AMIN,AMAX,AAMAX
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(*,*)'wave1_dat:' 
      OPEN(11,FILE='wave1_dat',STATUS='unknown')
      do i=1,nt,iskip
      write(11,2000)dt*(i-1),X1(I),x2(i),x3(i)
      end do
 2000 format(4e12.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
      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
