C 
C *** Wave Simulaion using Amplification Factor	for S and P wave
c     using the Haskell's method
C
      PROGRAM multi2
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=3000)
      DIMENSION DNS(ML),VP(ML),QP(ML),VS(ML),QS(ML),THK(ML),VE(ML)
      COMPLEX AMP(MOM)
      CHARACTER INFILE1*20
      REAL TIME(mom)
      REAL XR1(mom),XI1(mom),YR1(mom),YI1(mom),ZR1(mom),ZI1(mom)
      COMPLEX FX1(mom),FY1(mom),FZ1(mom)
C
      WRITE(*,*)
     *'"multi" - To obtain seismograms using grfltHF.f and'
      WRITE(*,*)'     the amplification factors of surface layers'
      WRITE(*,*)' : Soil Data must be "soil_dat"'
C
C * Read data *
c
      CALL SOILIN(DNS,VP,QP,VS,QS,THK,NQF,NL,ML,THETA,NSF)
c
      WRITE(*,*) 
      WRITE(*,1010) 
 1010 FORMAT(' Choose the input datafile (Enter Number 1 to 6):'/ 
     *' 1) grfault.dxyz (x,y,z-comp. of disp. from seismic sources)'/
     *' 2) grfault.dqvw (u,v,w-comp. of disp. from seismic sources)'/
     *' 3) grpoint.dxyz (x,y,z-comp. of disp. from point sources)'/
     *' 4) grpoint.sxyz (x,y,z-comp. of stress from point sources)'/
     *' 5) grpoint.dqvw (u,v,w-comp. of disp.  from point sources)'/
     *' 6) grpoint.sqvw (u,v,w-comp. of stress from point sources)')
      READ(*,*) IDF
C    
      IF(IDF.EQ.1)THEN
        INFILE1='grfault.dxyz'
      ELSE 
        IF(IDF.EQ.2)THEN
          INFILE1='grfault.dqvw'
        ELSE 
          IF(IDF.EQ.3)THEN
            INFILE1='grpoint.dxyz'
          ELSE
            IF(IDF.EQ.4)THEN
              INFILE1='grpoint.sxyz'
            ELSE 
              IF(IDF.EQ.5)THEN
                INFILE1='grpoint.dqvw'
              ELSE 
                IF(IDF.EQ.6)THEN
                  INFILE1='grpoint.sqvw'
                END IF
              END IF
            END IF
          END IF
        END IF
      END IF
c
      OPEN(11,FILE=INFILE1,STATUS='OLD')
