BearingEllipsoidal::inverseVincenty() - Code Metrics - Inspection of "update for new milestone" - mjaschen/phpgeo - Measure and Improve Code Quality continuously with Scrutinizer

Completed

Push — milestone/4.0 ( 4a7219...1803cd )

by Marcus

created 2020-02-08 09:21 UTC

BearingEllipsoidal::inverseVincenty() B

↳ Parent: BearingEllipsoidal

Complexity

Conditions	6
Paths	8

Size

Total Lines

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
dl	0
loc	77
rs	7.8795
c	0
b	0
f	0
cc	6
nc	8
nop	2

How to fix Long Method

<?php

declare(strict_types=1);

namespace Phpgeo\Bearing;

use InvalidArgumentException;
use Phpgeo\Point;
use Phpgeo\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.
     *
     * @param Point $point1
     * @param Point $point2
     *
     * @return float Bearing Angle
     */
    public function calculateBearing(Point $point1, Point $point2): float
    {
        return $this->inverseVincenty($point1, $point2)->getBearingInitial();
    }

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

    /**
     * Calculates a destination point for the given point, bearing angle,
     * and distance.
     *
     * @param Point $point
     * @param float $bearing the bearing angle between 0 and 360 degrees
     * @param float $distance the distance to the destination point in meters
     *
     * @return Point
     */
    public function calculateDestination(Point $point, float $bearing, float $distance): Point
    {
        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 Point $point
     * @param float $bearing
     * @param float $distance
     *
     * @return float
     *
     * @throws NotConvergingException
     */
    public function calculateDestinationFinalBearing(Point $point, float $bearing, float $distance): float
    {
        return $this->directVincenty($point, $bearing, $distance)->getBearingFinal();
    }

    /**
     * @param Point $point
     * @param float $bearing
     * @param float $distance
     *
     * @return DirectVincentyBearing
     *
     * @throws NotConvergingException
     */
    private function directVincenty(Point $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 Point(rad2deg($phi2), rad2deg($lambda2), $point->getEllipsoid()),
            rad2deg($alpha2)
        );
    }

    /**
     * @param Point $point1
     * @param Point $point2
     *
     * @return InverseVincentyBearing
     * @throws NotConvergingException
     */
    private function inverseVincenty(Point $point1, Point $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 Phpgeo\Bearing;
6
7			use InvalidArgumentException;
8			use Phpgeo\Point;
9			use Phpgeo\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			* @param Point $point1
29			* @param Point $point2
30			*
31			* @return float Bearing Angle
32			*/
33			public function calculateBearing(Point $point1, Point $point2): float
34			{
35			return $this->inverseVincenty($point1, $point2)->getBearingInitial();
36			}
37
38			/**
39			* Calculates the final bearing between the two points.
40			*
41			* @param Point $point1
42			* @param Point $point2
43			*
44			* @return float
45			*/
46			public function calculateFinalBearing(Point $point1, Point $point2): float
47			{
48			return $this->inverseVincenty($point1, $point2)->getBearingFinal();
49			}
50
51			/**
52			* Calculates a destination point for the given point, bearing angle,
53			* and distance.
54			*
55			* @param Point $point
56			* @param float $bearing the bearing angle between 0 and 360 degrees
57			* @param float $distance the distance to the destination point in meters
58			*
59			* @return Point
60			*/
61			public function calculateDestination(Point $point, float $bearing, float $distance): Point
62			{
63			return $this->directVincenty($point, $bearing, $distance)->getDestination();
64			}
65
66			/**
67			* Calculates the final bearing angle for a destination point.
68			* The method expects a starting point point, the bearing angle,
69			* and the distance to destination.
70			*
71			* @param Point $point
72			* @param float $bearing
73			* @param float $distance
74			*
75			* @return float
76			*
77			* @throws NotConvergingException
78			*/
79			public function calculateDestinationFinalBearing(Point $point, float $bearing, float $distance): float
80			{
81			return $this->directVincenty($point, $bearing, $distance)->getBearingFinal();
82			}
83
84			/**
85			* @param Point $point
86			* @param float $bearing
87			* @param float $distance
88			*
89			* @return DirectVincentyBearing
90			*
91			* @throws NotConvergingException
92			*/
93			private function directVincenty(Point $point, float $bearing, float $distance): DirectVincentyBearing
94			{
95			$phi1 = deg2rad($point->getLat());
96			$lambda1 = deg2rad($point->getLng());
97			$alpha1 = deg2rad($bearing);
98
99			$a = $point->getEllipsoid()->getA();
100			$b = $point->getEllipsoid()->getB();
101			$f = 1 / $point->getEllipsoid()->getF();
102
103			$sinAlpha1 = sin($alpha1);
104			$cosAlpha1 = cos($alpha1);
105
106			$tanU1 = (1 - $f) * tan($phi1);
107			$cosU1 = 1 / sqrt(1 + $tanU1 * $tanU1);
108			$sinU1 = $tanU1 * $cosU1;
109			$sigma1 = atan2($tanU1, $cosAlpha1);
110			$sinAlpha = $cosU1 * $sinAlpha1;
111			$cosSquAlpha = 1 - $sinAlpha * $sinAlpha;
112			$uSq = $cosSquAlpha * ($a * $a - $b * $b) / ($b * $b);
113			$A = 1 + $uSq / 16384 * (4096 + $uSq * (-768 + $uSq * (320 - 175 * $uSq)));
114			$B = $uSq / 1024 * (256 + $uSq * (-128 + $uSq * (74 - 47 * $uSq)));
115
116			$sigmaS = $distance / ($b * $A);
117			$sigma = $sigmaS;
118			$iterations = 0;
119
120			do {
121			$cos2SigmaM = cos(2 * $sigma1 + $sigma);
122			$sinSigma = sin($sigma);
123			$cosSigma = cos($sigma);
124			$deltaSigma = $B * $sinSigma
125			* ($cos2SigmaM + $B / 4
126			* ($cosSigma
127			* (-1 + 2 * $cos2SigmaM * $cos2SigmaM) - $B / 6
128			* $cos2SigmaM * (-3 + 4 * $sinSigma * $sinSigma)
129			* (-3 + 4 * $cos2SigmaM * $cos2SigmaM)
130			)
131			);
132			$sigmaS = $sigma;
133			$sigma = $distance / ($b * $A) + $deltaSigma;
134			} while (abs($sigma - $sigmaS) > 1e-12 && ++$iterations < 200);
135
136			if ($iterations >= 200) {
137			throw new NotConvergingException('Inverse Vincenty Formula did not converge');
138			}
139
140			$tmp = $sinU1 * $sinSigma - $cosU1 * $cosSigma * $cosAlpha1;
141			$phi2 = atan2(
142			$sinU1 * $cosSigma + $cosU1 * $sinSigma * $cosAlpha1,
143			(1 - $f) * sqrt($sinAlpha * $sinAlpha + $tmp * $tmp)
144			);
145			$lambda = atan2($sinSigma * $sinAlpha1, $cosU1 * $cosSigma - $sinU1 * $sinSigma * $cosAlpha1);
146			$C = $f / 16 * $cosSquAlpha * (4 + $f * (4 - 3 * $cosSquAlpha));
147			$L = $lambda
148			- (1 - $C) * $f * $sinAlpha
149			* ($sigma + $C * $sinSigma * ($cos2SigmaM + $C * $cosSigma * (-1 + 2 * $cos2SigmaM ** 2)));
150			$lambda2 = fmod($lambda1 + $L + 3 * M_PI, 2 * M_PI) - M_PI;
151
152			$alpha2 = atan2($sinAlpha, -$tmp);
153			$alpha2 = fmod($alpha2 + 2 * M_PI, 2 * M_PI);
154
155			return new DirectVincentyBearing(
156			new Point(rad2deg($phi2), rad2deg($lambda2), $point->getEllipsoid()),
157			rad2deg($alpha2)
158			);
159			}
160
161			/**
162			* @param Point $point1
163			* @param Point $point2
164			*
165			* @return InverseVincentyBearing
166			* @throws NotConvergingException
167			*/
168			private function inverseVincenty(Point $point1, Point $point2): InverseVincentyBearing
169			{
170			$φ1 = deg2rad($point1->getLat());
171			$φ2 = deg2rad($point2->getLat());
172			$λ1 = deg2rad($point1->getLng());
173			$λ2 = deg2rad($point2->getLng());
174
175			$a = $point1->getEllipsoid()->getA();
176			$b = $point1->getEllipsoid()->getB();
177			$f = 1 / $point1->getEllipsoid()->getF();
178
179			$L = $λ2 - $λ1;
180
181			$tanU1 = (1 - $f) * tan($φ1);
182			$cosU1 = 1 / sqrt(1 + $tanU1 * $tanU1);
183			$sinU1 = $tanU1 * $cosU1;
184			$tanU2 = (1 - $f) * tan($φ2);
185			$cosU2 = 1 / sqrt(1 + $tanU2 * $tanU2);
186			$sinU2 = $tanU2 * $cosU2;
187
188			$λ = $L;
189
190			$iterations = 0;
191
192			do {
193			$sinλ = sin($λ);
194			$cosλ = cos($λ);
195			$sinSqσ = ($cosU2 * $sinλ) * ($cosU2 * $sinλ)
196			+ ($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosλ) * ($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosλ);
197			$sinσ = sqrt($sinSqσ);
198
199			if ($sinσ == 0) {
200			new InverseVincentyBearing(0, 0, 0);
201			}
202
203			$cosσ = $sinU1 * $sinU2 + $cosU1 * $cosU2 * $cosλ;
204			$σ = atan2($sinσ, $cosσ);
205			$sinα = $cosU1 * $cosU2 * $sinλ / $sinσ;
206			$cosSqα = 1 - $sinα * $sinα;
207
208			$cos2σM = 0;
209			if ($cosSqα !== 0.0) {
210			$cos2σM = $cosσ - 2 * $sinU1 * $sinU2 / $cosSqα;
211			}
212
213			$C = $f / 16 * $cosSqα * (4 + $f * (4 - 3 * $cosSqα));
214			$λp = $λ;
215			$λ = $L + (1 - $C) * $f * $sinα * ($σ + $C * $sinσ * ($cos2σM + $C * $cosσ * (-1 + 2 * $cos2σM * $cos2σM)));
216			} while (abs($λ - $λp) > 1e-12 && ++$iterations < 200);
217
218			if ($iterations >= 200) {
219			throw new NotConvergingException('Inverse Vincenty Formula did not converge');
220			}
221
222			$uSq = $cosSqα * ($a * $a - $b * $b) / ($b * $b);
223			$A = 1 + $uSq / 16384 * (4096 + $uSq * (-768 + $uSq * (320 - 175 * $uSq)));
224			$B = $uSq / 1024 * (256 + $uSq * (-128 + $uSq * (74 - 47 * $uSq)));
225			$Δσ = $B * $sinσ
226			* ($cos2σM + $B / 4
227			* ($cosσ * (-1 + 2 * $cos2σM * $cos2σM) - $B / 6
228			* $cos2σM * (-3 + 4 * $sinσ * $sinσ)
229			* (-3 + 4 * $cos2σM * $cos2σM)
230			)
231			);
232
233			$s = $b * $A * ($σ - $Δσ);
234
235			$α1 = atan2($cosU2 * $sinλ, $cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosλ);
236			$α2 = atan2($cosU1 * $sinλ, -$sinU1 * $cosU2 + $cosU1 * $sinU2 * $cosλ);
237
238			$α1 = fmod($α1 + 2 * M_PI, 2 * M_PI);
239			$α2 = fmod($α2 + 2 * M_PI, 2 * M_PI);
240
241			$s = round($s, 3);
242
243			return new InverseVincentyBearing($s, rad2deg($α1), rad2deg($α2));
244			}
245			}
246

mjaschen / phpgeo

Push — milestone/4.0 ( 4a7219...1803cd )

BearingEllipsoidal::inverseVincenty() B

Complexity

Size

Duplication

Importance

How to fix Long Method

Long Method

Duplication Side-by-Side

Filter issues like