C
C *** 2-D SH HASKELL'S FAULT MODEL WITH SOLUSION IN FULL SPACE ***
C
      SUBROUTINE HKSHF1(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 H1S,VF
C
      RW=W/VRR
      CC=D*W/(4.D0*VS)
C
      RFX=XF-XO
      RFZ=ZF-ZO
      XN= DSIN(DA)
      ZN=-DCOS(DA)
      RF=DSQRT(RFX*RFX+RFZ*RFZ)
      RXZ=(RFX*XN+RFZ*ZN)/RF
C
      AS=OM*RF/VS
      CALL BESJN1(AS,B1S,Y1S)
      H1S=DCMPLX(B1S,-Y1S)
      XT=OM*RT/2.D0
      XW=OM*RW/2.D0
      ST=DSIN(XT)/XT
      SW=DSIN(XW)/XW
      VF=CC*H1S*RXZ*ST*SW*CDEXP(DCMPLX(0.D0,-XT-XW))
C
      RETURN
      END
