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