earth surface distance

what is the shortest distance between two points on the earths surface?

given two points (long1, lat1) and (long2, lat), dist1 finds the great circle distance assuming the earth is a sphere whilst dist2 gives a fairly accurate approximation for the shortest earth surface distance taking flattening into account.

#include <math.h>
#include <stdio.h>

#define PI 3.141592653589793238462643383279
#define RAD (PI/180)

double dist1(double long1, double lat1, double long2, double lat2)
{
    /* calculate great circle distance */
    double l1 = lat1*RAD;
    double l2 = lat2*RAD;
    double d;
    const double a = 6371; /* mean earth radius */

    d = sin(l1)*sin(l2)+cos(l1)*cos(l2)*cos((long1-long2)*RAD);
    return a*acos(d);
}

double dist2(double long1, double lat1, double long2, double lat2)
{
    /* approx distance between points on earth ellipsoid.
     * H.Andoyer from Atronomical Algorithms, Jean Meeus, second edition.
     */
    double f = (lat1 + lat2)/2*RAD;
    double g = (lat1 - lat2)/2*RAD;
    double l = (long1 - long2)/2*RAD;

    double sg = sin(g);
    double sl = sin(l);
    double sf = sin(f);

    double s, c, w, r, d, h1, h2;
    const double a = 6378.14; /* equator earth radius */
    const double fl = 1/298.257;  /* earth flattening */

    sg = sg*sg;
    sl = sl*sl;
    sf = sf*sf;

    s = sg*(1-sl)+(1-sf)*sl;
    c = (1-sg)*(1-sl)+sf*sl;

    w = atan(sqrt(s/c));
    r = sqrt(s*c)/w;
    d = 2*w*a;
    h1 = (3*r-1)/2/c;
    h2 = (3*r+1)/2/s;

    return d*(1+fl*(h1*sf*(1-sg)-h2*(1-sf)*sg));
}

main()
{
    double long1, lat1;
    double long2, lat2;
    double d1, d2;

    printf("place1: "); scanf("%lf %lf", &long1, &lat1);
    printf("place2: "); scanf("%lf %lf", &long2, &lat2);

    d1 = dist1(long1, lat1, long2, lat2);
    d2 = dist2(long1, lat1, long2, lat2);

    printf("great circle dist = %fKm\n", d1);
    printf("ellipsoid dist = %fKm\n", d2);
    return 0;
}