C
      PROGRAM mwvtran
C
      PARAMETER(MWV=5000,mob=100)
      CHARACTER INFILE*20
      REAL TIME(MWV),XR1(MWV),XR2(MWV),XR3(MWV)
      REAL R1(MWV,mob),R2(MWV,mob),R3(MWV,mob)
      COMPLEX*8 FX1(MWV),FX2(MWV),FX3(MWV)
C
      WRITE(*,*)
      WRITE(*,*)
      WRITE(*,*) 'PROGRAM MWVPLOT - To obtain filtered seismograms'
      WRITE(*,*) '------------------------------------------------------
     +-----------'
      WRITE(*,*)
C
C  Prompt for input RECON'ed file
C
      WRITE(*,1000) 
 1000 FORMAT(1X,'Enter name of the wave datafile')
      READ(*,1020) INFILE
 1020 format(a20)
C      
C *** Read input data file  ***

      OPEN(11,FILE=INFILE,STATUS='OLD')
        read(11,*)nob
        write(*,*)
        write(*,*)'Number of Observation Points =',nob
        if(nob.gt.mob)then
          write(*,*)'nob>mob !!!'
          stop
        end if

        i=1
 1      do j=1,nob
          READ(11,*,end=2)time(i),R1(I,j),R2(I,j),R3(I,j)
        end do
        i=i+1
        goto 1
 2    CLOSE(11)
c
      dt=time(2)-time(1)
      nt=i-1
      write(*,*)'dt,nt=',dt,nt
      if(nt.gt.mwv)then
        write(*,*)'nt>mwv !!!'
        stop
      end if
C
      NF=2
      DO 200 I=2,14
        NF=2*NF
        IF(NF.GE.NT)GOTO 10
  200 CONTINUE
      WRITE(6,*)' NT>2**14 !!!!!'
      STOP
   10 IFF=I
      WRITE(6,*)
      WRITE(6,*)' IFF,NF=',IFF,NF

      WRITE(6,*)
      WRITE(6,*)' Change +/-? (YES=1)?? :'
      READ(5,*)MYES
c
      IF(MYES.EQ.1)THEN
        WRITE(6,*)
        WRITE(6,*)' Change from + to - of FILE 1 ? (YES=1)?? :'
        READ(5,*)MY1
        WRITE(*,*)' Change from + to - of FILE 2 ? (YES=1)?? :'
        READ(5,*)MY2
        WRITE(*,*)' Change from + to - of FILE 3 ? (YES=1)?? :'
        READ(5,*)MY3
      END IF
c
      WRITE(*,*)
      WRITE(*,*)' BAND PASS FILTER? (YES=1)?? :'
      READ(5,*)NBD
C
      IF(NBD.EQ.1)THEN
        WRITE(6,*)
        WRITE(6,*)'  TRAPEZOIDAL BAND PATH FILTER'
        WRITE(6,*)'  ENTER OML1<OML2<OMH1<OMH2'
        READ(5,*)OML1,OML2,OMH1,OMH2
      end if
C
      WRITE(*,*)
      WRITE(*,*)' Integration(=1) or Differentiation(=2)? :'
      READ(5,*)NIT
C 
c *** Big-Loop for NOB ***
c
      do 900 iob=1,nob

      do i=1,nt
        xr1(i)=r1(i,iob)
        xr2(i)=r2(i,iob)
        xr3(i)=r3(i,iob)
      end do
c
      IF(MYES.EQ.1)THEN
         IF(MY1.EQ.1)THEN
           DO 100 I=1,NT
             XR1(I)=-XR1(I)
  100      CONTINUE
         END IF
         IF(MY2.EQ.1)THEN
           DO 110 I=1,NT
             XR2(I)=-XR2(I)
  110      CONTINUE
         END IF
         IF(MY3.EQ.1)THEN
           DO 120 I=1,NT
             XR3(I)=-XR3(I)
  120      CONTINUE
         END IF
      END IF
C
C *** FOURIER TRANSFORM ***
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
      FOM=0.0
      DOM=6.2831853070/DUR
      NOM=NF/2+1
      EOM=DOM*(NOM-1)
      if(nob.eq.1)then
        WRITE(6,*)
        WRITE(6,*)' DUR,DOM,NOM,EOM=',DUR,DOM,NOM,EOM
      end if
C     
C *** BAND PASS FILTER ***
C
      IF(NBD.EQ.1)THEN
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
c *** Integration or Differentiation
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.1)THEN
        DO 570 I = 1,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 *** 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
         R1(I,iob)=REAL(FX1(I))/DUR
         R2(I,iob)=REAL(FX2(I))/DUR
         R3(I,iob)=REAL(FX3(I))/DUR
  700 CONTINUE
c
 900  continue
C
c *** output ***
c
      WRITE(6,*)'Output wave data to wavem_dat:' 
      OPEN(11,FILE='wavem_dat',STATUS='unknown')
      write(11,*)nob,nob,nob,nob
      do i=1,nt
      do j=1,nob
      write(11,2100)time(i),R1(I,j),R2(I,j),R3(I,j)
      end do
      end do
 2100 format(f8.3,3e12.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.1415926540
      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
