  function rng = compass(lat0,lon0,lat1,lon1);
  %
  % function rng = compass(lat0,lon0,lat1,lon1);
  %
  % ROUTINE:     compass
  %
  % PURPOSE:     Given the geographical coordinates of an array of points
  %             on the on the earth's surface, compute the great-circle 
  %             distance and compass heading from the central point.
  %             Distances and azimuth angles about the central point are 
  %             provided in KM and Degrees.  
  %
  % INPUT:
  %  lat0       geographic latitude (degrees) of point(s) 1
  %  lon0       geographic longitude (degrees) of point(s) 1
  %  lat1       geographic latitudes (degrees) of point(s) 2
  %  lon1       geographic longitudes (degrees) of point(s) 2
  %
  % OUTPUT
  %  rng        great circle distances (km) from (lat0,lon0) to (lat1,lon1)
  %
  %  author:  Paul Ricchiazzi                            jan93
  %           Institute for Computational Earth System Science
  %           University of California, Santa Barbara
  %
  % re-written in matlab. Dave Tobin, CIMSS/SSEC 7/10/03
  %
  re = 6371.2;
  dtor = pi/180;

  t0 = (90-lat0)*dtor;
  t1 = (90-lat1)*dtor;
  p0 = lon0*dtor;
  p1 = lon1*dtor;
  zz = cos(t0).*cos(t1)+sin(t0).*sin(t1).*cos(p1-p0);   
  rng = re*acos(zz);


