C 
C *** Wave Simulaion using Amplification Factor	for S and P wave
c     using the Haskell's method
C
      PROGRAM multi3x
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=18500)
      CHARACTER INFILE1*30,head1*10,head2*10,head3*10,head4*10
      COMPLEX Samp(MOM,2),Pamp(MOM,2)
      COMPLEX FX1(mom),FY1(mom),FZ1(mom),Sdelay,Pdelay
      DIMENSION DNS(ML),VP(ML),QP(ML),FP(ML)
     *         ,VS(ML),QS(ML),FS(ML),THK(ML)
     * ,TIME(mom),XR1(mom),XI1(mom),YR1(mom),YI1(mom),ZR1(mom),ZI1(mom)
C
      WRITE(*,*)
      WRITE(*,*)' : Soil Data must be "soil_dat.csv"'
C
C  Prompt for wave file
C
      WRITE(*,1010)
 1010 FORMAT(/1X,'Enter Wave Data (H1,H2,V1; e.g., wave.csv)')
      READ(*,1020) INFILE1
 1020 FORMAT(A20)
C
      OPEN(11,FILE=INFILE1,status='old')
        read(11,*)head1,head2,head3,head4
        i=1
   10   read(11,*,end=11)time(i),xr1(i),yr1(i),zr1(i)
c        write(*,*)time(i),xr1(i),yr1(i),zr1(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)
        FY1(I)=CMPLX(yr1(I),0.0)
        FZ1(I)=CMPLX(zr1(I),0.0)
      end do
      do I=NT+1,NF
        FX1(I)=(0.0,0.0)
        FY1(I)=(0.0,0.0)
        FZ1(I)=(0.0,0.0)
      end do
C
      CALL WVFAST(NF,FX1,MOM,-1)
      CALL WVFAST(NF,FY1,MOM,-1)
      CALL WVFAST(NF,FZ1,MOM,-1)
      DUR=FLOAT(NF)*DT
      do I=1,NF
        FX1(I)=FX1(I)*DT
        FY1(I)=FY1(I)*DT
        FZ1(I)=FZ1(I)*DT
      end do
c
      DOM=6.2831853/DUR
      IOM=1
      NOM=NF/2+1
      EOM=DOM*NOM
      WRITE(*,*)
      WRITE(*,*)' DUR,DOM,NOM,EOM=',DUR,DOM,NOM,EOM
      WRITE(*,*)' Highest Frequency (hz)=',EOM/6.2831853
C
C * Read data *
c
      CALL SOILIN(DNS,VP,QP,FP,VS,QS,FS,THK,THETA,NL,ML,NIN)
c
c *** Calculation of Amplification Factor ***
c
      FOM=0.0
      CALL SBSHML(FOM,DOM,NOM,MOM,DNS,VS,QS,FS,THK,NL,ML,THETA,Samp)
      CALL SBUDML(FOM,DOM,NOM,MOM,DNS,VP,QP,FP,VS,QS,FS,THK
     *           ,NL,ML,THETA,Pamp)
c
      open(9,file='Amp-Factors.csv',status='unknown')
      write(9,*)'frequency (Hz),Samp1,Samp2,Pamp1,Pamp2'
      do i=1,NOM
        om=dom*(i-1)
        freq=om/6.2831853
        write(9,1140)freq,cabs(Samp(i,1)),cabs(Samp(i,2))
     *   ,cabs(Pamp(i,1)),cabs(Pamp(i,2))
        Sdelay=(1.0,0.0)
        Pdelay=(1.0,0.0)
