|
1
|
|
|
<?php |
|
2
|
|
|
/** |
|
3
|
|
|
* Calculation of bearing between two points using a |
|
4
|
|
|
* ellipsoidal model of the earth |
|
5
|
|
|
* |
|
6
|
|
|
* This class is based on the awesome work Chris Veness |
|
7
|
|
|
* has done. For more information visit the following URL. |
|
8
|
|
|
* |
|
9
|
|
|
* @see http://www.movable-type.co.uk/scripts/latlong-vincenty.html |
|
10
|
|
|
* |
|
11
|
|
|
* @author Marcus Jaschen <[email protected]> |
|
12
|
|
|
* @license https://opensource.org/licenses/GPL-3.0 GPL |
|
13
|
|
|
* @link https://github.com/mjaschen/phpgeo |
|
14
|
|
|
*/ |
|
15
|
|
|
|
|
16
|
|
|
namespace Location\Bearing; |
|
17
|
|
|
|
|
18
|
|
|
use Location\Coordinate; |
|
19
|
|
|
use Location\Exception\NotConvergingException; |
|
20
|
|
|
|
|
21
|
|
|
/** |
|
22
|
|
|
* Calculation of bearing between two points using a |
|
23
|
|
|
* ellipsoidal model of the earth |
|
24
|
|
|
* |
|
25
|
|
|
* @author Marcus Jaschen <[email protected]> |
|
26
|
|
|
* @license https://opensource.org/licenses/GPL-3.0 GPL |
|
27
|
|
|
* @link https://github.com/mjaschen/phpgeo |
|
28
|
|
|
*/ |
|
29
|
|
|
class BearingEllipsoidal implements BearingInterface |
|
30
|
|
|
{ |
|
31
|
|
|
/** |
|
32
|
|
|
* This method calculates the initial bearing between the |
|
33
|
|
|
* two points. |
|
34
|
|
|
* |
|
35
|
|
|
* @param \Location\Coordinate $point1 |
|
36
|
|
|
* @param \Location\Coordinate $point2 |
|
37
|
|
|
* |
|
38
|
|
|
* @return float Bearing Angle |
|
39
|
|
|
*/ |
|
40
|
|
|
public function calculateBearing(Coordinate $point1, Coordinate $point2) |
|
41
|
|
|
{ |
|
42
|
|
|
return $this->inverseVincenty($point1, $point2)['bearing_initial']; |
|
43
|
|
|
} |
|
44
|
|
|
|
|
45
|
|
|
/** |
|
46
|
|
|
* Calculates the final bearing between the two points. |
|
47
|
|
|
* |
|
48
|
|
|
* @param \Location\Coordinate $point1 |
|
49
|
|
|
* @param \Location\Coordinate $point2 |
|
50
|
|
|
* |
|
51
|
|
|
* @return float |
|
52
|
|
|
*/ |
|
53
|
|
|
public function calculateFinalBearing(Coordinate $point1, Coordinate $point2) |
|
54
|
|
|
{ |
|
55
|
|
|
return $this->inverseVincenty($point1, $point2)['bearing_final']; |
|
56
|
|
|
} |
|
57
|
|
|
|
|
58
|
|
|
/** |
|
59
|
|
|
* Calculates a destination point for the given point, bearing angle, |
|
60
|
|
|
* and distance. |
|
61
|
|
|
* |
|
62
|
|
|
* @param \Location\Coordinate $point |
|
63
|
|
|
* @param float $bearing the bearing angle between 0 and 360 degrees |
|
64
|
|
|
* @param float $distance the distance to the destination point in meters |
|
65
|
|
|
* |
|
66
|
|
|
* @return Coordinate |
|
67
|
|
|
*/ |
|
68
|
|
|
public function calculateDestination(Coordinate $point, $bearing, $distance) |
|
69
|
|
|
{ |
|
70
|
|
|
return $this->directVincenty($point, $bearing, $distance)['destination']; |
|
71
|
|
|
} |
|
72
|
|
|
|
|
73
|
|
|
/** |
|
74
|
|
|
* Calculates the final bearing angle for a destination point. |
|
75
|
|
|
* The method expects a starting point point, the bearing angle, |
|
76
|
|
|
* and the distance to destination. |
|
77
|
|
|
* |
|
78
|
|
|
* @param \Location\Coordinate $point |
|
79
|
|
|
* @param float $bearing |
|
80
|
|
|
* @param float $distance |
|
81
|
|
|
* |
|
82
|
|
|
* @return float |
|
83
|
|
|
*/ |
|
84
|
|
|
public function calculateDestinationFinalBearing(Coordinate $point, $bearing, $distance) |
|
85
|
|
|
{ |
|
86
|
|
|
return $this->directVincenty($point, $bearing, $distance)['bearing_final']; |
|
87
|
|
|
} |
|
88
|
|
|
|
|
89
|
|
|
private function directVincenty(Coordinate $point, $bearing, $distance) |
|
90
|
|
|
{ |
|
91
|
|
|
$φ1 = deg2rad($point->getLat()); |
|
92
|
|
|
$λ1 = deg2rad($point->getLng()); |
|
93
|
|
|
$α1 = deg2rad($bearing); |
|
94
|
|
|
|
|
95
|
|
|
$a = $point->getEllipsoid()->getA(); |
|
96
|
|
|
$b = $point->getEllipsoid()->getB(); |
|
97
|
|
|
$f = 1 / $point->getEllipsoid()->getF(); |
|
98
|
|
|
|
|
99
|
|
|
$sinα1 = sin($α1); |
|
100
|
|
|
$cosα1 = cos($α1); |
|
101
|
|
|
|
|
102
|
|
|
$tanU1 = (1 - $f) * tan($φ1); |
|
103
|
|
|
$cosU1 = 1 / sqrt(1 + $tanU1 * $tanU1); |
|
104
|
|
|
$sinU1 = $tanU1 * $cosU1; |
|
105
|
|
|
$σ1 = atan2($tanU1, $cosα1); |
|
106
|
|
|
$sinα = $cosU1 * $sinα1; |
|
107
|
|
|
$cosSqα = 1 - $sinα * $sinα; |
|
108
|
|
|
$uSq = $cosSqα * ($a * $a - $b * $b) / ($b * $b); |
|
109
|
|
|
$A = 1 + $uSq / 16384 * (4096 + $uSq * (- 768 + $uSq * (320 - 175 * $uSq))); |
|
110
|
|
|
$B = $uSq / 1024 * (256 + $uSq * (- 128 + $uSq * (74 - 47 * $uSq))); |
|
111
|
|
|
|
|
112
|
|
|
$σ = $distance / ($b * $A); |
|
113
|
|
|
$iterations = 0; |
|
114
|
|
|
|
|
115
|
|
|
do { |
|
116
|
|
|
$cos2σm = cos(2 * $σ1 + $σ); |
|
117
|
|
|
$sinσ = sin($σ); |
|
118
|
|
|
$cosσ = cos($σ); |
|
119
|
|
|
$Δσ = $B * $sinσ * ($cos2σm + $B / 4 * ($cosσ * (- 1 + 2 * $cos2σm * $cos2σm) - $B / 6 * $cos2σm * (- 3 + 4 * $sinσ * $sinσ) * (- 3 + 4 * $cos2σm * $cos2σm))); |
|
120
|
|
|
$σs = $σ; |
|
121
|
|
|
$σ = $distance / ($b * $A) + $Δσ; |
|
122
|
|
|
} while (abs($σ - $σs) > 1e-12 && ++ $iterations < 200); |
|
123
|
|
|
|
|
124
|
|
|
if ($iterations >= 200) { |
|
125
|
|
|
throw new NotConvergingException('Inverse Vincenty Formula did not converge'); |
|
126
|
|
|
} |
|
127
|
|
|
|
|
128
|
|
|
$tmp = $sinU1 * $sinσ - $cosU1 * $cosσ * $cosα1; |
|
129
|
|
|
$φ2 = atan2($sinU1 * $cosσ + $cosU1 * $sinσ * $cosα1, (1 - $f) * sqrt($sinα * $sinα + $tmp * $tmp)); |
|
130
|
|
|
$λ = atan2($sinσ * $sinα1, $cosU1 * $cosσ - $sinU1 * $sinσ * $cosα1); |
|
131
|
|
|
$C = $f / 16 * $cosSqα * (4 + $f * (4 - 3 * $cosSqα)); |
|
132
|
|
|
$L = $λ - (1 - $C) * $f * $sinα * ($σ + $C * $sinσ * ($cos2σm + $C * $cosσ * (- 1 + 2 * $cos2σm * $cos2σm))); |
|
133
|
|
|
$λ2 = fmod($λ1 + $L + 3 * M_PI, 2 * M_PI) - M_PI; |
|
134
|
|
|
|
|
135
|
|
|
$α2 = atan2($sinα, - $tmp); |
|
136
|
|
|
$α2 = fmod($α2 + 2 * M_PI, 2 * M_PI); |
|
137
|
|
|
|
|
138
|
|
|
return [ |
|
139
|
|
|
'destination' => new Coordinate(rad2deg($φ2), rad2deg($λ2), $point->getEllipsoid()), |
|
140
|
|
|
'bearing_final' => rad2deg($α2), |
|
141
|
|
|
]; |
|
142
|
|
|
} |
|
143
|
|
|
|
|
144
|
|
|
private function inverseVincenty(Coordinate $point1, Coordinate $point2) |
|
145
|
|
|
{ |
|
146
|
|
|
$φ1 = deg2rad($point1->getLat()); |
|
147
|
|
|
$φ2 = deg2rad($point2->getLat()); |
|
148
|
|
|
$λ1 = deg2rad($point1->getLng()); |
|
149
|
|
|
$λ2 = deg2rad($point2->getLng()); |
|
150
|
|
|
|
|
151
|
|
|
$a = $point1->getEllipsoid()->getA(); |
|
152
|
|
|
$b = $point1->getEllipsoid()->getB(); |
|
153
|
|
|
$f = 1 / $point1->getEllipsoid()->getF(); |
|
154
|
|
|
|
|
155
|
|
|
$L = $λ2 - $λ1; |
|
156
|
|
|
|
|
157
|
|
|
$tanU1 = (1 - $f) * tan($φ1); |
|
158
|
|
|
$cosU1 = 1 / sqrt(1 + $tanU1 * $tanU1); |
|
159
|
|
|
$sinU1 = $tanU1 * $cosU1; |
|
160
|
|
|
$tanU2 = (1 - $f) * tan($φ2); |
|
161
|
|
|
$cosU2 = 1 / sqrt(1 + $tanU2 * $tanU2); |
|
162
|
|
|
$sinU2 = $tanU2 * $cosU2; |
|
163
|
|
|
|
|
164
|
|
|
$λ = $L; |
|
165
|
|
|
|
|
166
|
|
|
$iterations = 0; |
|
167
|
|
|
|
|
168
|
|
View Code Duplication |
do { |
|
|
|
|
|
|
169
|
|
|
$sinλ = sin($λ); |
|
170
|
|
|
$cosλ = cos($λ); |
|
171
|
|
|
$sinSqσ = ($cosU2 * $sinλ) * ($cosU2 * $sinλ) |
|
172
|
|
|
+ ($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosλ) * ($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosλ); |
|
173
|
|
|
$sinσ = sqrt($sinSqσ); |
|
174
|
|
|
|
|
175
|
|
|
if ($sinσ == 0) { |
|
176
|
|
|
return 0; |
|
177
|
|
|
} |
|
178
|
|
|
|
|
179
|
|
|
$cosσ = $sinU1 * $sinU2 + $cosU1 * $cosU2 * $cosλ; |
|
180
|
|
|
$σ = atan2($sinσ, $cosσ); |
|
181
|
|
|
$sinα = $cosU1 * $cosU2 * $sinλ / $sinσ; |
|
182
|
|
|
$cosSqα = 1 - $sinα * $sinα; |
|
183
|
|
|
|
|
184
|
|
|
if ($cosSqα == 0.0) { |
|
185
|
|
|
$cos2σM = 0; |
|
186
|
|
|
} else { |
|
187
|
|
|
$cos2σM = $cosσ - 2 * $sinU1 * $sinU2 / $cosSqα; |
|
188
|
|
|
} |
|
189
|
|
|
|
|
190
|
|
|
$C = $f / 16 * $cosSqα * (4 + $f * (4 - 3 * $cosSqα)); |
|
191
|
|
|
$λp = $λ; |
|
192
|
|
|
$λ = $L + (1 - $C) * $f * $sinα * ($σ + $C * $sinσ * ($cos2σM + $C * $cosσ * (- 1 + 2 * $cos2σM * $cos2σM))); |
|
193
|
|
|
} while (abs($λ - $λp) > 1e-12 && ++ $iterations < 200); |
|
194
|
|
|
|
|
195
|
|
|
if ($iterations >= 200) { |
|
196
|
|
|
throw new NotConvergingException('Inverse Vincenty Formula did not converge'); |
|
197
|
|
|
} |
|
198
|
|
|
|
|
199
|
|
|
$uSq = $cosSqα * ($a * $a - $b * $b) / ($b * $b); |
|
200
|
|
|
$A = 1 + $uSq / 16384 * (4096 + $uSq * (- 768 + $uSq * (320 - 175 * $uSq))); |
|
201
|
|
|
$B = $uSq / 1024 * (256 + $uSq * (- 128 + $uSq * (74 - 47 * $uSq))); |
|
202
|
|
|
$Δσ = $B * $sinσ * ($cos2σM + $B / 4 * ($cosσ * (- 1 + 2 * $cos2σM * $cos2σM) - $B / 6 * $cos2σM * (- 3 + 4 * $sinσ * $sinσ) * (- 3 + 4 * $cos2σM * $cos2σM))); |
|
203
|
|
|
|
|
204
|
|
|
$s = $b * $A * ($σ - $Δσ); |
|
205
|
|
|
|
|
206
|
|
|
$α1 = atan2($cosU2 * $sinλ, $cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosλ); |
|
207
|
|
|
$α2 = atan2($cosU1 * $sinλ, - $sinU1 * $cosU2 + $cosU1 * $sinU2 * $cosλ); |
|
208
|
|
|
|
|
209
|
|
|
$α1 = fmod($α1 + 2 * M_PI, 2 * M_PI); |
|
210
|
|
|
$α2 = fmod($α2 + 2 * M_PI, 2 * M_PI); |
|
211
|
|
|
|
|
212
|
|
|
$s = round($s, 3); |
|
213
|
|
|
|
|
214
|
|
|
return [ |
|
215
|
|
|
'distance' => $s, |
|
216
|
|
|
'bearing_initial' => rad2deg($α1), |
|
217
|
|
|
'bearing_final' => rad2deg($α2), |
|
218
|
|
|
]; |
|
219
|
|
|
} |
|
220
|
|
|
} |
|
221
|
|
|
|
mjaschen /
phpgeo
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation.
You can also find more detailed suggestions in the “Code” section of your repository.