#include <math.h>

// computes the distance between points on earth's surface
double dist_earth(double latitude1, double longitude1,
                  double latitude2, double longitude2)  
{
  const double circ_earth = 24901.5;  // miles (at equator)
  double theta1 = longitude1 * M_PI / 180,
           phi1 = (90-latitude1) * M_PI / 180,
             x1 = cos(theta1) * sin(phi1),
             y1 = sin(theta1) * sin(phi1),
             z1 = cos(phi1),

         theta2 = longitude2 * M_PI / 180,
           phi2 = (90-latitude2) * M_PI / 180,
             x2 = cos(theta2) * sin(phi2),
             y2 = sin(theta2) * sin(phi2),
             z2 = cos(phi2),

              s = sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1)),
          alpha = 2 * asin(s/2),
         answer = circ_earth * alpha / (2 * M_PI);

  return answer;
}