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