PRO dis1,lon1,lat1,lon2,lat2,r1 pi2 = !pi/2. rad = !pi/180. grd = 180./!pi degtokm = 111.195 ; Latitude of North Pole latnp=pi2 lon0 = lon1*rad loo1 = lon2*rad lat0 = lat1*rad laa1 = lat2*rad ; Determing which quadrant siglat = lat2-lat1 siglon = lon2-lon1 IF siglat GE 0 THEN BEGIN IF siglon GE 0 THEN BEGIN quad = 1 ENDIF ELSE BEGIN quad = 4 ENDELSE ENDIF ELSE BEGIN IF siglon GE 0 THEN BEGIN quad = 2 ENDIF ELSE BEGIN quad = 3 ENDELSE ENDELSE ;print,quad c1 = loo1-lon0 if (c1 ne 0.) THEN BEGIN c1=1./tan(c1/2.) lla1=laa1-lat0 llb1=laa1+lat0 ;print,lla1,llb1 if (lla1 NE 0.) THEN BEGIN yxa1=c1*sin(lla1/2.)/cos(llb1/2.) yxb1=c1*cos(lla1/2.)/sin(llb1/2.) yxa1=atan(yxa1) yxb1=atan(yxb1) r1=tan(lla1/2.)*sin(yxb1)/sin(yxa1) r1=atan(r1)*2. b1 = r1 a1 = latnp-lat0 c1 = latnp-laa1 s=(a1+b1+c1)/2. r=abs(sin(s-a1)*sin(s-b1)*sin(s-c1)/sin(s)) r=sqrt(r) gama=2.*atan(r/sin(s-c1)) gama=gama*grd zeta = gama if quad eq 2 then gama = 180 + gama if quad eq 1 then gama = 180 + gama ;if quad eq 3 then print,gama IF quad EQ 3 THEN gama = 180.-gama ;if quad eq 4 then print,gama IF quad EQ 4 THEN gama = 180.-gama ENDIF ELSE BEGIN ; Applying the Right Spherical Triangle rule lla1=latnp-laa1 c1 = loo1-lon0 r1=sin(c1)*sin(lla1) r1=asin(r1) IF siglon GE 0 THEN gama = 90. IF siglon LT 0 THEN gama = 270. zeta = gama endelse ENDIF ELSE BEGIN ;In the Latitude Great Circle r1=abs(laa1-lat0) IF siglat GE 0 THEN gama = 0. IF siglat LT 0 THEN gama = 180. zeta = gama ENDELSE ; Converting to Degrees r1=r1*grd ; Converting to Km (picking just the positive) r1=abs(r1*degtokm) END