C
C *** 2-D SH HASKELL'S FAULT MODEL WITH SOLUSION IN HALF SPACE ***
C
      SUBROUTINE HKSHS1(OM,XF,ZF,XO,ZO,VS,W,DA,D,VRR,RT,VF)
C
C ** INPUT DATA **
C OM  : OMEGA
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 VS  : SHEAR WAVE VELOCITY
C W   : WIDTH OF FAULT
C DA  : DIP ANGLE OF FAULT (RADIAN)
C D   : DISPLACEMENT DIFFELENCE OF FAULT
C VRR : RUPTURE VELOCITY OF FAULT
C RT  : RISE TIME
C ** CALCULATE DATA **
C VF  : Y DISPLACEMENT AT (XO,ZO) BY FAULT
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 H1S1,H1S2,VF
C
      RW=W/VRR
      CC=D*W/(4.D0*VS)
C
      RFX=XF-XO
      RFZ1=ZF-ZO
      RFZ2=ZF+ZO
      XN= DSIN(DA)
      ZN=-DCOS(DA)
      RF1=DSQRT(RFX*RFX+RFZ1*RFZ1)
      RF2=DSQRT(RFX*RFX+RFZ2*RFZ2)
      RXZ1=(RFX*XN+RFZ1*ZN)/RF1
      RXZ2=(RFX*XN+RFZ2*ZN)/RF2
C
      AS1=OM*RF1/VS
      AS2=OM*RF2/VS
      CALL BESJN1(AS1,B1S1,Y1S1)
      CALL BESJN1(AS2,B1S2,Y1S2)
      H1S1=DCMPLX(B1S1,-Y1S1)
      H1S2=DCMPLX(B1S2,-Y1S2)
      XT=OM*RT/2.D0
      XW=OM*RW/2.D0
      ST=DSIN(XT)/XT
      SW=DSIN(XW)/XW
      VF=CC*(H1S1*RXZ1+H1S2*RXZ2)*ST*SW*CDEXP(DCMPLX(0.D0,-XT-XW))
C
      RETURN
      END
