C
C *** 2-D SH HASKELL'S FAULT MODEL WITH SOLUSION IN FULL SPACE ***
C
      SUBROUTINE HSKSHF(DOM,NOM,MOM,XF,ZF,XO,ZO,VS,W,DA,D,VRR,RT,VF)
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 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(I) : Y DISPLACEMENT OF Ith OMEGA AT (XO,ZO) BY FAULT
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 H1S,VF(MOM)
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
      DO 100 K=1,NOM
C
        OM=DOM*DFLOAT(K)
        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(K)=CC*H1S*RXZ*ST*SW*CDEXP(DCMPLX(0.D0,-XT-XW))
C
  100 CONTINUE
C
      RETURN
      END
