BearingEllipsoidal::inverseVincenty() - Code Metrics - mjaschen/phpgeo - Measure and Improve Code Quality continuously with Scrutinizer

BearingEllipsoidal::inverseVincenty() B
last analyzed 2021-03-04 10:04 UTC

↳ Parent: BearingEllipsoidal

Complexity

Conditions	6
Paths	8

Size

Total Lines

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
dl	0
loc	78
rs	7.8577
c	0
b	0
f	0
cc	6
nc	8
nop	2

How to fix Long Method

<?php

declare(strict_types=1);

namespace Location\Bearing;

use InvalidArgumentException;
use Location\Coordinate;
use Location\Exception\NotConvergingException;

/**
 * Calculation of bearing between two points using a
 * ellipsoidal model of the earth.
 *
 * This class is based on the awesome work Chris Veness
 * has done. For more information visit the following URL.
 *
 * @see http://www.movable-type.co.uk/scripts/latlong-vincenty.html
 *
 * @author Marcus Jaschen <[email protected]>
 */
class BearingEllipsoidal implements BearingInterface
{
    /**
     * This method calculates the initial bearing between the
     * two points.
     *
     * If the two points share the same location, the bearing
     * value will be 0.0.
     *
     * @param Coordinate $point1
     * @param Coordinate $point2
     *
     * @return float Bearing Angle
     */
    public function calculateBearing(Coordinate $point1, Coordinate $point2): float
    {
        if ($point1->hasSameLocation($point2)) {
            return 0.0;
        }

        return $this->inverseVincenty($point1, $point2)->getBearingInitial();
    }

    /**
     * Calculates the final bearing between the two points.
     *
     * @param Coordinate $point1
     * @param Coordinate $point2
     *
     * @return float
     */
    public function calculateFinalBearing(Coordinate $point1, Coordinate $point2): float
    {
        return $this->inverseVincenty($point1, $point2)->getBearingFinal();
    }

    /**
     * Calculates a destination point for the given point, bearing angle,
     * and distance.
     *
     * @param Coordinate $point
     * @param float $bearing the bearing angle between 0 and 360 degrees
     * @param float $distance the distance to the destination point in meters
     *
     * @return Coordinate
     */
    public function calculateDestination(Coordinate $point, float $bearing, float $distance): Coordinate
    {
        return $this->directVincenty($point, $bearing, $distance)->getDestination();
    }

    /**
     * Calculates the final bearing angle for a destination point.
     * The method expects a starting point point, the bearing angle,
     * and the distance to destination.
     *
     * @param Coordinate $point
     * @param float $bearing
     * @param float $distance
     *
     * @return float
     *
     * @throws NotConvergingException
     */
    public function calculateDestinationFinalBearing(Coordinate $point, float $bearing, float $distance): float
    {
        return $this->directVincenty($point, $bearing, $distance)->getBearingFinal();
    }

