C 
C *** Wave Simulaion using Amplification Factor	of Layered Structure
c     using the Haskell's method
C
      PROGRAM mwave
C
C INPUT DATA  
C ** FOM    : FIRST CIRCULAR FREQUENCY ( OMEGA )
C ** DOM    : INCREMENT OF CIRCULAR FREQUENCY ( DELTA OMEGA )
C ** NOM    : TOTAL NUMBER OF CIRCULAR FREQUENCY
C ** MOM    : NOM IN MAIN PROGRAM DIMENSION
C ** DNS(I) : DENSITY OF Ith LAYER
C ** VS (I) : VS OF Ith LAYER
C ** THK(I) : THICKNESS OF Ith LAYER
C ** VE (I) : VE=sqrt(E/dns) OF Ith LAYER (E: Young Coefficient)
C ** NL     : TOTAL NUMBER OF LAYERS
C ** ML     : NL IN MAIN PROGRAM DIMENSION
C ** THETA  : INCIDENT ANGLE OF SH WAVE ( RADIAN )
C OUTPUT DATA  
C ** AMP(I) : Amplifivarion FACTER ON SURFACE OF Ith OMEGA ( (I-1)*DOM+FOM )
C ** SPCA(I): SPECTRUM AMPLITUDE OF AMP(I)
C
      PARAMETER(ML=50,MOM=35800)
      DIMENSION DNS(ML),VP(ML),QP(ML),VS(ML),QS(ML),THK(ML)
      COMPLEX AMP(MOM)
      CHARACTER INFILE1*20
      REAL TIME(mom)
      REAL XR1(mom),XI1(mom)
      COMPLEX FX1(mom),delay
C
      WRITE(*,*)
      WRITE(*,*)' : Soil Data must be "soiln_dat"'
C
C  Prompt for wave file
C
      WRITE(*,1010) 
 1010 FORMAT(/1X,'Enter Acceleration Wave Data',)
      READ(*,1020) INFILE1
 1020 FORMAT(A20)
C
      OPEN(11,FILE=INFILE1,status='old')
        i=1
   10   read(11,*,end=11)time(i),xr1(i)
c        write(*,*)time(i),xr1(i)
        i=i+1
        goto 10
   11 CLOSE(11)
c
      NT=i-1
      dt=time(2)-time(1)
      write(*,*)'NT and DT=',NT,DT
c
C *** Find the Closest Number of 2**N ***
C
      NF=2
      do I=2,15
        NF=2*NF
        IF(NF.GE.NT)GOTO 30
      end do
      WRITE(*,*)' NT>2**15 !!!!!'
      STOP
   30 IF1=I
      WRITE(*,*)
      WRITE(*,*)' IF1,NF=',IF1,NF
c
C *** FOURIER TRANSFORM of the Original Wave ***
C
      do I=1,NT
        FX1(I)=CMPLX(xr1(I),0.0)
      end do
      do I=NT+1,NF
        FX1(I)=(0.0,0.0)
      end do
C
      CALL WVFAST(NF,FX1,mm,-1)
      DUR=dFLOAT(NF)*DT
      DOM=6.283185307/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
      do I=1,NF
        FX1(I)=FX1(I)*DT
      end do
c
c *** Calculation of Amplification Factor ***
c
      fom=0.0
      ndata=nom
      CALL SOILIN(DNS,VP,QP,VS,QS,THK,PQF,NL,ML,THETA,NSF)
      CALL MULTIF(FOM,DOM,NOM,MOM,DNS,VS,QS,THK,PQF,NL,ML,THETA,NSF,AMP)
      write(*,*)
      write(*,*)'S_amp: Amplification Factor for S wave'
      open(9,file='S_amp',status='unknown')
      do i=1,ndata
        om=dom*(i-1)
        write(9,1140)om/6.2831853,amp(i),cabs(amp(i))
        delay=(1.0,0.0)
        do l=1,nl-1
          delay=delay*cexp(cmplx(0.0,om*thk(l)/vs(l)))
        end do
        fx1(i)=fx1(i)*amp(i)*delay
        xr1(i)=real(fx1(i))
        xi1(i)=aimag(fx1(i))
      end do
      close(9)
 1140 format(4f8.4)
C     
C * DIFFERENTIATION *
C
      WRITE(*,*)
      WRITE(*,*)
     * ' You can integrate or differentiate the original wave.;'
      WRITE(*,*)
     * '   Do not change   (Enter -> 0),'
      WRITE(*,*)
     * '   Integrate       (Enter -> 1),'
      WRITE(*,*)
     * '   Differentiate   (Enter -> 2):'
      READ(*,*)NYES
C     
C * BAND PASS *
C
      WRITE(*,*)
      WRITE(*,*)' Band Pass Filter? (Yes=1 or No=0) :'
      READ(*,*)NBD
