C 
C *** AMPLITUDE FACTER OF MULTIPLE REFLECTION METHOD OF SH WAVE
C
      PROGRAM MSHML3
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 IN SUBROUTINE !!!
C ** DNS(I) : DENSITY OF Ith LAYER
C ** VS (I) : VS OF Ith LAYER
C ** DEP(I) : DEPTH ( THICKNESS ) OF Ith LAYER
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) : AMPLITUDE FACTER ON SURFACE OF Ith OMEGA ( (I-1)*DOM+FOM )
C ** SPCA(I): SPECTRUM AMPLITUDE OF AMP(I)
C
      PARAMETER(MOM=2000)
      complex amp(mom),fx1(mom),fy1(mom),fz1(mom),x(mom),y(mom),z(mom)
      REAL XR1(mom),XI1(mom),YR1(mom),YI1(mom),ZR1(mom),ZI1(mom)
      character name*20,fname*30
C
      WRITE(*,*)' Enter Soil Data File Name '
      READ(*,'(A)')FNAME

      OPEN (4,FILE=FNAME,status='old')
      READ(4,1000)NNN      
 1000 format(i10)
      write(*,*)'Number of Soil Data Files =',nnn
      READ(4,1100)FOM,DOM,NOM
      READ(4,1110)THET1
      READ(4,1120)NSF
      if(NSF.ne.0)then
        write(*,*)'nsf.eq.0 !!'
        stop
      end if
 1100 FORMAT(2F10.0,I10)
 1110 FORMAT(F10.0)
 1120 FORMAT(I10)
      close(4)
C
      write(*,*)
      write(*,*)'Acc (=1) or Vel(=2)?'
      read(*,*)nacc
      if(nacc.eq.1)then
        write(*,*)'Read Input Amplitude Data (from acc_amp)'
        OPEN (4,FILE='acc_amp',status='old')
      else
        write(*,*)'Read Input Amplitude Data (from vel_amp)'
        OPEN (4,FILE='vel_amp',status='old')
      end if

      i=1
   10 read(4,*,end=20)om,XR1(I),XI1(I),YR1(I),YI1(I),ZR1(I),ZI1(I)
      i=i+1
      goto 10
   20 n=i-1
      if(n.ne.nom)then
        write(*,*)'n.ne.nom !'
        stop
      end if
      close(4)
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,NOM
             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,NOM
             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,NOM
             ZR1(I)=-ZR1(I)
             ZI1(I)=-ZI1(I)
  184      CONTINUE
         END IF
C
      END IF
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*(NOM-1)
        WRITE(*,*)
        WRITE(*,*)
     *  '   Enter OML1,OML2,OMH1 and OMH2 (OML1<OML2<OMH1<OMH2): '
        READ(*,*)OML1,OML2,OMH1,OMH2
      END IF
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
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,NOM
           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))
  200 CONTINUE
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(NOM.LE.N)THEN
                NUM = N
                GO TO 30
           END IF
  300 CONTINUE
C
   30 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 30
      END IF

      IF(NUM*2.GT.MOM)THEN
        WRITE(*,*)
        WRITE(*,*)' NUM*2.GT.MOM !!!!!'
        STOP
      ENDIF
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)
      write(*,*)
      write(*,*)'NFAST,DUR,DT=',NFAST,DUR,DT
c
c *** Calculation of Amplified wave **+
c
      OPEN(7,FILE='spec_dat',status='old')
      OPEN(8,FILE='wave_dat',status='unknown')
      OPEN(9,FILE='wave1_dat',status='unknown')
      OPEN(10,FILE='wave_max',status='unknown')

      write(*,*)
      do 500 iii = 1, nnn

        read(7,'(a)')name
        DO 510 IOM=1,NOM
        read(7,2100)om,amp(IOM)
        x(iom)=amp(iom)*fx1(iom)
        y(iom)=amp(iom)*fy1(iom)
        z(iom)=amp(iom)*fz1(iom)
  510   CONTINUE
 2100   FORMAT(3f9.4)
C
        DO 520 J = NOM+1,NUM+1
         x(J) = CMPLX(0.0,0.0)
         y(J) = CMPLX(0.0,0.0)
         z(J) = CMPLX(0.0,0.0)
  520   CONTINUE
C
      DO 530 I=2,NFAST/2
           x(NFAST+2-I)=CMPLX(REAL(x(I)),-AIMAG(x(I)))
           y(NFAST+2-I)=CMPLX(REAL(y(I)),-AIMAG(y(I)))
           z(NFAST+2-I)=CMPLX(REAL(z(I)),-AIMAG(z(I)))
  530 CONTINUE
C
      CALL WVFAST(NFAST,x,MWV,1)
      CALL WVFAST(NFAST,y,MWV,1)
      CALL WVFAST(NFAST,z,MWV,1)
C
C * WAVE DATA *
C
      DO 540 I = 1,NFAST
         XR1(I)=acm*REAL(x(I))/DUR
         YR1(I)=acm*REAL(y(I))/DUR
         ZR1(I)=acm*REAL(z(I))/DUR
  540 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(8,'(a)')name
      do i=1,nfast
      write(8,3000)dt*(i-1),xr1(i),yr1(i),zr1(i)
      write(9,3000)dt*(i-1),xr1(i),yr1(i),zr1(i)
      end do
 3000 format(f8.4,3f8.2) 

      write(10,2100)aamax    
      write(*,*)'Max Amp =',aamax    
 3100 format(f8.2)

  500 continue

      CLOSE(7)
      CLOSE(8)
      CLOSE(9)
      CLOSE(10)
C
      STOP
      END
c
c
C
C *** FOURIER TRANSFORM OR INVERSE FOURIER TRANSFORM
C
      SUBROUTINE WVFAST(N,A,MOM,IND)
C
      COMPLEX A(MOM),TEMP,THETA
C
C * FFT and inverse 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=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

