program distance1
      real RS2,LS2,WS2,US2
      real xr(3),xr0(3),xl(3),xw(3),xu(3),a(3),b(3),xa(3),xb(3),
     *     xs(3),xh(3)
      character filen*80
      data eps/1.e-6/
      pai=4.*atan(1.0)
      pu=pai/180.
c Four verteces of earthquake fault
c xr(3):coordinates of the refence point
c xl(3):   "        of the vertex along strike direction
c xw(3):   "        of the vertex along dip direction
c xu(3):   "        of the fourth vertex
c
c xs(3):coordinates of the observation point.
c
c outer parameter of the fault: xr(3)(lat(deg),lon(deg),dep(km)),
c                                   length, width, strike, dip
        filen='distance1.prm'
        open(1,file=filen,status='old')
        read(1,*)(xr(i),i=1,3)
	read(1,*)al,aw
	read(1,*)strike,dip
c observation point (lat(deg),lon(deg))
        read(1,*)(xs(i),i=1,2)
	close(1)
	write(6,*)' Ref. Point    :',(xr(i),i=1,3)
	write(6,*)' Length & Width:',al,aw
	write(6,*)' Strike & Dip  :',strike, dip
	write(6,*)' Obs. Point    :',(xs(i),i=1,2)
        xs(3)=0.0                                 ! at the surface
c origin at xr(1:2) & at the surface, x- toward east, y- toward north
        if(abs((xs(1)-xr(1))**2+(xs(2)-xr(2))**2).gt.eps) then
         call distazi(xs(1),xs(2),xr(1),xr(2),dis_st,azi)
	else
	 xs(1)=xr(1)
	 xs(2)=xr(2)
	 dis_st=0.0
	 azi=0.0
	endif
c        write(6,*)xs(1),xs(2),xr(1),xr(2),dis_st,azi
        xr(1)=0.
        xr(2)=0.
c
        xs(1)=dis_st*sin(azi*pu)
        xs(2)=dis_st*cos(azi*pu)
c calculation of the four verteces of the fault plane
        xl(1)=xr(1)+al*sin(strike*pu)
        xl(2)=xr(2)+al*cos(strike*pu)
        xl(3)=xr(3)
c
        xw(1)=xr(1)+aw*cos(dip*pu)*cos(strike*pu)
        xw(2)=xr(2)-aw*cos(dip*pu)*sin(strike*pu)
        xw(3)=xr(3)+aw*sin(dip*pu)
c
        xu(1)=xw(1)+al*sin(strike*pu)
        xu(2)=xw(2)+al*cos(strike*pu)
        xu(3)=xw(3)
c
c        write(6,'(a5,3f8.1,a1)')' xs=(',(xs(i),i=1,3),')'
c        write(6,'(a5,3f6.1,a1)')' xr=(',(xr(i),i=1,3),')'
c        write(6,'(a5,3f6.1,a1)')' xl=(',(xl(i),i=1,3),')'
c        write(6,'(a5,3f6.1,a1)')' xw=(',(xw(i),i=1,3),')'
c        write(6,'(a5,3f6.1,a1)')' xu=(',(xu(i),i=1,3),')'
c+++++++++++calculation of ruture and JB distances++++++++++++++++
c
c (1)The shortest distance to the fault plane (Rupture distance or
c    fault distance)
c
c 1) the distance to the four verteces
      RS2=(xr(1)-xs(1))**2+(xr(2)-xs(2))**2+(xr(3)-xs(3))**2
      LS2=(xl(1)-xs(1))**2+(xl(2)-xs(2))**2+(xl(3)-xs(3))**2
      WS2=(xw(1)-xs(1))**2+(xw(2)-xs(2))**2+(xw(3)-xs(3))**2
      US2=(xu(1)-xs(1))**2+(xu(2)-xs(2))**2+(xu(3)-xs(3))**2
c 2) the shortest distance to the coseismic rupture
      tt=0.5-(RS2-LS2)/(2.*al*al)
      ss=0.5-(WS2-RS2)/(2.*aw*aw)
