Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Geodesic Distance in PHP
function HowFar( $p1, $p2 ){
// WGS-84 Ellipsoid
$a = 6378137;
$b = 6356752.3142;
$f = 1/289.257223563;
$L = $p2->GetLong() - $p1->GetLong();
$U1 = atan( ( 1 - $f ) * tan( $p1->GetLat() ) );
$U2 = atan( ( 1 - $f ) * tan( $p2->GetLat() ) );
//print “L:” . $L . “\n”;
//print “U1:” . $U1 . “\n”;
//print “U2:” . $U2 . “\n”;
$sinU1 = sin( $U1 );
$cosU1 = cos( $U1 );
//print “sinU1:” . $sinU1 . “\n”;
//print “cosU1:” . $cosU1 . “\n”;
$sinU2 = sin( $U2 );
$cosU2 = cos( $U2 );
//print “sinU2:” . $sinU2 . “\n”;
//print “cosU2:” . $cosU2 . “\n”;
$lambda = $L;
$lambdaP = 2 * pi();
$iterLimit = 20;
while( abs( $lambda - $lambdaP ) > 0000000000001 && –$iterLimit > 0 ){
$sinLambda = sin( $lambda );
//print “sinLambda:” . $sinLambda . “\n”;
$cosLambda = cos( $lambda );
//print “cosLambda:” . $cosLambda . “\n”;
$sinSigma = sqrt( ( $cosU2 * $sinLambda ) * ( $cosU2 * $sinLambda ) + ( $cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosLambda ) * ( $cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosLambda ) );
//print “sinSigma:” . $sinSigma . “\n”;
if( $sinSigma == 0 ){
return 0; //co-incident points
}
$cosSigma = $sinU1 * $sinU2 + $cosU1 * $cosU2 * $cosLambda;
//print “cosSigma:” . $cosSigma . “\n”;
$sigma = atan2( $sinSigma, $cosSigma );
//print “sigma:” . $sigma . “\n”;
$sinAlpha = $cosU1 * $cosU2 * $sinLambda / $sinSigma;
//print “sinAlpha:” . $sinAlpha . “\n”;
$cosSqAlpha = 1 - $sinAlpha * $sinAlpha;
//print “cosSqAlpha:” . $cosSqAlpha . “\n”;
$cos2SigmaM = $cosSigma - 2 * $sinU1 * $sinU2 / $cosSqAlpha;
//print “cos2SigmaM:” . $cos2SigmaM . “\n”;
if( !is_numeric( $cosSigmaM ) ){
$cosSigmaM = 0; // Equatorial line
}
$C = $f/16 * $cosSqAlpha * ( 4 + $f * ( 4 - 3 * $cosSqAlpha ) );
$lambdaP = $lambda;
$lambda = $L + ( 1 - $C ) * $f * $sinAlpha * ( $sigma + $C * $sinSigma * ( $cos2SigmaM + $C * $cosSigma * ( -1 + 2 * $cos2SigmaM * $cos2SigmaM ) ) );
}
if( $iterLimit == 0 ){
return false;
}
$uSq = $cosSqAlpha * ( $a * $a - $b * $b ) / ( $b * $b );
//print “uSq:” . $uSq . “\n”;
$A = 1 + $uSq / 16384 * ( 4096 + $uSq * ( -768 + $uSq * ( 320 - 175 * $uSq ) ) );
//print “A:” . $A. “\n”;
$B = $uSq / 1024 * ( 256 + $uSq * ( -128 + $uSq * ( 74 - 47 * $uSq ) ) );
//print “B:” . $B. “\n”;
$deltaSigma = $B * $sinSigma * ( $cos2SigmaM + $B / 4 * ( $cosSigma * ( -1 + 2 * $cos2SigmaM * $cos2SigmaM ) -
$B / 6 * $cos2SigmaM * ( -3 + 4 * $sinSigma * $sinSigma ) * ( -3 + 4 * $cos2SigmaM * $cos2SigmaM ) ) );
$s = $b * $A * ( $sigma - $deltaSigma );
return $s;
}
class CPoint{
var $m_nLat;
var $m_nLong;
var $m_sTs;
function CPoint( $lat, $long, $ts = “” ){
$this->m_nLat = $lat;
$this->m_nLong = $long;
$this->m_sTs = strtotime( $ts );
}
function GetLat(){
return $this->m_nLat;
}
function GetLong(){
return $this->m_nLong;
}
function GetTs(){
return $this->m_sTs;
}
function Degrees2Radian(){
$this->m_nLat = deg2rad( $this->m_nLat );
$this->m_nLong = deg2rad( $this->m_nLong );
}
}
$Start = new CPoint( 33.399518681, -115.4098444488 );
$Start->Degrees2Radian();
$End = new CPoint( 33.399843764, -115.4054863956 );
$End->Degrees2Radian();
print HowFar( $Start, $End );
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.