    /**
     * @param Coordinate $point
     * @param float $bearing
     * @param float $distance
     *
     * @return DirectVincentyBearing
     *
     * @throws NotConvergingException
     */
    private function directVincenty(Coordinate $point, float $bearing, float $distance): DirectVincentyBearing
    {
        $phi1 = deg2rad($point->getLat());
        $lambda1 = deg2rad($point->getLng());
        $alpha1 = deg2rad($bearing);

        $a = $point->getEllipsoid()->getA();
        $b = $point->getEllipsoid()->getB();
        $f = 1 / $point->getEllipsoid()->getF();

        $sinAlpha1 = sin($alpha1);
        $cosAlpha1 = cos($alpha1);

        $tanU1 = (1 - $f) * tan($phi1);
        $cosU1 = 1 / sqrt(1 + $tanU1 * $tanU1);
        $sinU1 = $tanU1 * $cosU1;
        $sigma1 = atan2($tanU1, $cosAlpha1);
        $sinAlpha = $cosU1 * $sinAlpha1;
        $cosSquAlpha = 1 - $sinAlpha * $sinAlpha;
        $uSq = $cosSquAlpha * ($a * $a - $b * $b) / ($b * $b);
        $A = 1 + $uSq / 16384 * (4096 + $uSq * (-768 + $uSq * (320 - 175 * $uSq)));
        $B = $uSq / 1024 * (256 + $uSq * (-128 + $uSq * (74 - 47 * $uSq)));

        $sigmaS = $distance / ($b * $A);
        $sigma = $sigmaS;
        $iterations = 0;

        do {
            $cos2SigmaM = cos(2 * $sigma1 + $sigma);
            $sinSigma = sin($sigma);
            $cosSigma = cos($sigma);
            $deltaSigma = $B * $sinSigma
                * ($cos2SigmaM + $B / 4
                    * ($cosSigma
                        * (-1 + 2 * $cos2SigmaM * $cos2SigmaM) - $B / 6
                        * $cos2SigmaM * (-3 + 4 * $sinSigma * $sinSigma)
                        * (-3 + 4 * $cos2SigmaM * $cos2SigmaM)
                    )
                );
            $sigmaS = $sigma;
            $sigma = $distance / ($b * $A) + $deltaSigma;
        } while (abs($sigma - $sigmaS) > 1e-12 && ++$iterations < 200);

        if ($iterations >= 200) {
            throw new NotConvergingException('Inverse Vincenty Formula did not converge');
        }

        $tmp = $sinU1 * $sinSigma - $cosU1 * $cosSigma * $cosAlpha1;
        $phi2 = atan2(
            $sinU1 * $cosSigma + $cosU1 * $sinSigma * $cosAlpha1,
            (1 - $f) * sqrt($sinAlpha * $sinAlpha + $tmp * $tmp)
        );
        $lambda = atan2($sinSigma * $sinAlpha1, $cosU1 * $cosSigma - $sinU1 * $sinSigma * $cosAlpha1);
        $C = $f / 16 * $cosSquAlpha * (4 + $f * (4 - 3 * $cosSquAlpha));
        $L = $lambda
            - (1 - $C) * $f * $sinAlpha
            * ($sigma + $C * $sinSigma * ($cos2SigmaM + $C * $cosSigma * (-1 + 2 * $cos2SigmaM ** 2)));
        $lambda2 = fmod($lambda1 + $L + 3 * M_PI, 2 * M_PI) - M_PI;

        $alpha2 = atan2($sinAlpha, -$tmp);
        $alpha2 = fmod($alpha2 + 2 * M_PI, 2 * M_PI);

        return new DirectVincentyBearing(
            new Coordinate(rad2deg($phi2), rad2deg($lambda2), $point->getEllipsoid()),
            rad2deg($alpha2)
        );
    }

