Skip to content

Commit

Permalink
added Hubeny, Andoyer-Lambert and Thomas Formulas. Fixed Vincenty.
Browse files Browse the repository at this point in the history
  • Loading branch information
jtejido committed Apr 27, 2018
1 parent ff9353f commit 7847f02
Show file tree
Hide file tree
Showing 14 changed files with 228 additions and 79 deletions.
15 changes: 11 additions & 4 deletions src/Geodesy/Constants/Constants.php
Original file line number Diff line number Diff line change
Expand Up @@ -6,12 +6,19 @@ class Constants
{

// Since earth has been considered spherical in most (Law of Cosine and Haversine) formulas here, we'll have the Earth's mean radius
CONST SPHERICAL_R = 6371.000;
CONST SPHERICAL_R = 6371000;

// WGS-84 semi-major axis
CONST WGS_R = 6378.137;
CONST A = 6378137;

// WGS-84 semi-minor axis
CONST B = 6356752.314245;

// WGS-84 flattening
CONST F = 1 / 298.257223563;

// WGS-84 eccentricity
CONST ECCENTRICITY = 0.0818191908426215;

// WGS-84 first eccentricity
CONST E = 8.1819190842622e-2;

}
7 changes: 3 additions & 4 deletions src/Geodesy/Conversion/BaseConversion.php
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@

use Geodesy\Constants\Constants;
use Geodesy\Unit\UnitInterface;
use Geodesy\Unit\KiloMetre;
use Geodesy\Unit\Metre;