c        do l=1,nl-1
c          Sdelay=Sdelay*cexp(cmplx(0.0,om*thk(l)/vs(l)))
c          Pdelay=Pdelay*cexp(cmplx(0.0,om*thk(l)/vp(l)))
c        end do
        fx1(i)=fx1(i)*Samp(i,NIN)*Sdelay
        xr1(i)=real(fx1(i))
        xi1(i)=aimag(fx1(i))
        fy1(i)=fy1(i)*Samp(i,NIN)*Sdelay
        yr1(i)=real(fy1(i))
        yi1(i)=aimag(fy1(i))
        fz1(i)=fz1(i)*Pamp(i,NIN)*Pdelay
        zr1(i)=real(fz1(i))
        zi1(i)=aimag(fz1(i))
      end do
      close(9)
 1140 format(f8.4,6(',',e12.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
      pi2=atan(1.0)*8.0
      OML1=0.0
      OMH2=99999.9
      IF(NBD.EQ.1)THEN
        f1=DOM*(IOM-1)/pi2
        f2=DOM*(NOM-1)/pi2
        WRITE(*,*)'   Data for TRAPEZOIDAL BAND PASS FILTER:'
        WRITE(*,*)'   INITIAL AND FINAL Frequencies =',f1,f2
        WRITE(*,*)
        WRITE(*,*)
     *  '   Enter fL1,fL2,fH1 and fH2 (fL1<fL2<fH1<fH2: Hz): '
        READ(*,*)fL1,fL2,fH1,fH2
        OML1=fL1*pi2
        OML2=fL2*pi2
        OMH1=fH1*pi2
        OMH2=fH2*pi2
      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,NOM
           OMEGA=FLOAT(I-1)*DOM
           IF(OMEGA.lt.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*(0.0,1.0)
             FY1(I)=FY1(I)*OMEGA*(0.0,1.0)
             FZ1(I)=FZ1(I)*OMEGA*(0.0,1.0)
           END IF
c        write(*,1500)omega,FX1(I),FY1(I),FZ1(I)
  200 CONTINUE
c
c * Output Fourier Spectrum data *
c
c
      write(*,*)
      write(*,*)'The Fourier Spectrum data are stored in "amp.csv"'
      open(10,file='amp1.csv',status='unknown')
        write(10,*)'frequency(Hz),',head2,',',head3,',',head4
        DO 250 I = 2,NOM
          OM=FLOAT(I-1)*DOM
          f=om/pi2
          x=cabs(FX1(I))
          y=cabs(FY1(I))
          z=cabs(FZ1(I))
          write(10,1500)f,x,y,z
  250   CONTINUE
      close(10)
 1500 format(f8.4,6(',',e12.4))
C
C * INVERSE FOURIER TRANSFORM  *
C
      DO 400 I=2,NF/2
           FX1(NF+2-I)=CMPLX(REAL(FX1(I)),-aIMAG(FX1(I)))
           FY1(NF+2-I)=CMPLX(REAL(FY1(I)),-aIMAG(FY1(I)))
           FZ1(NF+2-I)=CMPLX(REAL(FZ1(I)),-aIMAG(FZ1(I)))
  400 CONTINUE
C
      CALL WVFAST(NF,FX1,MOM,1)
      CALL WVFAST(NF,FY1,MOM,1)
      CALL WVFAST(NF,FZ1,MOM,1)
C
C * WAVE DATA *
C
      DO 450 I = 1,NF
         XR1(I)=REAL(FX1(I))/DUR
         YR1(I)=REAL(FY1(I))/DUR
         ZR1(I)=REAL(FZ1(I))/DUR
  450 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,NF
        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    =',NF
      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.csv"'
      WRITE(*,*)'You can plot the waveform using "wave.m" for Matlab'
      WRITE(*,*)
      OPEN(11,FILE='wave1.csv',STATUS='unknown')
        write(11,*)head1,',',head2,',',head3,',',head4
        do i=1,NF
         write(11,2000)TIME(I),XR1(I)-xr1(1),YR1(I)-yr1(1),ZR1(I)-zr1(1)
        end do
      CLOSE(11)
 2000 FORMAT(F10.5,3(',',E12.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.1415926
      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,FP,VS,QS,FS,THK,THETA,NL,ML,NIN)
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),Fp(ml)
     *         ,Vs(ML),Qs(ML),Fs(ml),thk(ml)
C
C *** SOIL DATA INPUT ***
C 
      write(*,*)
      write(*,*)'From Subroutine SOILIN:'
      OPEN(4,FILE='soil_dat.csv',STATUS='OLD')
        read(4,*)
        READ(4,*)THET1
        write(*,*)'  Incident Wave Angle (degree) =',THET1
        READ(4,*)NIN
        write(*,*)'  Input Wave Type =', NIN
        write(*,*)'    (=1: Data of Theoretical Input Data in Bedrock'
        write(*,*)'    ,=2: Data Observed at Boreholein Bedrock)'
        read(4,*)
        READ(4,*)NL
        WRITE(*,2200)NL
        read(4,*)
        DO 100 L=1,NL
        READ(4,*)i,DNS(L),VP(L),QP(L),FP(L),VS(L),QS(L),FS(L),THK(L)
        WRITE(*,2210)i,DNS(L),VP(L),QP(L),FP(L),VS(L),QS(L),FS(L),THK(L)
  100   CONTINUE
      Close(10)
 2200 FORMAT(/'   Total Number of Layers =',I5/,'   Layer,Density,Vp(rea
     *l),    Qp  ,  Fp  ,Vs(real),    Qs  ,  Fs  ,Thickness(m)')
 2210 FORMAT(I6,F9.3,2(F9.2,F9.3,F7.3),F12.3)
C
C *** CHANGE THETA FROM DEGREE TO RADIAN ***
C
      THETA=THET1*0.017453292
c
      RETURN
      END
c
C
C *** CUL. OF AMPLITUDE FACTER BY MULTIPLE REFLECTION METHOD OF SH WAVE
C
      SUBROUTINE SBSHML(FOM,DOM,NOM,MOM,DNS,VS,QS,FS,THK
     *                 ,NL,ML,THETA,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 ** DNS(I) : DENSITY OF Ith LAYER
C ** VS (I) : VS OF Ith LAYER
C ** QS(I)  : QUALITY FACTER at the zero frequency of Ith LAYER (Qs0)
c ** FS(I)  : Qs(f)=QS0*f**FS(I), and Qs(f) > Qs0
C ** THK(I) : 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
C OUTPUT DATA
C ** AMP(I,1) : AMPLITUDE FACTER ON SURFACE at Ith OMEGA, withe respect to the input wave
C ** AMP(I,2) : AMPLITUDE FACTER ON SURFACE at Ith OMEGA, withe respect to the bedrock motion
C
      DIMENSION DNS(ML),THK(ML),VS(ML),QS(ML),FS(ML)
      COMPLEX*8 AMP(MOM,2),VSN2,CINV2,A11,A12,A21,A22,VS1,VS2,VS21
     *      ,am,ex1,ex2,WN,WND,SA11,SA12,SA21,SA22,SB11,SB12,SB21,SB22
C
C *** CAL. OF AMPLITUDE FACTER ***
C
      VSN2=CMPLX(VS(NL)*VS(NL),0.0)
      TSIN=SIN(THETA)
      CINV2=TSIN*TSIN/VSN2
C
c *** Free Surface Coefficient ***
      SF=2.0
C
c *** Amplicication Factor using Haskell's Matrix ***
      IFOM=1
      IF(FOM.EQ.0.0)THEN
      AMP(1,1)=CMPLX(SF,0.0)
      AMP(1,2)=CMPLX(1.0,0.0)
      IFOM=2
      END IF
C
      DO 200 IO=IFOM,NOM
C
        OM=FLOAT(IO-1)*DOM+FOM
        freq=om/6.2831853
C
        A11=(1.0,0.0)
        A12=(0.0,0.0)
        A21=(0.0,0.0)
        A22=(1.0,0.0)
C
        DO 100 IL=NL-1,1,-1
C
          SA11=A11
          SA12=A12
          SA21=A21
          SA22=A22
c
          VS1=VS(IL)
          VS2=VS(IL+1)
          if(QS(IL).eq.0.0)goto 120
c
          Q1=QS(IL)*freq**FS(IL)
          if(Q1.lt.QS(IL))Q1=QS(IL)
          VS1=VS1*CMPLX(1.0,0.5/Q1)
          Q2=QS(IL+1)*freq*FS(IL+1)
          if(Q2.lt.QS(IL+1))Q2=QS(IL+1)
          VS2=VS2*CMPLX(1.0,0.5/Q2)
c
  120     am=DNS(IL)*VS1/(DNS(IL+1)*VS2)
          VS21=VS1*VS1
          WN=OM*CSQRT(1.0/VS21-CINV2)
          WND=WN*THK(IL)*cmplx(0.0,1.0)
          ex1=cexp(+WND)
          ex2=cexp(-WND)
C
           SB11=(1.0+am)*ex1/2.0
           SB12=(1.0-am)*ex2/2.0
           SB21=(1.0-am)*ex1/2.0
           SB22=(1.0+am)*ex2/2.0
C
          A11=SA11*SB11+SA12*SB21
          A12=SA11*SB12+SA12*SB22
          A21=SA21*SB11+SA22*SB21
          A22=SA21*SB12+SA22*SB22
C
  100   CONTINUE
C
        AMP(IO,1)=SF/(A11+A12)
        AMP(IO,2)=2.0/(A11+A12+A21+A22)
c	write(*,*)'om,amp2',om,AMP(IO,2)
c	if(io.eq.20)stop
C
  200 CONTINUE
C
      RETURN
      END
c
C
C *** CUL. OF AMPLITUDE FACTER BY MULTIPLE REFLECTION METHOD OF UD WAVE
C
      SUBROUTINE SBUDML(FOM,DOM,NOM,MOM,DNS,VP,QP,FP,VS,QS,FS,THK
     *                 ,NL,ML,THETA,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 ** DNS(I) : DENSITY OF Ith LAYER
C ** VS (I) : VS OF Ith LAYER
C ** THK(I) : THICKNESS OF Ith LAYER
c ** Qp(I)  : QUALITY FACTER of P Wave at the zero frequency of Ith LAYER (Qp0)
c ** Fp(I)  : Qp(f)=Qp0*f**Fp(I), and Qp(f) > Qp0
c ** Qs(I)  : QUALITY FACTER of S Wave at the zero frequency of Ith LAYER (Qs0)
c ** Fs(I)  : Qs(f)=QS0*f**FS(I), and Qs(f) > Qs0
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,1) : AMPLITUDE FACTER ON SURFACE at Ith OMEGA, withe respect to the input wave
C ** AMP(I,2) : AMPLITUDE FACTER ON SURFACE at Ith OMEGA, withe respect to the bedrock motion
C
      DIMENSION DNS(ML),THK(ML),VP(ML),QP(ML),FP(ML)
     *         ,VS(ML),QS(ML),FS(ML)
      COMPLEX*8 AMP(MOM,2),VSN2,CINV2,A11,A12,A21,A22,VS1,VS2,VP1,VP2
     *       ,am,ex1,ex2,WN,WND,SA11,SA12,SA21,SA22,SB11,SB12,SB21,SB22
     *       ,aa,vsp1,vsp2,vv2
C
C *** CAL. OF AMPLITUDE FACTER ***
C
      aa=3.0/(1.0-VS(NL)*VS(NL)/(VP(NL)*VP(NL)))
      VSN2=CMPLX(VS(NL)*VS(NL),0.0)*aa
c	write(1,*)'NL,VSP=',NL,csqrt(VSN2)
      TSIN=SIN(THETA)
      CINV2=TSIN*TSIN/VSN2
C
c *** Free Surface Coefficient ***
      SF=2.0
C
c *** Amplicication Factor using Haskell's Matrix ***
      IFOM=1
      IF(FOM.EQ.0.0)THEN
      AMP(1,1)=CMPLX(SF,0.0)
      AMP(1,2)=CMPLX(1.0,0.0)
      IFOM=2
      END IF
C
      DO 200 IO=IFOM,NOM
C
        OM=FLOAT(IO-1)*DOM+FOM
        freq=om/6.2831853
C
        A11=(1.0,0.0)
        A12=(0.0,0.0)
        A21=(0.0,0.0)
        A22=(1.0,0.0)
C
        DO 100 IL=NL-1,1,-1
C
          SA11=A11
          SA12=A12
          SA21=A21
          SA22=A22
c
          VS1=VS(IL)
          VS2=VS(IL+1)
          VP1=VP(IL)
          VP2=VP(IL+1)
          if(QS(IL).eq.0.0)goto 120

          Q1=QS(IL)*freq**FS(IL)
          if(Q1.lt.QS(IL))Q1=QS(IL)
          VS1=VS1*CMPLX(1.0,0.5/Q1)
          Q2=QS(IL+1)*freq*FS(IL+1)
          if(Q2.lt.QS(IL+1))Q2=QS(IL+1)
          VS2=VS2*CMPLX(1.0,0.5/Q2)
c
          Q1=QP(IL)*freq**FP(IL)
          if(Q1.lt.QP(IL))Q1=QP(IL)
          VP1=VP1*CMPLX(1.0,0.5/Q1)
          Q2=QP(IL+1)*freq*FP(IL+1)
          if(Q2.lt.QP(IL+1))Q2=QP(IL+1)
          VP2=VP2*CMPLX(1.0,0.5/Q2)
C
  120     aa=csqrt(3.0/(1.0-VS1*VS1/(VP1*VP1)))
          vsp1=VS1*aa
          aa=csqrt(3.0/(1.0-VS2*VS2/(VP2*VP2)))
          vsp2=VS2*aa
c
          am=dns(IL)*vsp1/(dns(IL+1)*vsp2)
          vv2=vsp1*vsp1
          WN=OM*CSQRT(1.0/vv2-CINV2)
          WND=WN*THK(IL)*cmplx(0.0,1.0)
          ex1=cexp(+WND)
          ex2=cexp(-WND)
C
           SB11=(1.0+am)*ex1/2.0
           SB12=(1.0-am)*ex2/2.0
           SB21=(1.0-am)*ex1/2.0
           SB22=(1.0+am)*ex2/2.0
C
          A11=SA11*SB11+SA12*SB21
          A12=SA11*SB12+SA12*SB22
          A21=SA21*SB11+SA22*SB21
          A22=SA21*SB12+SA22*SB22
C
  100   CONTINUE
C
        AMP(IO,1)=SF/(A11+A12)
        AMP(IO,2)=2.0/(A11+A12+A21+A22)
C
  200 CONTINUE
C
      RETURN
      END