    /**
     * @param Coordinate $point1
     * @param Coordinate $point2
     *
     * @return InverseVincentyBearing
     *
     * @throws NotConvergingException
     */
    private function inverseVincenty(Coordinate $point1, Coordinate $point2): InverseVincentyBearing
    {
        $φ1 = deg2rad($point1->getLat());
        $φ2 = deg2rad($point2->getLat());
        $λ1 = deg2rad($point1->getLng());
        $λ2 = deg2rad($point2->getLng());

        $a = $point1->getEllipsoid()->getA();
        $b = $point1->getEllipsoid()->getB();
        $f = 1 / $point1->getEllipsoid()->getF();

        $L = $λ2 - $λ1;

        $tanU1 = (1 - $f) * tan($φ1);
        $cosU1 = 1 / sqrt(1 + $tanU1 * $tanU1);
        $sinU1 = $tanU1 * $cosU1;
        $tanU2 = (1 - $f) * tan($φ2);
        $cosU2 = 1 / sqrt(1 + $tanU2 * $tanU2);
        $sinU2 = $tanU2 * $cosU2;

        $λ = $L;

        $iterations = 0;

        do {
            $sinλ = sin($λ);
            $cosλ = cos($λ);
            $sinSqσ = ($cosU2 * $sinλ) * ($cosU2 * $sinλ)
                + ($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosλ) * ($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosλ);
            $sinσ = sqrt($sinSqσ);

            if ($sinσ == 0) {
                new InverseVincentyBearing(0, 0, 0);
            }

            $cosσ = $sinU1 * $sinU2 + $cosU1 * $cosU2 * $cosλ;
            $σ = atan2($sinσ, $cosσ);
            $sinα = $cosU1 * $cosU2 * $sinλ / $sinσ;
            $cosSqα = 1 - $sinα * $sinα;

            $cos2σM = 0;
            if ($cosSqα !== 0.0) {
                $cos2σM = $cosσ - 2 * $sinU1 * $sinU2 / $cosSqα;
            }

            $C = $f / 16 * $cosSqα * (4 + $f * (4 - 3 * $cosSqα));
            $λp = $λ;
            $λ = $L + (1 - $C) * $f * $sinα
                * ($σ + $C * $sinσ * ($cos2σM + $C * $cosσ * (-1 + 2 * $cos2σM * $cos2σM)));
        } while (abs($λ - $λp) > 1e-12 && ++$iterations < 200);

        if ($iterations >= 200) {
            throw new NotConvergingException('Inverse Vincenty Formula did not converge');
        }

        $uSq = $cosSqα * ($a * $a - $b * $b) / ($b * $b);
        $A = 1 + $uSq / 16384 * (4096 + $uSq * (-768 + $uSq * (320 - 175 * $uSq)));
        $B = $uSq / 1024 * (256 + $uSq * (-128 + $uSq * (74 - 47 * $uSq)));
        $Δσ = $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)
                )
            );

        $s = $b * $A * ($σ - $Δσ);

        $α1 = atan2($cosU2 * $sinλ, $cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosλ);
        $α2 = atan2($cosU1 * $sinλ, -$sinU1 * $cosU2 + $cosU1 * $sinU2 * $cosλ);

        $α1 = fmod($α1 + 2 * M_PI, 2 * M_PI);
        $α2 = fmod($α2 + 2 * M_PI, 2 * M_PI);

        $s = round($s, 3);

        return new InverseVincentyBearing($s, rad2deg($α1), rad2deg($α2));
    }
}