c 
      if(ss.lt.0.) then
       if(tt.gt.1.) then
        disrup=sqrt(RS2)
       else if(tt.ge.0..and.tt.le.1.) then
        disrup=sqrt(LS2-(al*tt)**2)
       else
        disrup=sqrt(LS2)
       endif
      else if(ss.ge.0..and.ss.le.1.) then
       if(tt.gt.1.) then
        disrup=sqrt(RS2-(aw*ss)**2)
       else if(tt.ge.0..and.tt.le.1.) then
        do i=1,3
         a(i)=xl(i)-xr(i)
         b(i)=xw(i)-xr(i)
        enddo
        aa=a(2)*b(3)-a(3)*b(2)
        bb=a(3)*b(1)-a(1)*b(3)
        cc=a(1)*b(2)-a(2)*b(1)
        dd=-(aa*xr(1)+bb*xr(2)+cc*xr(3))
        disrup=abs(aa*xs(1)+bb*xs(2)+cc*xs(3)+dd)
     *                   /sqrt(aa**2+bb**2+cc**2)               
       else
        disrup=sqrt(LS2-(aw*ss)**2)
       endif
      else
       if(tt.gt.1.) then
        disrup=sqrt(WS2)
       else if(tt.ge.0..and.tt.le.1.) then
        disrup=sqrt(US2-(al*tt)**2)
       else
        disrup=sqrt(US2)
       endif
      endif
c(2) The shortest distance to the surface projection of the fault plane 
c    (Joyner-Boore distance)
c 1) the distance to the four verteces along the surface
      RS2=(xr(1)-xs(1))**2+(xr(2)-xs(2))**2
      LS2=(xl(1)-xs(1))**2+(xl(2)-xs(2))**2
      WS2=(xw(1)-xs(1))**2+(xw(2)-xs(2))**2
      US2=(xu(1)-xs(1))**2+(xu(2)-xs(2))**2
c 2) the shortest distance to the coseismic rupture
      tt=0.5-(RS2-LS2)/(2.*al*al)
      ss=0.5-(US2-LS2)/(2.*(aw*cos(dip*pu))**2)
c 
      if(ss.lt.0.) then
       if(tt.gt.1.) then
        dis_jb=sqrt(RS2)
       else if(tt.ge.0..and.tt.le.1.) then
        dis_jb=sqrt(LS2-(al*tt)**2)
       else
        dis_jb=sqrt(LS2)
       endif
      else if(ss.ge.0..and.ss.le.1.) then
       if(tt.gt.1.) then
        dis_jb=sqrt(RS2-(aw*cos(dip*pu)*ss)**2)
       else if(tt.ge.0..and.tt.le.1.) then
        dis_jb=0.0              
       else
        dis_jb=sqrt(LS2-(aw*cos(dip*pu)*ss)**2)
       endif
      else
       if(tt.gt.1.) then
        dis_jb=sqrt(WS2)
       else if(tt.ge.0..and.tt.le.1.) then
        dis_jb=sqrt(US2-(al*tt)**2)
       else
        dis_jb=sqrt(US2)
       endif
      endif
      write(6,*)'Rupture  (Fault) dist.=',disrup,
     *          'Joyner - Boore_s dist.=',dis_jb
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      stop
      end
c
      subroutine distazi(evlar,evlor,stlar,stlor,delta,azi)
c delta:distance in Km
c azi  :azimuth from (stlar stlor) to (evlar evlor) in degree
      pai=4.*atan(1.0)
      pu=pai/180.
c coordinates of event
      evla=evlar*pu
      evlo=evlor*pu
c coordinates of station
      stla=stlar*pu
      stlo=stlor*pu
c conversion from geographic latitude to geocentric latitude
      evla=evla-0.00335975881*sin(2.*evla)
      stla=stla-0.00335975881*sin(2.*stla)
c
      ae=cos(evla)*cos(evlo)
      as=cos(stla)*cos(stlo)
      be=cos(evla)*sin(evlo)
      bs=cos(stla)*sin(stlo)
      ce=sin(evla)
      cs=sin(stla)
      a=sqrt((ae-as)**2+(be-bs)**2+(ce-cs)**2)/2.0
      at=atan(a/sqrt(-a*a+1.0))*2.0
c distance
      delta=at*6369.0
c azimuth from station to event in degree
      as=sin(evlo-stlo)*cos(evla)/sin(at)
      ac=(sin(evla)-sin(stla)*cos(at))/(cos(stla)*sin(at))
      azi=atan(abs(as/ac))
      if(as.gt.0.and.ac.ge.0.) then
      else if(as.le.0.and.ac.gt.0.) then
       azi=azi+pai/2.
      else if(as.lt.0.and.ac.le.0.) then
       azi=azi+pai
      else if(as.ge.0.and.ac.lt.0.) then
       azi=azi+3.*pai/2.
      endif
       azi=atan2(as,ac)
      if(azi.lt.0.) azi=azi+2.*pai
      azi=azi/pu
      return
      end