C
      OML1=0.0
      OMH2=99999.9
      IF(NBD.EQ.1)THEN
        WRITE(*,*)'   Data for TRAPEZOIDAL BAND PASS FILTER:'
        WRITE(*,*)'   INITIAL AND FINAL OMs =',0,DOM*(NDATA-1)
        WRITE(*,*)
        WRITE(*,*)
     *  '   Enter OML1,OML2,OMH1 and OMH2 (OML1<OML2<OMH1<OMH2): '
        READ(*,*)OML1,OML2,OMH1,OMH2
      END IF
C
C * Write data to a complex array *
C
      FX1(1) = CMPLX(0.0,0.0)
      DO 200 I = 2,NDATA
           OMEGA=FLOAT(I-1)*DOM
           AA=1.0
           IF(NBD.EQ.1)THEN
             IF((OMEGA.LE.OML1).or.(OMEGA.GE.OMH2))then
	       aa=0.0
             else
               IF(OMEGA.LT.OML2)AA=(OMEGA-OML1)/(OML2-OML1)
               IF(OMEGA.GT.OMH1)AA=(OMH2-OMEGA)/(OMH2-OMH1)
	     end if
           end if
    	   FX1(I) = AA*CMPLX(XR1(I),XI1(I))

           IF(NYES.EQ.1)THEN
             FX1(I)=FX1(I)/(OMEGA*CMPLX(0.0,1.0))
           END IF
           IF(NYES.EQ.2)THEN
             FX1(I)=FX1(I)*OMEGA*CMPLX(0.0,1.0)
           END IF
  200 CONTINUE
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 250 I = 2,NDATA
        OMEGA=FLOAT(I-1)*DOM
        freq=omega/6.283185307
        write(10,1500)freq,FX1(I),cabs(FX1(I))
  250 CONTINUE
 1500 format(f8.4,6e12.4)
      close(10)
C
C * INVERSE FOURIER TRANSFORM  *
C
      DO 400 I=2,NF/2
           FX1(NF+2-I)=CMPLX(REAL(FX1(I)),-AIMAG(FX1(I)))
  400 CONTINUE
C
      CALL WVFAST(NF,FX1,MOM,1)
C
C * WAVE DATA *
C
      DO 450 I = 1,NF
         TIME(I) = FLOAT(I-1)*DT
         XR1(I)=REAL(FX1(I))/DUR
  450 CONTINUE
C
C * MAX AND MIN AMPLITUDE *
C
      XMAX=0.0 
      XMIN=0.0 
C
      DO 550 I=1,NF
        IF(XR1(I).GT.XMAX)XMAX=XR1(I)
        IF(XR1(I).LT.XMIN)XMIN=XR1(I)
  550 CONTINUE
C
      AMAX=XMAX
      AMIN=XMIN
C
      WRITE(*,*)
      WRITE(*,*)' Nmber of data points    =',NF
      WRITE(*,*)' Duration and Delta Time =',DUR,DT
      WRITE(*,*)' Initial and final times =',0.0,DUR-DT
      WRITE(*,*)
      WRITE(*,*)' Min and Max amplitudes  =',AMIN,AMAX
C      
C * Output data file  *
C
      WRITE(*,*)
      WRITE(*,*)'Output data are stored in "mwave_dat"'
      WRITE(*,*)
      OPEN(11,FILE='mwave_dat',STATUS='unknown')
c      WRITE(11,*)NF,DT,AMIN,AMAX
      do i=1,NT
      write(11,2000)TIME(I),XR1(I)
      end do
      CLOSE(11)
 2000 FORMAT(F8.3,3E12.4)
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 POINTS (=2**I)
C * A(I): FOURIER COEFFICIENTS ( N/2+1 ) OR WAVE DATA (REAL PART)
C    ( FOURIER SPECTRUM = FOURIER COEFFICIENTS * DURATION/2.0 )
C * MOM : N IN MAIN DIMMENSION
C * IND : = -1 : TRANS.,   = 1 : INVERSE TRANS.
C
      COMPLEX A(MOM),TEMP,THETA
C
C * FFT and inverse 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
c
C
C *** READ OF SOIL DATA FROM DATA FILE
C 
      SUBROUTINE SOILIN(DNS,VP,QP,VS,QS,THK,PQF,NL,ML,THETA,NSF)
C
C ** THETA  : INCIDENT ANGLE OF SH WAVE ( DEGREE --> RADIAN )
C ** NSF    : =1:OUTPUT POINT IS ON SURFACE, =0:OUTPUT POINT IS IN SOIL
C
      DIMENSION DNS(ML),vp(ml),qp(ml),VS(ML),qs(ML),thk(ml)