1			<?php
2
3			declare(strict_types=1);
4
5			namespace Location\Bearing;
6
7			use InvalidArgumentException;
8			use Location\Coordinate;
9			use Location\Exception\NotConvergingException;
10
11			/**
12			* Calculation of bearing between two points using a
13			* ellipsoidal model of the earth.
14			*
15			* This class is based on the awesome work Chris Veness
16			* has done. For more information visit the following URL.
17			*
18			* @see http://www.movable-type.co.uk/scripts/latlong-vincenty.html
19			*
20			* @author Marcus Jaschen <[email protected]>
21			*/
22			class BearingEllipsoidal implements BearingInterface
23			{
24			/**
25			* This method calculates the initial bearing between the
26			* two points.
27			*
28			* If the two points share the same location, the bearing
29			* value will be 0.0.
30			*
31			* @param Coordinate $point1
32			* @param Coordinate $point2
33			*
34			* @return float Bearing Angle
35			*/
36			public function calculateBearing(Coordinate $point1, Coordinate $point2): float
37			{
38			if ($point1->hasSameLocation($point2)) {
39			return 0.0;
40			}
41
42			return $this->inverseVincenty($point1, $point2)->getBearingInitial();
43			}
44
45			/**
46			* Calculates the final bearing between the two points.
47			*
48			* @param Coordinate $point1
49			* @param Coordinate $point2
50			*
51			* @return float
52			*/
53			public function calculateFinalBearing(Coordinate $point1, Coordinate $point2): float
54			{
55			return $this->inverseVincenty($point1, $point2)->getBearingFinal();
56			}
57
58			/**
59			* Calculates a destination point for the given point, bearing angle,
60			* and distance.
61			*
62			* @param 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, float $bearing, float $distance): Coordinate
69			{
70			return $this->directVincenty($point, $bearing, $distance)->getDestination();
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 Coordinate $point
79			* @param float $bearing
80			* @param float $distance
81			*
82			* @return float
83			*
84			* @throws NotConvergingException
85			*/
86			public function calculateDestinationFinalBearing(Coordinate $point, float $bearing, float $distance): float
87			{
88			return $this->directVincenty($point, $bearing, $distance)->getBearingFinal();
89			}
90
91			/**
92			* @param Coordinate $point
93			* @param float $bearing
94			* @param float $distance
95			*
96			* @return DirectVincentyBearing
97			*
98			* @throws NotConvergingException
99			*/
100			private function directVincenty(Coordinate $point, float $bearing, float $distance): DirectVincentyBearing
101			{
102			$phi1 = deg2rad($point->getLat());
103			$lambda1 = deg2rad($point->getLng());
104			$alpha1 = deg2rad($bearing);
105
106			$a = $point->getEllipsoid()->getA();
107			$b = $point->getEllipsoid()->getB();
108			$f = 1 / $point->getEllipsoid()->getF();
109
110			$sinAlpha1 = sin($alpha1);
111			$cosAlpha1 = cos($alpha1);
112
113			$tanU1 = (1 - $f) * tan($phi1);
114			$cosU1 = 1 / sqrt(1 + $tanU1 * $tanU1);
115			$sinU1 = $tanU1 * $cosU1;
116			$sigma1 = atan2($tanU1, $cosAlpha1);
117			$sinAlpha = $cosU1 * $sinAlpha1;
118			$cosSquAlpha = 1 - $sinAlpha * $sinAlpha;
119			$uSq = $cosSquAlpha * ($a * $a - $b * $b) / ($b * $b);
120			$A = 1 + $uSq / 16384 * (4096 + $uSq * (-768 + $uSq * (320 - 175 * $uSq)));
121			$B = $uSq / 1024 * (256 + $uSq * (-128 + $uSq * (74 - 47 * $uSq)));
122
123			$sigmaS = $distance / ($b * $A);
124			$sigma = $sigmaS;
125			$iterations = 0;
126
127			do {
128			$cos2SigmaM = cos(2 * $sigma1 + $sigma);
129			$sinSigma = sin($sigma);
130			$cosSigma = cos($sigma);
131			$deltaSigma = $B * $sinSigma
132			* ($cos2SigmaM + $B / 4
133			* ($cosSigma
134			* (-1 + 2 * $cos2SigmaM * $cos2SigmaM) - $B / 6
135			* $cos2SigmaM * (-3 + 4 * $sinSigma * $sinSigma)
136			* (-3 + 4 * $cos2SigmaM * $cos2SigmaM)
137			)
138			);
139			$sigmaS = $sigma;
140			$sigma = $distance / ($b * $A) + $deltaSigma;
141			} while (abs($sigma - $sigmaS) > 1e-12 && ++$iterations < 200);
142
143			if ($iterations >= 200) {
144			throw new NotConvergingException('Inverse Vincenty Formula did not converge');
145			}
146
147			$tmp = $sinU1 * $sinSigma - $cosU1 * $cosSigma * $cosAlpha1;
148			$phi2 = atan2(
149			$sinU1 * $cosSigma + $cosU1 * $sinSigma * $cosAlpha1,
150			(1 - $f) * sqrt($sinAlpha * $sinAlpha + $tmp * $tmp)
151			);
152			$lambda = atan2($sinSigma * $sinAlpha1, $cosU1 * $cosSigma - $sinU1 * $sinSigma * $cosAlpha1);
153			$C = $f / 16 * $cosSquAlpha * (4 + $f * (4 - 3 * $cosSquAlpha));
154			$L = $lambda
155			- (1 - $C) * $f * $sinAlpha
156			* ($sigma + $C * $sinSigma * ($cos2SigmaM + $C * $cosSigma * (-1 + 2 * $cos2SigmaM ** 2)));
157			$lambda2 = fmod($lambda1 + $L + 3 * M_PI, 2 * M_PI) - M_PI;
158
159			$alpha2 = atan2($sinAlpha, -$tmp);
160			$alpha2 = fmod($alpha2 + 2 * M_PI, 2 * M_PI);
161
162			return new DirectVincentyBearing(
163			new Coordinate(rad2deg($phi2), rad2deg($lambda2), $point->getEllipsoid()),
164			rad2deg($alpha2)
165			);
166			}
167
168			/**
169			* @param Coordinate $point1
170			* @param Coordinate $point2
171			*
172			* @return InverseVincentyBearing
173			*
174			* @throws NotConvergingException
175			*/
176			private function inverseVincenty(Coordinate $point1, Coordinate $point2): InverseVincentyBearing
177			{
178			$φ1 = deg2rad($point1->getLat());
179			$φ2 = deg2rad($point2->getLat());
180			$λ1 = deg2rad($point1->getLng());
181			$λ2 = deg2rad($point2->getLng());
182
183			$a = $point1->getEllipsoid()->getA();
184			$b = $point1->getEllipsoid()->getB();
185			$f = 1 / $point1->getEllipsoid()->getF();
186
187			$L = $λ2 - $λ1;
188
189			$tanU1 = (1 - $f) * tan($φ1);
190			$cosU1 = 1 / sqrt(1 + $tanU1 * $tanU1);
191			$sinU1 = $tanU1 * $cosU1;
192			$tanU2 = (1 - $f) * tan($φ2);
193			$cosU2 = 1 / sqrt(1 + $tanU2 * $tanU2);
194			$sinU2 = $tanU2 * $cosU2;
195
196			$λ = $L;
197
198			$iterations = 0;
199
200			do {
201			$sinλ = sin($λ);
202			$cosλ = cos($λ);
203			$sinSqσ = ($cosU2 * $sinλ) * ($cosU2 * $sinλ)
204			+ ($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosλ) * ($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosλ);
205			$sinσ = sqrt($sinSqσ);
206
207			if ($sinσ == 0) {
208			new InverseVincentyBearing(0, 0, 0);
209			}
210
211			$cosσ = $sinU1 * $sinU2 + $cosU1 * $cosU2 * $cosλ;
212			$σ = atan2($sinσ, $cosσ);
213			$sinα = $cosU1 * $cosU2 * $sinλ / $sinσ;
214			$cosSqα = 1 - $sinα * $sinα;
215
216			$cos2σM = 0;
217			if ($cosSqα !== 0.0) {
218			$cos2σM = $cosσ - 2 * $sinU1 * $sinU2 / $cosSqα;
219			}
220
221			$C = $f / 16 * $cosSqα * (4 + $f * (4 - 3 * $cosSqα));
222			$λp = $λ;
223			$λ = $L + (1 - $C) * $f * $sinα
224			* ($σ + $C * $sinσ * ($cos2σM + $C * $cosσ * (-1 + 2 * $cos2σM * $cos2σM)));
225			} while (abs($λ - $λp) > 1e-12 && ++$iterations < 200);
226
227			if ($iterations >= 200) {
228			throw new NotConvergingException('Inverse Vincenty Formula did not converge');
229			}
230
231			$uSq = $cosSqα * ($a * $a - $b * $b) / ($b * $b);
232			$A = 1 + $uSq / 16384 * (4096 + $uSq * (-768 + $uSq * (320 - 175 * $uSq)));
233			$B = $uSq / 1024 * (256 + $uSq * (-128 + $uSq * (74 - 47 * $uSq)));
234			$Δσ = $B * $sinσ
235			* ($cos2σM + $B / 4
236			* ($cosσ * (-1 + 2 * $cos2σM * $cos2σM) - $B / 6
237			* $cos2σM * (-3 + 4 * $sinσ * $sinσ)
238			* (-3 + 4 * $cos2σM * $cos2σM)
239			)
240			);
241
242			$s = $b * $A * ($σ - $Δσ);
243
244			$α1 = atan2($cosU2 * $sinλ, $cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosλ);
245			$α2 = atan2($cosU1 * $sinλ, -$sinU1 * $cosU2 + $cosU1 * $sinU2 * $cosλ);
246
247			$α1 = fmod($α1 + 2 * M_PI, 2 * M_PI);
248			$α2 = fmod($α2 + 2 * M_PI, 2 * M_PI);
249
250			$s = round($s, 3);
251
252			return new InverseVincentyBearing($s, rad2deg($α1), rad2deg($α2));
253			}
254			}
255

mjaschen / phpgeo

BearingEllipsoidal::inverseVincenty() B last analyzed 2021-03-04 10:04 UTC

Complexity

Size

Duplication

Importance

How to fix Long Method

Long Method

Duplication Side-by-Side

Filter issues like

BearingEllipsoidal::inverseVincenty() B
last analyzed 2021-03-04 10:04 UTC