C
C *** INPUT OF FAULT RARAMETERS FROM FILE ***
C
      SUBROUTINE HSKIN(XF,ZF,DNS,VP,VS,W,DA,D,VRR,RT)
C
C** DATA **
C XF    : FAULT LOCATION ON X-AXIS 
C ZF    : FAULT LOCATION ON Z-AXIS 
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
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER FNAME*30
C
C *** INPUT OF FAULT PARAMETERS FROM FILE DATA ***
C
      WRITE(*,*)
      WRITE(*,1000)
   10 WRITE(*,1010)
      READ(*,'(A)')FNAME
      IF(FNAME(1:1).EQ.' ')GOTO 10
 1000 FORMAT( ' From SUBROUTINE H2DIN :')
 1010 FORMAT(/' Enter File Name of 2-D SH Haskell''s Fault Model :',\)

      OPEN (7,FILE=FNAME)
        READ(7,1200)XF,ZF
        READ(7,1200)W,DA,D
        READ(7,1200)VRR,RT
        READ(7,1200)DNS,VP,VS
      CLOSE(7)
 1200 FORMAT(3F10.0)
C
C *** DEGREE => RADIAN ***
C
      DA1=DA
      DA=DA*3.141592654D0/180.D0
C
C *** DATA ECHO ***
C
      WRITE(*,2000)XF,ZF
      WRITE(*,2010)W,D
      WRITE(*,2020)DA1,DA
      WRITE(*,2030)VRR,RT
      WRITE(*,2040)DNS,VP,VS
 2000 FORMAT(/' Fault Location (XF,ZF)         =',2F10.3)
 2010 FORMAT( ' Width and Dislocation          =',F10.3,F10.5)
 2020 FORMAT( ' Dip Angle (deg) and (rad)      =',2F10.5)
 2030 FORMAT( ' Rise Time and Rupture Velocity =',F10.3,F10.5)
 2040 FORMAT( ' Density, Vp and Vs             =',F10.3,2F10.5)
C
      RETURN
      END
C
C
C
C *** 2-D SV HASKELL'S FAULT MODEL ***
C
      SUBROUTINE HSKSV(OM,XF,ZF,XO,ZO,VP,VS,W,DA,D,VRR,RT,UF,WF)
C
C ** INPUT DATA **
C OM  : 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 (RADIAN)
C D   : DISPLACEMENT DIFFELENCE OF FAULT
C VRR : RUPTURE VELOCITY OF FAULT
C RT  : RISE TIME
C ** CALCULATE DATA **
C UF  : X DISPLACEMENT OF Ith OMEGA AT (XO,ZO) BY FAULT
C WF  : 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,WF,CCC1,H0P,H0S,H2P,H2S,CCC2
C
      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
      AS=OM*RF/VS
      AP=OM*RF/VP
C
      CALL BESJN0(AS,B0S,Y0S)
      CALL BESJN0(AP,B0P,Y0P)
      CALL BESJN1(AS,B1S,Y1S)
      CALL BESJN1(AP,B1P,Y1P)
C
      H0S=DCMPLX(B0S,-Y0S)
      H0P=DCMPLX(B0P,-Y0P)
      H1S=DCMPLX(B1S,-Y1S)
      H1P=DCMPLX(B1P,-Y1P)
      H2S=H1S*2.D0/AS-H0S
      H2P=H1P*2.D0/AP-H0P
C
      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
C
      CCC1=CC*ST*SW*CDEXP(DCMPLX(0.D0,-XT-XW))
      UF=CCC1*FX
      WF=CCC1*FZ
C
  100 CONTINUE
C
      RETURN
      END
