|
| 1 | +<?php |
| 2 | + |
| 3 | +/* |
| 4 | + * This file is part of the Symfony package. |
| 5 | + * |
| 6 | + * (c) Fabien Potencier <[email protected]> |
| 7 | + * |
| 8 | + * For the full copyright and license information, please view the LICENSE |
| 9 | + * file that was distributed with this source code. |
| 10 | + */ |
| 11 | + |
| 12 | +namespace Symfony\UX\Map\Distance; |
| 13 | + |
| 14 | +use Symfony\UX\Map\Point; |
| 15 | + |
| 16 | +/** |
| 17 | + * Vincenty formula-based distance calculator. |
| 18 | + * |
| 19 | + * This calculator is more accurate than the Haversine formula, but slower. |
| 20 | + * |
| 21 | + * @author Simon André <[email protected]> |
| 22 | + */ |
| 23 | +final readonly class VincentyDistanceCalculator implements DistanceCalculatorInterface |
| 24 | +{ |
| 25 | + /** |
| 26 | + * WS-84 ellipsoid parameters. |
| 27 | + */ |
| 28 | + // Major Axis in meters |
| 29 | + private const A = 6378137.0; |
| 30 | + // Flattening |
| 31 | + private const F = 1 / 298.257223563; |
| 32 | + // Minor Axis in meters |
| 33 | + private const B = 6356752.314245; |
| 34 | + |
| 35 | + public function calculateDistance(Point $point1, Point $point2): float |
| 36 | + { |
| 37 | + $phi1 = deg2rad($point1->getLatitude()); |
| 38 | + $phi2 = deg2rad($point2->getLatitude()); |
| 39 | + $lambda1 = deg2rad($point1->getLongitude()); |
| 40 | + $lambda2 = deg2rad($point2->getLongitude()); |
| 41 | + |
| 42 | + $L = $lambda2 - $lambda1; |
| 43 | + $U1 = atan((1 - self::F) * tan($phi1)); |
| 44 | + $U2 = atan((1 - self::F) * tan($phi2)); |
| 45 | + $sinU1 = sin($U1); |
| 46 | + $cosU1 = cos($U1); |
| 47 | + $sinU2 = sin($U2); |
| 48 | + $cosU2 = cos($U2); |
| 49 | + |
| 50 | + $lambda = $L; |
| 51 | + $iterLimit = 100; |
| 52 | + do { |
| 53 | + $sinLambda = sin($lambda); |
| 54 | + $cosLambda = cos($lambda); |
| 55 | + $sinSigma = sqrt(($cosU2 * $sinLambda) ** 2 |
| 56 | + + ($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosLambda) ** 2); |
| 57 | + |
| 58 | + if (0.0 === $sinSigma) { |
| 59 | + return 0.0; |
| 60 | + } |
| 61 | + |
| 62 | + $cosSigma = $sinU1 * $sinU2 + $cosU1 * $cosU2 * $cosLambda; |
| 63 | + $sigma = atan2($sinSigma, $cosSigma); |
| 64 | + $sinAlpha = $cosU1 * $cosU2 * $sinLambda / $sinSigma; |
| 65 | + $cosSqAlpha = 1 - $sinAlpha * $sinAlpha; |
| 66 | + $cos2SigmaM = (0.0 === $cosSqAlpha) ? 0.0 : $cosSigma - 2 * $sinU1 * $sinU2 / $cosSqAlpha; |
| 67 | + $C = self::F / 16 * $cosSqAlpha * (4 + self::F * (4 - 3 * $cosSqAlpha)); |
| 68 | + |
| 69 | + $lambdaPrev = $lambda; |
| 70 | + $lambda = $L + (1 - $C) * self::F * $sinAlpha |
| 71 | + * ($sigma + $C * $sinSigma * ($cos2SigmaM + $C * $cosSigma * (-1 + 2 * $cos2SigmaM * $cos2SigmaM))); |
| 72 | + } while (abs($lambda - $lambdaPrev) > 1e-12 && --$iterLimit > 0); |
| 73 | + |
| 74 | + if (0 === $iterLimit) { |
| 75 | + throw new \RuntimeException('Vincenty formula failed to converge.'); |
| 76 | + } |
| 77 | + |
| 78 | + $uSq = $cosSqAlpha * (self::A * self::A - self::B * self::B) / (self::B * self::B); |
| 79 | + $Acoeff = 1 + $uSq / 16384 * (4096 + $uSq * (-768 + $uSq * (320 - 175 * $uSq))); |
| 80 | + $Bcoeff = $uSq / 1024 * (256 + $uSq * (-128 + $uSq * (74 - 47 * $uSq))); |
| 81 | + $deltaSigma = $Bcoeff * $sinSigma * ($cos2SigmaM + $Bcoeff / 4 * ($cosSigma * (-1 + 2 * $cos2SigmaM * $cos2SigmaM) |
| 82 | + - $Bcoeff / 6 * $cos2SigmaM * (-3 + 4 * $sinSigma * $sinSigma) * (-3 + 4 * $cos2SigmaM * $cos2SigmaM))); |
| 83 | + $distance = self::B * $Acoeff * ($sigma - $deltaSigma); |
| 84 | + |
| 85 | + return $distance; |
| 86 | + } |
| 87 | +} |
0 commit comments