; ======================================================
  function distance,lat1,lon1,lat2,lon2
; ======================================================
; Calculates the distance between the two points 
;    [lat1,lon1] and [lat2,lon2] in km.
; ------------------------------------------------------

re=40000./2./!pi
dtr=double(!pi/180.)

dlat=abs(lat1-lat2)*dtr
dlon=abs(lon1-lon2)*dtr
avglat=.5*(lat1+lat2)*dtr

dxf=re*dlon*cos(avglat) 
dyf=re*dlat
d=sqrt(dxf*dxf+dyf*dyf)

return,d
end