c
      READ(11,1100)DOM,IOM,NDATA
 1100 FORMAT(1X,F10.0,2I10)
      WRITE(*,*)
      WRITE(*,*)' DOM,IOM,NDATA =',DOM,IOM,NDATA
      NDATA=NDATA+1
      WRITE(*,*)' New NDATA including OM=0 =',NDATA
      READ(11,1110)NOB
 1110 FORMAT(1X,I10///)
      WRITE(*,*)
      WRITE(*,*)' Total number of receivers =',NOB
      WRITE(*,1120) 
 1120 FORMAT(1X,' Enter the receiver number to output : ')
      READ(*,*)NOB1
      WRITE(*,*)
      WRITE(*,*)' XR,XI,YR,YI,ZR,ZI='
      DO 150 J=1,NOB1
      DO 150 I = IOM+1,NDATA
	 READ(11,1130) fx1(I),fy1(I),fz1(I)
         IF(J.NE.NOB1)GOTO 150
c	 WRITE(*,1130) fx1(I),fy1(I),fz1(I)
  150 CONTINUE
 1130 FORMAT(1X,6E12.5)
C
      CLOSE(11)
c
c *** Calculation of Amplification Factor ***
c
      fom=0.0
      nom=ndata
      CALL MULTIF(FOM,DOM,NOM,MOM,DNS,VS,QS,THK,NQF,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
        write(9,1140)dom*(i-1),amp(i)
        fx1(i)=fx1(i)*amp(i)
        fy1(i)=fy1(i)*amp(i)
        xr1(i)=real(fx1(i))
        xi1(i)=aimag(fx1(i))
        yr1(i)=real(fy1(i))
        yi1(i)=aimag(fy1(i))
      end do
      close(9)
 1140 format(3f8.4)
c
      do i=1,nl
        ram=dns(i)*(vp(i)*vp(i)-2.0*vs(i)*vs(i))
        sm=dns(i)*vs(i)*vs(i)
        E=(3.0*ram+2*sm)*sm/(ram+sm)
        VE(i)=sqrt(E/dns(i))
      end do
      CALL MULTIF(FOM,DOM,NOM,MOM,DNS,VE,QP,THK,NQF,NL,ML,THETA,NSF,AMP)
      write(*,*)
      write(*,*)'P_amp: Amplification Factor for P wave'
      open(9,file='P_amp',status='unknown')
      do i=1,ndata
        write(9,1140)dom*(i-1),amp(i)
        fz1(i)=fz1(i)*amp(i)
        zr1(i)=real(fz1(i))
        zi1(i)=aimag(fz1(i))
      end do
      close(9)
C     
C   * change of signs *
C
      WRITE(*,*)
      WRITE(*,*)' Change signs of amplitudes? (Yes=1 or No=0) :'
      READ(*,*)MYES
C 
      IF(MYES.EQ.1)THEN
C
        WRITE(*,*)'   Change signs of X? (Yes=1 or No=0) :'
        READ(*,*)MY1
C
         IF(MY1.EQ.1)THEN
           DO 180 I=1,NDATA
             XR1(I)=-XR1(I)
             XI1(I)=-XI1(I)
  180      CONTINUE
         END IF
C
        WRITE(*,*)'   Change signs of Y? (Yes=1 or No=0) :'
        READ(*,*)MY2
C
         IF(MY2.EQ.1)THEN
           DO 182 I=1,NDATA
             YR1(I)=-YR1(I)
             YI1(I)=-YI1(I)
  182      CONTINUE
         END IF
C
        WRITE(*,*)'   Change signs of Z? (Yes=1 or No=0) :'
        READ(*,*)MY3
C
         IF(MY3.EQ.1)THEN
           DO 184 I=1,NDATA
             ZR1(I)=-ZR1(I)
             ZI1(I)=-ZI1(I)
  184      CONTINUE
         END IF
C
      END IF
C     
C * DIFFERENTIATION *
C
      WRITE(*,*)
      WRITE(*,*)
     * ' You can transfer the displacements to vel. or acc.;'
      WRITE(*,*)
     * '   Do not transfer                              (Enter 0),'
      WRITE(*,*)
     * '   Transfer from displacements to velocities    (Enter 1),'
      WRITE(*,*)
     * '   Transfer from displacements to accelerations (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 =',DOM*IOM,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)
      FY1(1) = CMPLX(0.0,0.0)
      FZ1(1) = CMPLX(0.0,0.0)
      DO 200 I = 2,NDATA
           OMEGA=FLOAT(I-1)*DOM
           IF(OMEGA.LE.OML1)GOTO 200
           IF(OMEGA.GE.OMH2)GOTO 200
           AA=1.0
           IF(NBD.EQ.1)THEN
             IF(OMEGA.LT.OML2)AA=(OMEGA-OML1)/(OML2-OML1)
             IF(OMEGA.GT.OMH1)AA=(OMH2-OMEGA)/(OMH2-OMH1)
           END IF
	   FX1(I) = AA*CMPLX(XR1(I),XI1(I))
	   FY1(I) = AA*CMPLX(YR1(I),YI1(I))
	   FZ1(I) = AA*CMPLX(ZR1(I),ZI1(I))
           IF(NYES.EQ.1)THEN
             FX1(I)=FX1(I)*OMEGA*(0.0,1.0)
             FY1(I)=FY1(I)*OMEGA*(0.0,1.0)
             FZ1(I)=FZ1(I)*OMEGA*(0.0,1.0)
           END IF
           IF(NYES.EQ.2)THEN
             FX1(I)=-FX1(I)*OMEGA*OMEGA
             FY1(I)=-FY1(I)*OMEGA*OMEGA
             FZ1(I)=-FZ1(I)*OMEGA*OMEGA
           END IF
  200 CONTINUE
c
c * Output Fourier Spectrum data *
c
      write(*,*)
      write(*,*)'The Fourier Spectrum data are stored in "amp_dat"'
      write(*,*)'The amplitude of "amp_dat" is *100'
      open(10,file='amp_dat',status='unknown')
      DO 250 I = 2,NDATA
        OMEGA=FLOAT(I-1)*DOM
        write(10,1500)omega,FX1(I)*100.,FY1(I)*100.,FZ1(I)*100.
  250 CONTINUE
 1500 format(f8.4,6e12.4)
      close(10)
C
C * Find closest (but larger) 2**J value to the input data dimension *
C
      WRITE(*,*)
      DO 300 J = 1,999
           N = 2**J
           IF(NDATA.LE.N)THEN
                NUM = N
                GO TO 20
           END IF
  300 CONTINUE
C
   20 IF(NUM.GT.MOM)THEN
        WRITE(*,*)' NUM > MOM !!!'
        STOP
      END IF

      WRITE(*,*)' N and J of N=2**J =',NUM,J
      WRITE(*,*)' Increase J for smoother seismograms? (Yes=1 or No=0)'
      READ(*,*)IYES
      IF(IYES.EQ.1)THEN
        J=J+1
        NUM=NUM*2
        GO TO 20
      END IF
C
C * Pad the remaining spots (up to N=2**J) of the input array with zeroes *
C
      IF(NUM*2.GT.MOM)THEN
      WRITE(*,*)
      WRITE(*,*)' NUM*2.GT.MOM !!!!!'
      WRITE(*,*)' Increase MOM in the parameter statement!'
      STOP
      ENDIF
C
      DO 350 J = NDATA+1,NUM+1
           FX1(J) = CMPLX(0.0,0.0)
           FY1(J) = CMPLX(0.0,0.0)
           FZ1(J) = CMPLX(0.0,0.0)
  350 CONTINUE
C
C * Define some variables (First freq., duration of seismogram, number
C     of IFT data points, total number of output points, delta time
C     of seismogram)
C
      PI = ACOS(-1.0)
c      FOM=0.0
      DUR = (2.0*PI)/DOM
      NFAST = NUM*2
      DT = DUR/FLOAT(NFAST)
C
C * INVERSE FOURIER TRANSFORM  *
C
      DO 400 I=2,NFAST/2
           FX1(NFAST+2-I)=CMPLX(REAL(FX1(I)),-AIMAG(FX1(I)))
           FY1(NFAST+2-I)=CMPLX(REAL(FY1(I)),-AIMAG(FY1(I)))
           FZ1(NFAST+2-I)=CMPLX(REAL(FZ1(I)),-AIMAG(FZ1(I)))
  400 CONTINUE
C
      CALL WVFAST(NFAST,FX1,MOM,1)
      CALL WVFAST(NFAST,FY1,MOM,1)
      CALL WVFAST(NFAST,FZ1,MOM,1)
C
C * WAVE DATA *
C
      DO 450 I = 1,NFAST
         TIME(I) = FLOAT(I-1)*DT
         XR1(I)=REAL(FX1(I))/DUR
         YR1(I)=REAL(FY1(I))/DUR
         ZR1(I)=REAL(FZ1(I))/DUR
  450 CONTINUE
C
      WRITE(*,*)
      WRITE(*,1550) 
 1550 FORMAT(1X,'Change from (m) to (cm)? (yes=1)    : ')
      READ(*,*)ICM
      ACM=1.0
      IF(ICM.EQ.1)ACM=100.0
C
      DO 462 I=1,NFAST+1
        XR1(I)=ACM*XR1(I)
        YR1(I)=ACM*YR1(I)
        ZR1(I)=ACM*ZR1(I)
  462 CONTINUE
C
C * MAX AND MIN AMPLITUDE *
C
      XMAX=0.0 
      YMAX=0.0 
      ZMAX=0.0 
      XMIN=0.0 
      YMIN=0.0 
      ZMIN=0.0 
C
      DO 550 I=1,NFAST
        IF(XR1(I).GT.XMAX)XMAX=XR1(I)
        IF(XR1(I).LT.XMIN)XMIN=XR1(I)
        IF(YR1(I).GT.YMAX)YMAX=YR1(I)
        IF(YR1(I).LT.YMIN)YMIN=YR1(I)
        IF(ZR1(I).GT.ZMAX)ZMAX=ZR1(I)
        IF(ZR1(I).LT.ZMIN)ZMIN=ZR1(I)
  550 CONTINUE
C
      AMAX=XMAX
      AMIN=XMIN
      IF(YMAX.GT.AMAX)AMAX=YMAX
      IF(ZMAX.GT.AMAX)AMAX=ZMAX
      IF(YMIN.LT.AMIN)AMIN=YMIN
      IF(ZMIN.LT.AMIN)AMIN=ZMIN
C
      AAMAX=AMAX
      IF(ABS(AMIN).GT.AAMAX)AAMAX=ABS(AMIN)
C
      WRITE(*,*)
      WRITE(*,*)' Nmber of data points    =',NFAST
      WRITE(*,*)' Duration and Delta Time =',DUR,DT
      WRITE(*,*)' Initial and final times =',0.0,DUR-DT
      WRITE(*,*)
      WRITE(*,*)' Min X, Y and Z          =',XMIN,YMIN,ZMIN
      WRITE(*,*)' Max X, Y abd Z          =',XMAX,YMAX,ZMAX
      WRITE(*,*)' Min and Max amplitudes  =',AMIN,AMAX
      WRITE(*,*)' Absolute Max amplitude  =',AAMAX
C      
C * Output data file  *
C
      WRITE(*,*)
      WRITE(*,*)'Output data are stored in "wave1_dat"'
      WRITE(*,*)
      OPEN(11,FILE='wave1_dat',STATUS='unknown')
c      WRITE(11,*)NFAST,DT,AMIN,AMAX
      do i=1,NFAST
      write(11,2000)TIME(I),XR1(I),YR1(I),ZR1(I)
      end do
      CLOSE(11)
 2000 FORMAT(F8.2,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,NQF,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='soil_dat',STATUS='OLD')
      READ(4,1010)THET1
      READ(4,1020)NSF
      READ(4,1030)NL,NQF
      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(2I10)
 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,NQF
     *                 ,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 ** NQF    : = 0 : QF IN,   = 1 : QF NOT IN
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
      VVN2=DCMPLX(VV(NL)*VV(NL),0.0)
      IF(NQF.EQ.1)VVN2=VVN2*CMPLX(1.0,1.0/QF(NL))
      TSIN=SIN(THETA)
      CINV2=TSIN*TSIN/VVN2
C
      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
C
        R11=(1.0,0.0)
        R12=(0.0,0.0)
        R21=(0.0,0.0)
        R22=(1.0,0.0)
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(NQF.EQ.1)VV2=VV2*CMPLX(1.0,1.0/QF(IL))
          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
