# Great Circle Distance question

13.2k Views

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

i.e.

``````<?php
\$theta = \$lon1 - \$lon2;
\$dist = acos(\$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. ;) • 2
• Great question and great answers!!! Thanks guys!!

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;

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;
}

\$d = \$distance * \$degreesPerMeter * \$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();
}

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 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;
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;

} /* GreatCirclePos() */
``````

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

• 1
• Thanks very much, I converted this to PHP and it works perfectly. I will post the PHP code below.

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.