-
Notifications
You must be signed in to change notification settings - Fork 5
/
VincentyFormula.php
117 lines (73 loc) · 3.22 KB
/
VincentyFormula.php
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
<?php
namespace Geodesy\Distance;
use Geodesy\Location\LatLong;
/**
* 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 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
*/
class VincentyFormula extends BaseDistance implements DistanceInterface {
public function __construct(LatLong $source, LatLong $destination) {
parent::__construct($source, $destination);
}
protected function distance(): float{
$lat1 = $this->getSourceLatitude();
$lat2 = $this->getDestinationLatitude();
$long1 = $this->getSourceLongitude();
$long2 = $this->getDestinationLongitude();
$a = $this->getSemiMajorAxis();
$b = $this->getSemiMinorAxis();
$f = $this->getInverseFlattening();
$L = $long2 - $long1;
$tanU1 = (1 - $f) * tan($lat1);
$tanU2 = (1 - $f) * tan($lat2);
$cosU1 = 1 / sqrt((1 + pow($tanU1, 2)));
$sinU1 = $tanU1 * $cosU1;
$cosU2 = 1 / sqrt((1 + pow($tanU2, 2)));
$sinU2 = $tanU2 * $cosU2;
$lambda = $L;
$lambdaPrime = 0;
$iterations = 1000;
$antimeridian = abs($L) > pi();
do {
$sinLambda = sin($lambda);
$cosLambda = cos($lambda);
$sinSqSigma = pow($cosU2 * $sinLambda, 2) + pow($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosLambda, 2);
if ($sinSqSigma == 0) {
break;
}
$sinSigma = sqrt($sinSqSigma);
$cosSigma = $sinU1 * $sinU2 + $cosU1 * $cosU2 * $cosLambda;
$sigma = atan2($sinSigma, $cosSigma);
$sinAlpha = $cosU1 * $cosU2 * $sinLambda / $sinSigma;
$cosSqAlpha = 1 - pow($sinAlpha, 2);
$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 - $lambdaPrime) > 1e-12 && --$iterations >= 0);
if ($iterations == 0) {
throw new \Exception('Failed to converge');
}
$uSq = $cosSqAlpha * (pow($a / $b, 2) - 1);
/**
we'll use Helmert's expansion for A and B
$k1 = (sqrt(1 + $uSq) - 1) / (sqrt(1 + $uSq) + 1);
$A = 1 + ((1/4) * pow($k1, 2)) / (1 - $k1);
$B = $k1 * (1 - ((3/8) * pow($k1, 2)));
*/
$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 * pow($cos2SigmaM, 2)) -
$B / 6 * $cos2SigmaM * (-3 + 4 * pow($sinSigma, 2)) * (-3 + 4 * pow($cos2SigmaM, 2))));
$s = $b * $A * ($sigma - $deltaSigma);
return round($s, 3);
}
}