Great Circle Distance question

To answer my own question just so it's here for anyone curious, a PHP class as converted from the C function provided by Chad Birch:

class GreatCircle {    /*     * Find a point a certain distance and vector away from an initial point     * converted from c function found at: http://sam.ucsd.edu/sio210/propseawater/ppsw_c/gcdist.c     *      * @param int distance in meters     * @param double direction in degrees i.e. 0 = North, 90 = East, etc.     * @param double lon starting longitude     * @param double lat starting latitude     * @return array ('lon' => $lon, 'lat' => $lat)     */    public static function getPositionByDistance($distance, $direction, $lon, $lat)    {        $metersPerDegree = 111120.00071117;        $degreesPerMeter = 1.0 / $metersPerDegree;        $radiansPerDegree = pi() / 180.0;        $degreesPerRadian = 180.0 / pi();        if ($distance > $metersPerDegree*180)        {            $direction -= 180.0;            if ($direction < 0.0)            {                $direction += 360.0;            }            $distance = $metersPerDegree * 360.0 - $distance;        }        if ($direction > 180.0)        {            $direction -= 360.0;        }        $c = $direction * $radiansPerDegree;        $d = $distance * $degreesPerMeter * $radiansPerDegree;        $L1 = $lat * $radiansPerDegree;        $lon *= $radiansPerDegree;        $coL1 = (90.0 - $lat) * $radiansPerDegree;        $coL2 = self::ahav(self::hav($c) / (self::sec($L1) * self::csc($d)) + self::hav($d - $coL1));        $L2   = (pi() / 2) - $coL2;        $l    = $L2 - $L1;        $dLo = (cos($L1) * cos($L2));        if ($dLo != 0.0)        {            $dLo  = self::ahav((self::hav($d) - self::hav($l)) / $dLo);        }        if ($c < 0.0)         {            $dLo = -$dLo;        }        $lon += $dLo;        if ($lon < -pi())        {            $lon += 2 * pi();        }        elseif ($lon > pi())        {            $lon -= 2 * pi();        }        $xlat = $L2 * $degreesPerRadian;        $xlon = $lon * $degreesPerRadian;        return array('lon' => $xlon, 'lat' => $xlat);    }    /*     * copy the sign     */    private static function copysign($x, $y)    {        return ((($y) < 0.0) ? - abs($x) : abs($x));    }       /*     * not greater than 1     */    private static function ngt1($x)    {        return (abs($x) > 1.0 ? self::copysign(1.0 , $x) : ($x));    }       /*     * haversine     */    private static function hav($x)    {        return ((1.0 - cos($x)) * 0.5);    }    /*     * arc haversine     */    private static function ahav($x)    {        return acos(self::ngt1(1.0 - ($x * 2.0)));    }    /*     * secant     */    private static function sec($x)    {        return (1.0 / cos($x));    }    /*     * cosecant     */    private static function csc($x)    {        return (1.0 / sin($x));    }}

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).*/voidGreatCirclePos(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

It would be harder to back out the Haversine fomula, then generate your own, I would think.

First the angle generated from the Earth core by traveling a "straight" line on the surface (you think it is straight, but it is curving).

Angle in Radians = Arc Length / radius.Angle = ArcLen/6371 km

Latitude should be easy, just the "vertical" (north/south) component of your angle.

Lat1 + Cos(bearing) * Angle

Longitude divisions vary by latitude. So that becomes harder. You would use:

Sin(bearing) * Angle (with East defined as negative) to find the angle in longitude direction, but converting back to actual longitude at that latitude would be more difficult.