abstract class BaseConversion
{
Expand All @@ -15,16 +15,15 @@ abstract class BaseConversion

/**
* https://en.wikipedia.org/wiki/Geographic_coordinate_conversion
* There have been known methods for two-way conversion but what we use is Olson's method as this is computationally inexpensive
* and pretty much accurate for majority of cases(except, of course, if you aim for ultimate precision).
* There have been known methods for two-way conversion but we'll use Olson's method as this is computationally inexpensive.
* Olson, D. K. (1996). "Converting earth-Centered, Earth-Fixed Coordinates to Geodetic Coordinates,"
* IEEE Transactions on Aerospace and Electronic Systems, Vol. 32, No. 1, January 1996, pp. 473-476
*/

public function __construct()
{
$this->constants = new Constants;
$this->unit = new KiloMetre;
$this->unit = new Metre;
}

public abstract function convert();
Expand Down
10 changes: 5 additions & 5 deletions src/Geodesy/Conversion/ECEF2LLA.php
Original file line number Diff line number Diff line change
Expand Up @@ -25,17 +25,17 @@ public function convert(){
$y = $this->ecef->getY();
$z = $this->ecef->getZ();

$esq = pow($this->constants::E, 2);
$asq = pow($this->constants::WGS_R, 2);
$esq = pow($this->constants::ECCENTRICITY, 2);
$asq = pow($this->constants::A, 2);
$b = sqrt( $asq * (1 - $esq) );
$bsq = pow($b, 2);
$ep = sqrt( ($asq - $bsq) / $bsq );
$p = sqrt( pow($x, 2) + pow($y, 2) );
$th = atan2($this->constants::WGS_R * $z, $b * $p);
$th = atan2($this->constants::A * $z, $b * $p);

$long = atan2($y,$x);
$lat = atan2( ($z + pow($ep,2) * $b * pow(sin($th), 3) ), ($p - $esq * $this->constants::WGS_R * pow(cos($th),3)) );
$N = $this->constants::WGS_R / ( sqrt(1 - $esq * pow(sin($lat),2)) );
$lat = atan2( ($z + pow($ep,2) * $b * pow(sin($th), 3) ), ($p - $esq * $this->constants::A * pow(cos($th),3)) );
$N = $this->constants::A / ( sqrt(1 - $esq * pow(sin($lat),2)) );
$alt = $p / cos($lat) - $N;

$long = ($long % (2*pi())) == 0 ? $long : $long % (2*pi());
Expand Down
4 changes: 2 additions & 2 deletions src/Geodesy/Conversion/LLA2ECEF.php
Original file line number Diff line number Diff line change
Expand Up @@ -27,9 +27,9 @@ public function convert(){

$alt = $this->latlong->getAltitude();

$esq = pow($this->constants::E, 2);
$esq = pow($this->constants::ECCENTRICITY, 2);

$n = $this->constants::WGS_R / sqrt( 1 - $esq * pow(sin( $lat ), 2) );
$n = $this->constants::A / sqrt( 1 - $esq * pow(sin( $lat ), 2) );

$this->ecef->setX($this->getUnit()->convert(( $n + $alt ) * cos( $lat ) * cos( $long ))); //ECEF x

Expand Down
39 changes: 39 additions & 0 deletions src/Geodesy/Distance/AndoyerLambert.php
Original file line number Diff line number Diff line change
@@ -0,0 +1,39 @@
<?php

namespace Geodesy\Distance;

use Geodesy\Location\LatLong;

class AndoyerLambert extends BaseDistance implements DistanceInterface
{


public function __construct(LatLong $source, LatLong $destination)
{
parent::__construct($source, $destination);
}


public function distance()
{
$lat1 = deg2rad($this->source->getLatitude());
$lat2 = deg2rad($this->destination->getLatitude());
$long1 = deg2rad($this->source->getLongitude());
$long2 = deg2rad($this->destination->getLongitude());
$a = $this->constants::A;
$b = $this->constants::B;
$f = $this->constants::F;
$long_diff = $long2 - $long1;
$rLat1=atan($b*tan($lat1)/$a);
$rLat2=atan($b*tan($lat2)/$a);
$spherical_distance = acos(sin($rLat1) * sin($rLat2) + cos($rLat1) * cos($rLat2) * cos($long_diff));
$cosineSd = cos($spherical_distance/2.0);
$sinSd = sin($spherical_distance/2.0);
$c = (sin($spherical_distance) - $spherical_distance) * (sin($rLat1)+sin($rLat2))*(sin($rLat1)+sin($rLat2)) / pow($cosineSd, 2);
$s = (sin($spherical_distance) + $spherical_distance) * (sin($rLat1)-sin($rLat2))*(sin($rLat1)-sin($rLat2)) / pow($sinSd,2);
$delta = $f / 8.0 * ($c-$s);

return $a*($spherical_distance+$delta);
}

}
4 changes: 2 additions & 2 deletions src/Geodesy/Distance/BaseDistance.php
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@

use Geodesy\Location\LatLong;
use Geodesy\Unit\UnitInterface;
use Geodesy\Unit\KiloMetre;
use Geodesy\Unit\Metre;
use Geodesy\Constants\Constants;

abstract class BaseDistance
Expand All @@ -22,7 +22,7 @@ public function __construct(LatLong $source, LatLong $destination)
{
$this->source = $source;
$this->destination = $destination;
$this->unit = new KiloMetre;
$this->unit = new Metre;
$this->constants = new Constants;
}

Expand Down
48 changes: 48 additions & 0 deletions src/Geodesy/Distance/HubenyFormula.php
Original file line number Diff line number Diff line change
@@ -0,0 +1,48 @@
<?php

namespace Geodesy\Distance;

use Geodesy\Location\LatLong;

class HubenyFormula extends BaseDistance implements DistanceInterface
{


public function __construct(LatLong $source, LatLong $destination)
{
parent::__construct($source, $destination);
}


public function distance()
{
$lat1 = deg2rad($this->source->getLatitude());
$lat2 = deg2rad($this->destination->getLatitude());
$long1 = deg2rad($this->source->getLongitude());
$long2 = deg2rad($this->destination->getLongitude());
$a = $this->constants::A;
$b = $this->constants::B;
$f = $this->constants::F;
$f2 = ($b * $b) / ($a * $a);
$e2 = 1.0 - $f2;

$lat_diff = ($lat1 - $lat2);
$long_diff = ($long1 - $long2);
$lat_average = 0.5 * ($lat1 + $lat2);
$sinlat_average = sin($lat_average);
$coslat_average = cos($lat_average);
$w2 = 1.0 - $sinlat_average * $sinlat_average * $e2;
$w = sqrt($w2);
$meridian = $a * $f2 / ($w2 * $w);
$n = $a / $w;

return sqrt(
$lat_diff * $lat_diff * $meridian * $meridian +
$long_diff * $long_diff * $n * $n * $coslat_average * $coslat_average
);


return $a*($spherical_distance+$delta);
}

}
61 changes: 61 additions & 0 deletions src/Geodesy/Distance/ThomasFormula.php
Original file line number Diff line number Diff line change
@@ -0,0 +1,61 @@
<?php

namespace Geodesy\Distance;

use Geodesy\Location\LatLong;

class ThomasFormula extends BaseDistance implements DistanceInterface
{


public function __construct(LatLong $source, LatLong $destination)
{
parent::__construct($source, $destination);
}


public function distance()
{

$A = $this->constants::A;
$F = $this->constants::F;
$ac = deg2rad($this->source->getLatitude());
$ad = -deg2rad($this->source->getLongitude());
$ac = atan((1.0-$F) * tan($ac));

$ae = deg2rad($this->destination->getLatitude());
$af = -deg2rad($this->destination->getLongitude());
$ae = atan((1.0-$F) * tan($ae));


$i = $af - $ad;
$j = ($ae + $ac) / 2.0;
$ak = ($ae - $ac) / 2.0;
$h = cos($ak) * cos($ak) - sin($j) * sin($j);

$l = sin($ak) * sin($ak) + $h * sin($i/2.0) * sin($i/2.0);
$d = 2.0 * atan(sqrt($l/(1.0-$l)));

$u = 2.0 * sin($j) * sin($j) * cos($ak) * cos($ak) / (1.0 - $l);
$v = 2.0 * sin($ak) * sin($ak) * cos($j) * cos($j) / $l;

$x = $u + $v;
$y = $u - $v;
$t = $d / sin($d);
$dd = 4.0 * $t * $t;
$e =2.0 * cos($d);
$a = $dd * $e;

$c = $t - ($a - $e) / 2.0;
$n1 = $x * ($a + $c * $x);
$b = 2.0 * $dd;
$n2 = $y * ($b + $e * $y);
$n3 = $dd * $x * $y;

$d1d=$F * ($t * $x - $y)/4.0;
$d2d=$F * $F * ($n1 - $n2 + $n3)/64.0;

return $A * ($t - $d1d + $d2d) * sin($d);
}

}
106 changes: 49 additions & 57 deletions src/Geodesy/Distance/VincentyFormula.php
Original file line number Diff line number Diff line change
Expand Up @@ -15,69 +15,63 @@ public function __construct(LatLong $source, LatLong $destination)
/**
* Vincenty's formula is using an accurate ellipsoidal model of the earth as opposed to a pressumably perfectly spherical shape.
* Vincenty’s solution for the distance between points on an ellipsoidal earth model is accurate to within a millimeter.
* While it doesn't appear to take into account the seismic activity of the earth (yes, where we are right now isn't the same next month
* or year, and so on), it depends on standards to give out a precise latitude and longitude data (http://itrf.ign.fr/ITRF_solutions/2014)
* This is a long formula (we use the indirect):
* While it doesn't appear to take into account the seismic activity of the earth (yes, where we are right now isn't the same next year
* or decade, and so on), it depends on standards to give out a precise latitude and longitude data (http://itrf.ign.fr/ITRF_solutions/2014)
* This is a long formula:
* https://en.wikipedia.org/wiki/Vincenty%27s_formulae
*/
public function distance()
{

$lat1 = deg2rad($this->source->getLatitude());
$lat2 = deg2rad($this->destination->getLatitude());
$long1 = deg2rad($this->source->getLongitude());
$long2 = deg2rad($this->destination->getLongitude());
$a = 6378137;
$b = 6356752.314245;
$f = 1 / 298.257223563;
$a = $this->constants::A;
$b = $this->constants::B;
$f = $this->constants::F;
$L = $long2 - $long1;
$U1 = atan((1 - $f) * tan($lat1));
$U2 = atan((1 - $f) * tan($lat2));
$sinU1 = sin($U1);
$cosU1 = cos($U1);
$sinU2 = sin($U2);
$cosU2 = cos($U2);
$tanU1 = (1-$f) * tan($lat1);
$tanU2 = (1-$f) * tan($lat2);
$cosU1 = 1 / sqrt((1 + ($tanU1*$tanU1) ));
$sinU1 = $tanU1 * $cosU1;
$cosU2 = 1 / sqrt((1 + ($tanU2*$tanU2) ));
$sinU2 = $tanU2 * $cosU2;

$lambda = $L;
$iterLimit = 100;
$lambdaPrime = 0;
$iterations = 1000;
$antimeridian = abs($L) > pi();
do
{
$sinLambda = sin($lambda);
$cosLambda = cos($lambda);
$sinSigma = sqrt(pow(($cosU2 * $sinLambda), 2) + pow(($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosLambda),2));
if ($sinSigma == 0)
{
return 0;
}

$cosSigma = $sinU1 * $sinU2 + $cosU1 * $cosU2 * $cosLambda;
$sigma = atan2($sinSigma, $cosSigma);
$sinAlpha = ($cosU1 * $cosU2 * $sinLambda) / $sinSigma;
$cosSqAlpha = 1 - $sinAlpha * $sinAlpha;
$cos2SigmaM = $cosSigma - 2 * $sinU1 * $sinU2 / $cosSqAlpha;

if (is_nan($cos2SigmaM)) {
$cos2SigmaM = 0;
}

$C = $f / 16 * $cosSqAlpha * (4 + $f * (4 - 3 * $cosSqAlpha));
$lambda_prime = $lambda;
$lambda = $L + (1 - $C) * $f * $sinAlpha
* ($sigma + $C * $sinSigma
* ($cos2SigmaM + $C * $cosSigma
* (-1 + 2 * $cos2SigmaM * $cos2SigmaM)
)
);
$sinLambda = sin($lambda);
$cosLambda = cos($lambda);
$sinSqSigma = ($cosU2*$sinLambda) * ($cosU2*$sinLambda) + ($cosU1*$sinU2-$sinU1*$cosU2*$cosLambda) * ($cosU1*$sinU2-$sinU1*$cosU2*$cosLambda);
if ($sinSqSigma == 0) {
break;
}
$sinSigma = sqrt($sinSqSigma);
$cosSigma = $sinU1*$sinU2 + $cosU1*$cosU2*$cosLambda;
$sigma = atan2($sinSigma, $cosSigma);
$sinAlpha = $cosU1 * $cosU2 * $sinLambda / $sinSigma;
$cosSqAlpha = 1 - $sinAlpha*$sinAlpha;
$cos2SigmaM = ($cosSqAlpha != 0) ? ($cosSigma - 2*$sinU1*$sinU2/$cosSqAlpha) : 0;
$C = $f/16*$cosSqAlpha*(4+$f*(4-3*$cosSqAlpha));
$lambdaPrime = $lambda;
$lambda = $L + (1-$C) * $f * $sinAlpha * ($sigma + $C*$sinSigma*($cos2SigmaM+$C*$cosSigma*(-1+2*$cos2SigmaM*$cos2SigmaM)));
$iterationCheck = $antimeridian ? abs($lambda)-pi() : abs($lambda);
if ($iterationCheck > pi()) {
throw new \Exception('lambda is greater than pi');
}

} while (abs($lambda - $lambda_prime) > 1e12 && --$iterLimit > 0);
} while (abs($lambda-$lambdaPrime) > 1e-12 && --$iterations >= 0);

if ($iterLimit == 0)
if ($iterations == 0)
{
return 0; // Failed to converge
throw new \Exception('Failed to converge');
}

$uSq = $cosSqAlpha * (pow($a, 2)- pow($b, 2)) / ($b * $b);
$A = 1 + $uSq / 16384
* (4096 + $uSq * (-768 + $uSq * (320 - 175 * $uSq)));
$B = $uSq / 1024 * (256 + $uSq * (-128 + $uSq * (74 - 47 * $uSq)));
$uSq = $cosSqAlpha * ($a*$a - $b*$b) / ($b*$b);

/**
* Alternatively, you can use Helbert's expansion for A and B
Expand All @@ -87,17 +81,15 @@ public function distance()
* $B = $k1 * (1 - ((3/8) * pow($k1, 2)));
*/

$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 / 1000;
}
$A = 1 + $uSq/16384*(4096+$uSq*(-768+$uSq*(320-175*$uSq)));
$B = $uSq/1024 * (256+$uSq*(-128+$uSq*(74-47*$uSq)));
$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;
}
}
Loading

0 comments on commit 7847f02

Please sign in to comment.