C
C *** CUL. OF AMPLITUDE SPECTRUM 
C      BY MULTIPLE REFLECTION METHOD OF SH WAVE
C 
      SUBROUTINE WVSHML(FOM,DOM,NOM,MOM,DNS,VS,DEP,QF,NQF
     *                 ,NL,ML,THETA,NSF,FOUR)
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 ** VS (I) : VS OF Ith LAYER
C ** DEP(I) : DEPTH ( 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 ** FOUR(I): FOURIER SPECTRUM
C OUTPUT
C ** FOUR(I): AMPLITUDE FACTER * FOUR(I)
C
      DIMENSION DNS(ML),DEP(ML),VS(ML),QF(ML)
      COMPLEX*8 FOUR(MOM),VSN2,CINV2,R11,R12,R21,R22,VS2,SM,WN,WND,WNS
     *   ,COWN,SIWN,SA11,SA12,SA21,SA22,SB11,SB12,SB21,SB22,WNN,WNSN,AMP
C
C *** CAL. OF AMPLITUDE FACTER ***
C
      VSN2=CMPLX(VS(NL)*VS(NL),0.0)
      IF(NQF.EQ.1)VSN2=VSN2*CMPLX(1.0,1.0/QF(NL))
      TSIN=SIN(THETA)
      CINV2=TSIN*TSIN/VSN2
C
      SF=1.0
      IF(NSF.EQ.1)SF=2.0
C
      IFOM=1
      IF(FOM.EQ.0.0)THEN
      FOUR(1)=CMPLX(SF,0.0)*FOUR(1)
      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
          VS2=CMPLX(VS(IL)*VS(IL),0.0)
          IF(NQF.EQ.1)VS2=VS2*CMPLX(1.0,1.0/QF(IL))
          SM=DNS(IL)*VS2   
          WN=OM*CSQRT(1.0/VS2-CINV2)
          WND=WN*DEP(IL)
          WNS=WN*SM
          COWN=CCOS(WND)
          SIWN=CSIN(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*CSQRT(1.0/VSN2-CINV2)
        WNSN=WNN*DNS(NL)*VSN2
C
        AMP=SF/(R11-R21*CMPLX(0.0,1.0)/WNSN)
        FOUR(IO)=AMP*FOUR(IO)
C
  200 CONTINUE
C
      RETURN
      END