C
C *** SOIL DATA INPUT ***
C 
      OPEN(4,FILE='soiln_dat',STATUS='OLD')
      READ(4,1010)THET1
      READ(4,1020)NSF
      READ(4,1030)NL,PQF
      DO 100 IL=1,NL
      READ(4,1040)DNS(IL),VP(IL),QP(IL),VS(IL),QS(IL),THK(IL)
  100 CONTINUE
 1010 FORMAT(F10.0)
 1020 FORMAT(I10)
 1030 FORMAT(I10,F10.0)
 1040 FORMAT(6f8.0)
C
C *** CHANGE THETA FROM DEGREE TO RADIAN ***
C
      THETA=THET1*0.017453292520e0
c
      RETURN
      END
c
C
C *** CUL. OF AMPLITUDE FACTER BY MULTIPLE REFLECTION METHOD OF S or P WAVE
C 
      SUBROUTINE MULTIF(FOM,DOM,NOM,MOM,DNS,VV,QF,THK,PQF
     *                 ,NL,ML,THETA,NSF,AMP)
C
C INPUT DATA  
C ** FOM    : FIRST CIRCULAR FREQUENCY ( OMEGA )
C ** DOM    : INCREMENT OF CIRCULAR FREQUENCY ( DELTA OMEGA )
C ** NOM    : TOTAL NUMBER OF CIRCULAR FREQUENCY
C ** MOM    : NOM IN MAIN PROGRAM DIMENSION
C           : IF MOM > 2**13(=8192), CAN NOT USE !!!
C ** DNS(I) : DENSITY OF Ith LAYER
C ** VV (I) : VS(=sqrt(G/dns)) or VE(=sqrt(E/dns)) OF Ith LAYER
C ** THK(I) : THICKNESS OF Ith LAYER
C ** QF(I)  : QUALITY FACTER OF Ith LAYER
C ** PQF    : Q=Qs*(freq)**PQF
C ** NL     : TOTAL NUMBER OF LAYERS
C ** ML     : NL IN MAIN PROGRAM DIMENSION
C ** THETA  : INCIDENT ANGLE OF SH WAVE ( RADIAN )
C ** NSF    : = 1:OUTPUT POINT IS ON SURFACE,  = 0:IN SOIL
C OUTPUT DATA  
C ** AMP(I) : AMPLITUDE FACTER ON SURFACE OF Ith OMEGA ( (I-1)*DOM+FOM )
C
      DIMENSION DNS(ML),THK(ML),VV(ML),QF(ML)
      COMPLEX AMP(MOM),VVN2,CINV2,R11,R12,R21,R22,VV2,SM,WN,WND,WNS
     *      ,COWN,SIWN,SA11,SA12,SA21,SA22,SB11,SB12,SB21,SB22,WNN,WNSN
C
C *** CAL. OF AMPLITUDE FACTER ***
C
      TSIN=SIN(THETA)
      SF=1.0
      IF(NSF.EQ.1)SF=2.0
C
      IFOM=1
      IF(FOM.EQ.0.0)THEN
      AMP(1)=CMPLX(SF,0.0)
      IFOM=2
      END IF
C
      DO 200 IO=IFOM,NOM
C
        OM=FLOAT(IO-1)*DOM+FOM
        freq=om/6.2831853
C
        R11=(1.0,0.0)
        R12=(0.0,0.0)
        R21=(0.0,0.0)
        R22=(1.0,0.0)

        VVN2=DCMPLX(VV(NL)*VV(NL),0.0)
        IF(PQF.ne.999.0)then
           Q=QF(NL)*(freq)**PQF
           VVN2=VVN2*CMPLX(1.0,1.0/Q)
        end if
        CINV2=TSIN*TSIN/VVN2
C
        DO 100 IL=NL-1,1,-1
C
          SA11=R11
          SA12=R12
          SA21=R21
          SA22=R22
C
          VV2=CMPLX(VV(IL)*VV(IL),0.0)
          IF(PQF.ne.999.0)then
            Q=QF(IL)*(freq)**PQF
            VV2=VV2*CMPLX(1.0,1.0/Q)
          end if
          SM=DNS(IL)*VV2   
          WN=OM*SQRT(1.0/VV2-CINV2)
          WND=WN*THK(IL)
          WNS=WN*SM
          COWN=COS(WND)
          SIWN=SIN(WND)
C
          SB11=COWN
          SB12=SIWN/WNS
          SB21=-WNS*SIWN
          SB22=COWN
C
          R11=SA11*SB11+SA12*SB21
          R12=SA11*SB12+SA12*SB22
          R21=SA21*SB11+SA22*SB21
          R22=SA21*SB12+SA22*SB22
C
  100   CONTINUE
C
        WNN=OM*SQRT(1.0/VVN2-CINV2)
        WNSN=WNN*DNS(NL)*VVN2
C
        AMP(IO)=SF/(R11-R21*CMPLX(0.0,1.0)/WNSN)
C
  200 CONTINUE
C
      RETURN
      END
c
