C
C *** 2-D SV HASKELL'S FAULT MODEL ***
C
      SUBROUTINE HASKSV(DOM,NOM,MOM,XF,ZF,XO,ZO
     *                 ,VP,VS,W,DA,D,VRR,RT,UF,WF)
C
C ** MARAMETER **
C MOM : NOM IN MAIN PROGRAM DIMENSION
C     : IF MOM > 2**13(=8192), CAN NOT USE IN SUBROUTINE !!!
C ** INPUT DATA **
C DOM : INCREMENT OF CIRCULAR FREQUENCY ( DELTA OMEGA )
C NOM : TOTAL NUMBER OF CIRCULAR FREQUENCY
C XF  : FAULT LOCATION ON X-AXIS 
C ZF  : FAULT LOCATION ON Z-AXIS 
C XO  : LOCATION OF OBSERVATION POINT ON X-AXIS 
C ZO  : LOCATION OF OBSERVATION POINT ON Z-AXIS 
C VP  : COMPRESSIONAL WAVE VELOCITY
C VS  : SHEAR WAVE VELOCITY
C W   : WIDTH OF FAULT
C DA  : DIP ANGLE OF FAULT (DEGREE--->RADIAN IN THIS SUBROUTINE)
C D   : DISPLACEMENT DIFFELENCE OF FAULT
C VRR : RUPTURE VELOCITY OF FAULT
C RT  : RISE TIME
C ** CALCULATE DATA **
C UF(I) : X DISPLACEMENT OF Ith OMEGA AT (XO,ZO) BY FAULT
C WF(I) : Z DISPLACEMENT OF Ith OMEGA AT (XO,ZO) BY FAULT
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 H1S,H1P,GXXX,GXXZ,GZZX,GZZZ,GXZX,GXZZ,FX,FZ
     *          ,UF(MOM),WF(MOM),CCC1,H0P,H0S,H2P,H2S,CCC2
C
      PAI=3.141592654D0
      T=2.D0*PAI/DOM
      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)
C
      RFX=XF-XO
      RFZ=ZF-ZO
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,NOM
C
        OM=DOM*DFLOAT(K)
        AS=OM*RF/VS
        AP=OM*RF/VP
        Y1S=BESL21(AS)
        Y1P=BESL21(AP)
        B1S=BESL11(AS)
        B1P=BESL11(AP)
        H1S=DCMPLX(B1S,-Y1S)
        H1P=DCMPLX(B1P,-Y1P)
        Y0S=BESL20(AS)
        Y0P=BESL20(AP)
        B0S=BESL10(AS)
        B0P=BESL10(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*RT/2.D0
        XW=OM*RW/2.D0
        ST=DSIN(XT)/XT
        SW=DSIN(XW)/XW
        CCC1=CC*ST*SW*CDEXP(DCMPLX(0.D0,-XT-XW))
        UF(K)=CCC1*FX
        WF(K)=CCC1*FZ
C
  100 CONTINUE
C
      STOP
      END
