Geodesic Distance in PHP

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