C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 H1S,H1P,GXXX,GXXZ,GZZX,GZZZ,GXZX,GXZZ,FX,FZ
     *          ,UF(2000),WF(2000),CCC1,H0P,H0S,H2P,H2S,CCC2
     *          ,UFS(2000),WFS(2000),EX
      DIMENSION OM(2000),UD(2000),UV(2000),UA(2000)
     *         ,DT(2000),WD(2000),WV(2000),WA(2000)
     *         ,RC(0:1000),RX(100),DKL(1000)
C
      WRITE(6,1000)
      WRITE(6,1200)
C     READ(5,1100)VS
      WRITE(6,1250)
C     READ(5,1100)VP
      WRITE(6,1300)
C     READ(5,1100)RX0
      WRITE(6,1350)
C     READ(5,1100)RZ0
      WRITE(6,1280)
      READ(5,1150)NR
      DO 10 I=1,NR
      WRITE(6,1360)I
      READ(5,1100)RX(I)
   10 CONTINUE
      WRITE(6,1400)
C     READ(5,1100)W
      WRITE(6,1450)
C     READ(5,1100)DA
      WRITE(6,1500)
C     READ(5,1100)D
      WRITE(6,1550)
C     READ(5,1100)VRR
      WRITE(6,1650)
C     READ(5,1100)RT
      WRITE(6,1700)
      READ(5,1100)DDO
      WRITE(6,1750)
      READ(5,1150)KK
      WRITE(6,1800)
C     READ(5,1100)DDT
      WRITE(6,1850)
