# Great Circle Distance question

Posted on

Great Circle Distance question – Here in this article, we will share some of the most common and frequently asked about PHP problem in programming with detailed answers and code samples. There’s nothing quite so frustrating as being faced with PHP errors and being unable to figure out what is preventing your website from functioning as it should like php and latitude-longitude . If you have an existing PHP-based website or application that is experiencing performance issues, let’s get thinking about Great Circle Distance question.

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. ðŸ˜‰

Solution :

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

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).