function d = distsphr(xlon1,xlat1,xlon2,xlat2)

      x1 = cos(xlon1)*cos(xlat1);
      y1 = sin(xlon1)*cos(xlat1);
      z1 = sin(xlat1);

      x2 = cos(xlon2)*cos(xlat2);
      y2 = sin(xlon2)*cos(xlat2);
      z2 = sin(xlat2);

      alpha = acos( (x1*x2) + (y1*y2) + (z1*z2) );
      d = abs(alpha);
