Great Circle Distance question

Question

I am familiar with the formula to calculate the Great Circle Distance between two points.

i.e.

<?php
$theta = $lon1 - $lon2; 
$dist = sin(deg2rad($lat1)) * sin(deg2rad($lat2)) +  cos(deg2rad($lat1)) * cos(deg2rad($lat2)) * cos(deg2rad($theta)); 
$dist = acos($dist); 
$dist = rad2deg($dist); 
//convert degrees to distance depending on units desired
?>

What I need though, is the reverse of this. Given a starting point, a distance, and a simple cardinal NSEW direction, to calculate the position of the destination point. It's been a long time since I was in a math class. ;)

Answer 1

Here's a C implementation that I found, should be fairly straightforward to translate to PHP:

#define KmPerDegree         111.12000071117
#define DegreesPerKm        (1.0/KmPerDegree)
#define PI                  M_PI
#define TwoPI               (M_PI+M_PI)
#define HalfPI              M_PI_2
#define RadiansPerDegree    (PI/180.0)
#define DegreesPerRadian    (180.0/PI)
#define copysign(x,y)       (((y)<0.0)?-fabs(x):fabs(x))
#define NGT1(x)             (fabs(x)>1.0?copysign(1.0,x):(x))
#define ArcCos(x)           (fabs(x)>1?quiet_nan():acos(x))
#define hav(x)              ((1.0-cos(x))*0.5)              /* haversine */
#define ahav(x)             (ArcCos(NGT1(1.0-((x)*2.0))))   /* arc haversine */
#define sec(x)              (1.0/cos(x))                    /* secant */
#define csc(x)              (1.0/sin(x))                    /* cosecant */

/*
**  GreatCirclePos() --
**
**  Compute ending position from course and great-circle distance.
**
**  Given a starting latitude (decimal), the initial great-circle
**  course and a distance along the course track, compute the ending
**  position (decimal latitude and longitude).
**  This is the inverse function to GreatCircleDist).
*/
void
GreatCirclePos(dist, course, slt, slg, xlt, xlg)
    double  dist;   /* -> great-circle distance (km) */
    double  course; /* -> initial great-circle course (degrees) */
    double  slt;    /* -> starting decimal latitude (-S) */
    double  slg;    /* -> starting decimal longitude(-W) */
    double  *xlt;   /* <- ending decimal latitude (-S) */
    double  *xlg;   /* <- ending decimal longitude(-W) */
{
    double  c, d, dLo, L1, L2, coL1, coL2, l;

    if (dist > KmPerDegree*180.0) {
        course -= 180.0;
        if (course < 0.0) course += 360.0;
        dist    = KmPerDegree*360.0-dist;
    }
    if (course > 180.0) course -= 360.0;
    c    = course*RadiansPerDegree;
    d    = dist*DegreesPerKm*RadiansPerDegree;
    L1   = slt*RadiansPerDegree;
    slg *= RadiansPerDegree;
    coL1 = (90.0-slt)*RadiansPerDegree;
    coL2 = ahav(hav(c)/(sec(L1)*csc(d))+hav(d-coL1));
    L2   = HalfPI-coL2;
    l    = L2-L1;
    if ((dLo=(cos(L1)*cos(L2))) != 0.0)
        dLo  = ahav((hav(d)-hav(l))/dLo);
    if (c < 0.0) dLo = -dLo;
    slg += dLo;
    if (slg < -PI)
        slg += TwoPI;
    else if (slg > PI)
        slg -= TwoPI;

    *xlt = L2*DegreesPerRadian;
    *xlg = slg*DegreesPerRadian;

} /* GreatCirclePos() */

Source: http://sam.ucsd.edu/sio210/propseawater/ppsw_c/gcdist.c