C     READ(5,1150)NN
      WRITE(6,1900)
      READ(5,1160)NF
 1000 FORMAT(/' 2-D.Rayleigh Waves in 2-L from 2-D. Fault MODEL'
     *      //' Basin edge is assumed at ( 0, 0 )'
     *      //' NEED PAPER 132 LINESIZE by TERM LINESIZE(132) !!!'
     *      //' NEED MSL2 !!!'//)
 1100 FORMAT(F10.0)
 1150 FORMAT(I3)
 1160 FORMAT(I1)
 1200 FORMAT(/' VS      (  M / SEC ) ? = ')
 1250 FORMAT( ' VP      (  M / SEC ) ? = ')
 1280 FORMAT( ' NUM .OF OUTPUT O.P. ( I3 !!!!) ? = ')
 1300 FORMAT( ' ; X-COMP. OF FAULT LOCATION ( M ) : RFX ? = ')
 1350 FORMAT( ' ; Z-COMP. OF FAULT LOCATION ( M ) : RFZ ? = ')
 1360 FORMAT(I3,' ; X-COMP. OF O.P. LOCATION ( M ) : RX ? = ')
 1400 FORMAT( ' WIDTH OF FAULT ( M )   : W    ? = ')
 1450 FORMAT( ' DIP ANGLE ( DEGREE )   : DA   ? = ')
 1500 FORMAT( ' DISP. OF FAULT ( M )   : D    ? = ')
 1550 FORMAT( ' RUPTURE VEL. ( M/SEC ) : VRR  ? = ')
 1650 FORMAT( ' RISE TIME ( SEC )      : RT   ? = ')
 1700 FORMAT( ' STEP OF OMEGA ( RAD/SEC )      : DDO ? = ')
 1750 FORMAT( ' NUM. OF OMEGA STEP ( I3 !!!! ) : KK  ? = ')
 1800 FORMAT( ' STEP OF TIME ( SEC )           : DDT ? = ')
 1850 FORMAT( ' NUM. OF TIME STEP ( I3 !!!! )  : NN  ? = ')
 1900 FORMAT( ' FILE OUT...1      ( I1 !!!! )  : NF  ? = ')
      VS=3000.D0
      VP=5500.D0
      D=1.D0
      W=3141.5927D0
      RT=3.1415927D0
      VRR=1000.D0
      DA=80.D0
      RX0=-4696.155D0
      RZ0= 3472.9635D0
C     NR=1
      DDT=0.D0
      NN=0
      WRITE(6,2000)VS,VP,D,W,RT
      WRITE(6,2050)VRR,DA,NR
      WRITE(6,2100)DDO,KK,DDT,NN
      WRITE(6,2150)RX0,RZ0
      WRITE(6,2160)(RX(I),I=1,NR)
 2000 FORMAT(//'  VS =',F10.4,'  VP =',F10.4,'  D =',F8.4,'   W =',F10.4
     *        ,' RT =',F8.4)
 2050 FORMAT(' VRR =',F10.3,' DA  =',F8.3,' NR  =',I4)
 2100 FORMAT(' DDO = ',F8.4,'  KK = ',I4,'  DDT = ',F8.4,'  NN = ',I4)
 2150 FORMAT(' RX0,RZ0  = ',7F11.3)
 2160 FORMAT(' RX  = ',7F11.3)
C
        WRITE(7,7005)DDO,KK
        WRITE(7,7000)RFX,RFZ,RX(IR)
        WRITE(7,*)
        WRITE(7,*)
        WRITE(7,*)
 7000 FORMAT(3F15.4)
 7005 FORMAT(F15.8,I15)
C
      PAI=3.141592654D0
      T=2.D0*PAI/DDO
      DA=DA*PAI/180.0D0
      SIN2=DSIN(2.D0*DA)
      COS2=DCOS(2.D0*DA)
      VSP3=VS*VS*VS/(VP*VP*VP)
      RW=W/VRR
      CC=D*W/(4.D0*VS)
      WRITE(6,2500)T,RW
 2500 FORMAT(/'  T (SEC)=',F8.4,'   RW (=W/VR; SEC)=',F9.4)
C
      RFX=RX0
      RFZ=RZ0
C
      RF=DSQRT(RFX*RFX+RFZ*RFZ)
      RXRF=RFX/RF
      RZRF=RFZ/RF
      RFX2=RXRF*RXRF
      RFZ2=RZRF*RZRF
      RFXZ=RXRF*RZRF
C
      DO 100 K=1,KK
      OM(K)=DDO*DFLOAT(K)
      AS=OM(K)*RF/VS
      AP=OM(K)*RF/VP
      Y1S=$DMY1M(AS)
      Y1P=$DMY1M(AP)
      B1S=$DMJ1M(AS)
      B1P=$DMJ1M(AP)
      H1S=DCMPLX(B1S,-Y1S)
      H1P=DCMPLX(B1P,-Y1P)
      Y0S=$DMY0M(AS)
      Y0P=$DMY0M(AP)
      B0S=$DMJ0M(AS)
      B0P=$DMJ0M(AP)
      H0S=DCMPLX(B0S,-Y0S)
      H0P=DCMPLX(B0P,-Y0P)
      H2S=H1S*2.D0/AS-H0S
      H2P=H1P*2.D0/AP-H0P
      CCC2=VSP3*H2P/AP-H2S/AS
      R34X=RXRF*(3.D0-4.D0*RFX2)
      R34Z=RZRF*(3.D0-4.D0*RFZ2)
      R14X=RXRF*(1.D0-4.D0*RFZ2)
      R14Z=RZRF*(1.D0-4.D0*RFX2)
      CCC1=VSP3*RFX2*H1P+(1.D0-RFX2)*H1S
      GXXX=RXRF*CCC1+R34X*CCC2
      GXXZ=RZRF*CCC1+R14Z*CCC2
      CCC1=VSP3*RFZ2*H1P+(1.D0-RFZ2)*H1S
      GZZX=RXRF*CCC1+R14X*CCC2
      GZZZ=RZRF*CCC1+R34Z*CCC2
      CCC1=RFXZ*(VSP3*H1P-H1S)
      GXZX=RXRF*CCC1+R14Z*CCC2
      GXZZ=RZRF*CCC1+R14X*CCC2
      FX=SIN2*(GXXX-GXZZ)-COS2*(GXXZ+GXZX)
      FZ=SIN2*(GXZX-GZZZ)-COS2*(GXZZ+GZZX)
      XT=OM(K)*RT/2.D0
      XW=OM(K)*RW/2.D0
      ST=DSIN(XT)/XT
      SW=DSIN(XW)/XW
      CCC1=CC*ST*SW*CDEXP(DCMPLX(0.D0,-XT-XW))
      UFS(K)=CCC1*FX
      WFS(K)=CCC1*FZ
  100 CONTINUE
C
      CALL RAYC(VS,VP,DDO,KK,RC,DKL)
C
      DO 600 IR=1,NR
C
      DO 150 K=1,KK
      EX=CDEXP(DCMPLX(0.D0,-OM(K)*RX(IR)/RC(K)))
      WF(K)=2.D0*WFS(K)*EX
      UF(K)=WF(K)*DKL(K)*DCMPLX(0.D0,1.D0)
  150 CONTINUE
C
      WRITE(6,3000)IR,RX(IR)
C     DO 200 K=1,KK
C     WRITE(6,3100)OM(K),2.D0*PAI/OM(K),UF(K),WF(K),CDABS(UF(K))
C    *         ,CDABS(WF(K)),CDABS(UF(K))*OM(K),CDABS(WF(K))*OM(K)
C 200 CONTINUE
C
      IF(NF.EQ.1) THEN
        DO 300 K=1,KK,5
          WRITE(7,7010)(UF(K+J),J=0,4)
  300   CONTINUE
        DO 310 K=1,KK,5
          WRITE(7,7010)(WF(K+J),J=0,4)
  310   CONTINUE
      END IF
C
 3000 FORMAT(//'  IR =',I4,'   RX =',F11.3/
     *        ,'  OMEGA , PERIOD,  UF(REAL) ,  UF(IMAG) ,  WF(REAL) ,'
     *,'  WF(IMAG) ,  UF(DABS) ,  WF(DABS) ,UF(DABS)*OM,WF(DABS)*OM'/)
 3100 FORMAT(2F8.3,8D12.4)
 7010 FORMAT(10F8.5)
C
      STOP
      END
