C
C *** EIGEN VALUE EQUATION OF RAYLEIGH WAVE IN CASE OF CC=VS(NL)
C
      REAL FUNCTION FCFRN*8 (OM,ML,NL,DNS,VP,VS,DEP)
C
C INPUT DATA
C OM     : CIRCULAR FREQUENCY ( OMEGA ) FROM COMMON
C OUTPUT DATA
C FNCFRN : EIGEN VALUE EQUATION OF RAYLEIGH WAVE
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION VP(ML),VS(ML),DEP(ML),DNS(ML)
      COMPLEX*16 AA(4,4),SA(4,4),SB(4,4)
     *          ,RAN,B11,B13,RL,RK,RN,RM
C
C *** WAVE NUMBER ***
C
      BN=VS(NL)
      WN=OM/BN
C
C *** INITIALIZATION OF AA MATRIX ***
C
      DO 100 I=1,4
        DO 110 J=1,4
          AA(I,J)=(0.0D0,0.0D0)
  110   CONTINUE
        AA(I,I)=(1.0D0,0.0D0)
  100 CONTINUE
C
C *** AA MATRIX ***
C
      DO 200 IL=NL-1,1,-1
C
        DO 210 I=1,4
          DO 220 J=1,4
            SA(I,J)=AA(I,J)
  220     CONTINUE
  210   CONTINUE
C
C    ** CQ, STQ AND SDQ **
C
        BM=VS(IL)
        CB2=BN*BN/(BM*BM)
C
        RBM=DSQRT(CB2-1.0D0)
        QM=WN*RBM*DEP(IL)
        CQ=DCOS(QM)
        SQ=DSIN(QM)
        STQ=RBM*SQ
        SDQ=SQ/RBM
C
C    ** CP, STP AND SDP **
C
        AM=VP(IL)
        CA2=BN*BN/(AM*AM)
        IF(AM-BN)10,20,30
C
   10     RAM=DSQRT(CA2-1.0D0)
          PM=WN*RAM*DEP(IL)
          CP=DCOS(PM)
          SP=DSIN(PM)
          STP=RAM*SP
          SDP=SP/RAM
          GOTO 40
C
   20     CP=1.0D0
          STP=0.0D0
          SDP=WN*DEP(IL)
          GOTO 40
C
   30     RAM=DSQRT(1.0D0-CA2)
          PM=WN*RAM*DEP(IL)
          EX1=DEXP(PM)
          EX2=DEXP(-PM)
          CP=(EX1+EX2)/2.0D0
          EX12=(EX1-EX2)/2.0D0
          STP=-RAM*EX12
          SDP=EX12/RAM
C
C    ** C11-C16, D11-D15 AND SB MATRIX **
C
   40   GM=2.D0/CB2
C
        C11=GM
        C12=GM-1.D0
        C13=GM*GM
        C14=C12*C12
        C16=DNS(IL)*BN*BN
        C15=1.D0/C16
C
        D11=CP-CQ
        D12=SDP+STQ
        D13=STP+SDQ
        D14=C14*SDP+C13*STQ
        D15=C13*STP+C14*SDQ
C
        SB(1,1)= DCMPLX(C11*CP-C12*CQ,0.0D0)
        SB(1,2)= DCMPLX(0.0D0,C12*SDP+C11*STQ)
        SB(1,3)=-DCMPLX(C15*D11,0.D0)
        SB(1,4)= DCMPLX(0.D0,C15*D12)
        SB(2,1)=-DCMPLX(0.0D0,C11*STP+C12*SDQ)
        SB(2,2)= DCMPLX(-C12*CP+C11*CQ,0.D0)
        SB(2,3)= DCMPLX(0.D0,C15*D13)
        SB(2,4)= SB(1,3)
        SB(3,1)= DCMPLX(C16*C11*C12*D11,0.0D0)
        SB(3,2)= DCMPLX(0.D0,C16*D14)
        SB(3,3)= SB(2,2)
        SB(3,4)= SB(1,2)
        SB(4,1)= DCMPLX(0.D0,C16*D15)
        SB(4,2)= SB(3,1)
        SB(4,3)= SB(2,1)
        SB(4,4)= SB(1,1)
C
C    ** AA MATRIX **
C
        DO 250 I=1,4
          DO 260 J=1,4
            AA(I,J)=(0.D0,0.D0)
            DO 270 K=1,4
              AA(I,J)=AA(I,J)+SA(I,K)*SB(K,J)
  270       CONTINUE
  260     CONTINUE
  250   CONTINUE
C
  200 CONTINUE
C
C *** GN,RAN,B11-B24 AND RL-RM ***
C
      GN=2.D0*VS(NL)*VS(NL)/(BN*BN)
      RAN=-DCMPLX(0.D0,DSQRT(1.0D0-BN*BN/(VP(NL)*VP(NL))))
C
      B11=GN*RAN
      B12=GN-1.D0
      B14=1.D0/(DNS(NL)*BN*BN)
      B13=-RAN*B14
      B21=-B12
      B23=B14
C
      RL=B11*AA(1,1)+B12*AA(2,1)+B13*AA(3,1)+B14*AA(4,1)
      RK=B11*AA(1,2)+B12*AA(2,2)+B13*AA(3,2)+B14*AA(4,2)
      RN=B21*AA(1,1)+B23*AA(3,1)
      RM=B21*AA(1,2)+B23*AA(3,2)
C
C *** FNCFRN ***
C
      FCFRN=REAL(RN*RK-RL*RM)
C                     
      RETURN
      END
