Geodesic::inverseStart() - Code Metrics - dvdoug/PHPCoord - Measure and Improve Code Quality continuously with Scrutinizer

Geodesic::inverseStart() C
last analyzed 2025-01-31 22:38 UTC

↳ Parent: Geodesic

Complexity

Conditions	14
Paths	120

Size

Total Lines	66
Code Lines	48

Duplication

Lines	0
Ratio	0 %

Code Coverage

Tests	42
CRAP Score	14.0581

Importance

Changes	1
Bugs	1	Features	0

Metric	Value
cc	14
eloc	48
c	1
b	1
f	0
nc	120
nop	7
dl	0
loc	66
ccs	42
cts	45
cp	0.9333
crap	14.0581
rs	6.1

How to fix Long Method Complexity

<?php

/**
 * PHPCoord.
 *
 * @author Doug Wright
 */
declare(strict_types=1);

namespace PHPCoord\Geometry;

use PHPCoord\CoordinateOperation\GeographicValue;
use PHPCoord\Datum\Ellipsoid;
use PHPCoord\UnitOfMeasure\Length\Length;
use PHPCoord\UnitOfMeasure\Length\Metre;

use function abs;
use function atan2;
use function cos;
use function count;
use function deg2rad;
use function fmod;
use function hypot;
use function intdiv;
use function max;
use function min;
use function range;
use function round;
use function sin;
use function sqrt;

use const M_PI;
use const PHP_FLOAT_EPSILON;
use const PHP_FLOAT_MIN;

/**
 * This is a cut-down and "PHP native" translation of the distance calculation code found inside GeographicLib,
 * Original author/copyright Charles Karney (2011-2022) released under MIT license.
 * https://geographiclib.sourceforge.io/.
 */
class Geodesic
{
    private float $a;
    private float $b;
    private float $f;
    private float $f1;
    private float $e2;
    private float $ep2;
    private float $n;

    /**
     * @var float[]
     */
    private array $a3Coefficients;

    /**
     * @var float[]
     */
    private array $c3Coefficients;

    public function __construct(Ellipsoid $ellipsoid)
    {
        $this->a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
        $this->b = $ellipsoid->getSemiMinorAxis()->asMetres()->getValue();
        $this->f = $ellipsoid->getFlattening();
        $this->f1 = 1 - $this->f;
        $this->e2 = $ellipsoid->getEccentricitySquared();
        $this->ep2 = $this->e2 / $this->f1 ** 2;
        $this->n = $this->f / (2 - $this->f);
        $this->a3Coefficients = $this->a3Coefficients();
        $this->c3Coefficients = $this->c3Coefficients();
    }

    /**
     * @return float[]
     */
    private function a3Coefficients(): array
    {
        $A3x = [];
        $coeff = [
            -3, 128,
            -2, -3, 64,
            -1, -3, -1, 16,
            3, -1, -2, 8,
            1, -1, 2,
            1, 1,
        ];
        $o = 0;
        $k = 0;
        foreach (range(5, 0, -1) as $j) {
            $m = min(5 - $j, $j);
            $A3x[$k] = polyVal($m, $coeff, $o, $this->n) / $coeff[$o + $m + 1];
            ++$k;
            $o += $m + 2;
        }

        return $A3x;
    }

    /**
     * @return float[]
     */
    private function c3Coefficients(): array
    {
        $C3x = [];

        $coeff = [
            3, 128,
            2, 5, 128,
            -1, 3, 3, 64,
            -1, 0, 1, 8,
            -1, 1, 4,
            5, 256,
            1, 3, 128,
            -3, -2, 3, 64,
            1, -3, 2, 32,
            7, 512,
            -10, 9, 384,
            5, -9, 5, 192,
            7, 512,
            -14, 7, 512,
            21, 2560,
        ];
        $o = 0;
        $k = 0;
        foreach (range(1, 5) as $l) {
            foreach (range(5, $l, -1) as $j) {
                $m = min(5 - $j, $j);
                $C3x[$k] = polyVal($m, $coeff, $o, $this->n) / $coeff[$o + $m + 1]; // @phpstan-ignore binaryOp.invalid
                ++$k;
                $o += $m + 2;
            }
        }

        return $C3x;
    }

    /**
     * @return float[]
     */
    private function c3(float $eps): array
    {
        $c = [];
        $mult = 1;
        $o = 0;
        foreach (range(1, 5) as $l) {
            $m = 5 - $l;
            $mult *= $eps;
            $c[$l] = $mult * polyVal($m, $this->c3Coefficients, $o, $eps);
            $o += $m + 1;
        }

        return $c;
    }

    /**
     * @return float[]
     */
    private function lengths(
        float $eps,
        float $sig12,
        float $ssig1,
        float $csig1,
        float $dn1,
        float $ssig2,
        float $csig2,
        float $dn2,
    ): array {
        $C1a = c1($eps);
        $C2a = c2($eps);
        $A1 = a1($eps);
        $A1 = 1 + $A1;
        $A2 = a2($eps);
        $A2 = 1 + $A2;
        $m0x = $A1 - $A2;
        $B1 = (sinCosSeries($ssig2, $csig2, $C1a) - sinCosSeries($ssig1, $csig1, $C1a));
        $B2 = (sinCosSeries($ssig2, $csig2, $C2a) - sinCosSeries($ssig1, $csig1, $C2a));
        $J12 = $m0x * $sig12 + ($A1 * $B1 - $A2 * $B2);
        $s12b = $A1 * ($sig12 + $B1);
        $m12b = ($dn2 * ($csig1 * $ssig2) - $dn1 * ($ssig1 * $csig2) - $csig1 * $csig2 * $J12);

        return [$s12b, $m12b];
    }

    /**
     * @return float[]
     */
    private function inverseStart(
        float $sbet1,
        float $cbet1,
        float $sbet2,
        float $cbet2,
        float $lam12,
        float $slam12,
        float $clam12,
    ): array {
        $sig12 = -1;
        $dnm = 1;
        $sbet12 = $sbet2 * $cbet1 - $cbet2 * $sbet1;
        $cbet12 = $cbet2 * $cbet1 + $sbet2 * $sbet1;
        $sbet12a = $sbet2 * $cbet1;
        $sbet12a += $cbet2 * $sbet1;

        $shortline = ($cbet12 >= 0 && $sbet12 < 0.5 && $cbet2 * $lam12 < 0.5);
        if ($shortline) {
            $sbetm2 = ($sbet1 + $sbet2) ** 2;
            $sbetm2 /= $sbetm2 + ($cbet1 + $cbet2) ** 2;
            $dnm = sqrt(1 + $this->ep2 * $sbetm2);
            $omg12 = $lam12 / ($this->f1 * $dnm);
            $somg12 = sin($omg12);
            $comg12 = cos($omg12);
        } else {
            $somg12 = $slam12;
            $comg12 = $clam12;
        }
        $salp1 = $cbet2 * $somg12;
        $calp1 = $comg12 >= 0 ?
          $sbet12 + $cbet2 * $sbet1 * $somg12 ** 2 / (1 + $comg12) : $sbet12a - $cbet2 * $sbet1 * $somg12 ** 2 / (1 - $comg12);

        $ssig12 = hypot($salp1, $calp1);
        $csig12 = $sbet1 * $sbet2 + $cbet1 * $cbet2 * $comg12;

        if ($shortline && $ssig12 < PHP_FLOAT_EPSILON) {
            $sig12 = atan2($ssig12, $csig12);
        } elseif (!(abs($this->n) >= 0.1 || $csig12 >= 0 || $ssig12 >= 6 * abs($this->n) * M_PI * $cbet1 ** 2)) {
            $lam12x = atan2(-$slam12, -$clam12);
            $k2 = $sbet1 ** 2 * $this->ep2;
            $eps = $k2 / (2 * (1 + sqrt(1 + $k2)) + $k2);
            $lamscale = $this->f * $cbet1 * polyVal(5, $this->a3Coefficients, 0, $eps) * M_PI;
            $betscale = $lamscale * $cbet1;
            $x = $lam12x / $lamscale;
            $y = $sbet12a / $betscale;

            if ($y > -PHP_FLOAT_EPSILON && $x > -1 - PHP_FLOAT_EPSILON) {
                $salp1 = min(1.0, -$x);
                $calp1 = -sqrt(1 - $salp1 ** 2);
            } else {
                $k = astroid($x, $y);
                $omg12a = $lamscale * ($this->f >= 0 ? -$x * $k / (1 + $k) : -$y * (1 + $k) / $k);
                $somg12 = sin($omg12a);
                $comg12 = -cos($omg12a);
                $salp1 = $cbet2 * $somg12;
                $calp1 = $sbet12a - $cbet2 * $sbet1 * $somg12 ** 2 / (1 - $comg12);
            }
        }
        if (!($salp1 <= 0)) {
            [$salp1, $calp1] = normalise($salp1, $calp1);
        } else {
            $salp1 = 1;
            $calp1 = 0;
        }

        return [$sig12, $salp1, $calp1, $dnm];
    }

    /**
     * @return float[]
     */
    private function lambda(
        float $sbet1,
        float $cbet1,
        float $dn1,
        float $sbet2,
        float $cbet2,
        float $dn2,
        float $salp1,
        float $calp1,
        float $slam120,
        float $clam120,
    ): array {
        if ($sbet1 == 0 && $calp1 == 0) {
            $calp1 = -PHP_FLOAT_MIN;
        }

        $salp0 = $salp1 * $cbet1;
        $calp0 = hypot($calp1, $salp1 * $sbet1);

        $ssig1 = $sbet1;
        $somg1 = $salp0 * $sbet1;
        $csig1 = $comg1 = $calp1 * $cbet1;
        [$ssig1, $csig1] = normalise($ssig1, $csig1);

        $calp2 = $cbet2 !== $cbet1 || abs($sbet2) !== -$sbet1 ?
            sqrt(($calp1 * $cbet1) ** 2 + ($cbet1 < -$sbet1 ?
                    ($cbet2 - $cbet1) * ($cbet1 + $cbet2) :
                    ($sbet1 - $sbet2) * ($sbet1 + $sbet2))) /
            $cbet2 : abs($calp1);

        $ssig2 = $sbet2;
        $somg2 = $salp0 * $sbet2;
        $csig2 = $comg2 = $calp2 * $cbet2;
        [$ssig2, $csig2] = normalise($ssig2, $csig2);

        $sig12 = atan2(
            max(0.0, $csig1 * $ssig2 - $ssig1 * $csig2),
            $csig1 * $csig2 + $ssig1 * $ssig2
        );

        $somg12 = max(0.0, $comg1 * $somg2 - $somg1 * $comg2);
        $comg12 = $comg1 * $comg2 + $somg1 * $somg2;
        $eta = atan2(
            $somg12 * $clam120 - $comg12 * $slam120,
            $comg12 * $clam120 + $somg12 * $slam120
        );

        $k2 = $calp0 ** 2 * $this->ep2;
        $eps = $k2 / (2 * (1 + sqrt(1 + $k2)) + $k2);
        $C3a = $this->c3($eps);
        $B312 = (sinCosSeries($ssig2, $csig2, $C3a) - sinCosSeries($ssig1, $csig1, $C3a));
        $domg12 = -$this->f * polyVal(5, $this->a3Coefficients, 0, $eps) * $salp0 * ($sig12 + $B312);
        $lam12 = $eta + $domg12;

        if ($calp2 == 0) {
            $dlam12 = -2 * $this->f1 * $dn1 / $sbet1;
        } else {
            [, $dlam12] = $this->lengths(
                $eps,
                $sig12,
                $ssig1,
                $csig1,
                $dn1,
                $ssig2,
                $csig2,
                $dn2,
            );
            $dlam12 *= $this->f1 / ($calp2 * $cbet2);
        }

        return [$lam12, $sig12, $ssig1, $csig1, $ssig2, $csig2, $eps, $dlam12];
    }

    public function distance(GeographicValue $from, GeographicValue $to): Length
    {
        $lon12 = abs(angleDiff($from->getLongitude()->asDegrees()->getValue(), $to->getLongitude()->asDegrees()->getValue()));
        $lam12 = deg2rad($lon12);
        [$slam12, $clam12] = sinCos($lon12);
        $lon12s = (180 - $lon12);

        $lat1 = $from->getLatitude()->asDegrees()->getValue();
        $lat2 = $to->getLatitude()->asDegrees()->getValue();
        $swapp = abs($lat1) < abs($lat2) ? -1 : 1;
        if ($swapp < 0) {
            [$lat2, $lat1] = [$lat1, $lat2];
        }
        $latsign = copySign(1, -$lat1);
        $lat1 *= $latsign;
        $lat2 *= $latsign;

        [$sbet1, $cbet1] = sinCos($lat1);
        $sbet1 *= $this->f1;
        [$sbet1, $cbet1] = normalise($sbet1, $cbet1);
        $cbet1 = max(PHP_FLOAT_MIN, $cbet1);

        [$sbet2, $cbet2] = sinCos($lat2);
        $sbet2 *= $this->f1;
        [$sbet2, $cbet2] = normalise($sbet2, $cbet2);
        $cbet2 = max(PHP_FLOAT_MIN, $cbet2);

        $dn1 = sqrt(1 + $this->ep2 * $sbet1 ** 2);
        $dn2 = sqrt(1 + $this->ep2 * $sbet2 ** 2);

        $meridian = $lat1 == -90 || $slam12 == 0;

        if ($meridian) {
            $calp1 = $clam12;
            $calp2 = 1.0;

            $ssig1 = $sbet1;
            $csig1 = $calp1 * $cbet1;
            $ssig2 = $sbet2;
            $csig2 = $calp2 * $cbet2;

            $sig12 = atan2(
                max(0.0, $csig1 * $ssig2 - $ssig1 * $csig2),
                $csig1 * $csig2 + $ssig1 * $ssig2
            );

            [$s12x] = $this->lengths(
                $this->n,
                $sig12,
                $ssig1,
                $csig1,
                $dn1,
                $ssig2,
                $csig2,
                $dn2,
            );

            if ($sig12 < 1) {
                $s12x *= $this->b;
            }
        } elseif ($sbet1 == 0 && $lon12s >= $this->f * 180) {
            $s12x = $this->a * $lam12;
        } else {
            [$sig12, $salp1, $calp1, $dnm] = $this->inverseStart(
                $sbet1,
                $cbet1,
                $sbet2,
                $cbet2,
                $lam12,
                $slam12,
                $clam12,
            );

            if ($sig12 >= 0) {
                $s12x = $sig12 * $this->b * $dnm;
            } else {
                $numit = 0;
                $tripn = $tripb = false;
                $salp1a = PHP_FLOAT_MIN;
                $calp1a = 1.0;
                $salp1b = PHP_FLOAT_MIN;
                $calp1b = -1.0;

                while ($numit < 20) {
                    [$v, $sig12, $ssig1, $csig1, $ssig2, $csig2,
                        $eps, $dv] = $this->lambda(
                            $sbet1,
                            $cbet1,
                            $dn1,
                            $sbet2,
                            $cbet2,
                            $dn2,
                            $salp1,
                            $calp1,
                            $slam12,
                            $clam12,
                        );
                    if ($tripb || !(abs($v) >= ($tripn ? 8 : 1) * PHP_FLOAT_EPSILON)) {
                        break;
                    }
                    if ($v > 0 && ($calp1 / $salp1 > $calp1b / $salp1b)) {
                        $salp1b = $salp1;
                        $calp1b = $calp1;
                    } elseif ($v < 0 && ($calp1 / $salp1 < $calp1a / $salp1a)) {
                        $salp1a = $salp1;
                        $calp1a = $calp1;
                    }

                    ++$numit;
                    if ($dv > 0) {
                        $dalp1 = -$v / $dv;
                        $sdalp1 = sin($dalp1);
                        $cdalp1 = cos($dalp1);
                        $nsalp1 = $salp1 * $cdalp1 + $calp1 * $sdalp1;
                        if ($nsalp1 > 0 && abs($dalp1) < M_PI) {
                            $calp1 = $calp1 * $cdalp1 - $salp1 * $sdalp1;
                            $salp1 = $nsalp1;
                            [$salp1, $calp1] = normalise($salp1, $calp1);
                            $tripn = abs($v) <= 16 * PHP_FLOAT_EPSILON;
                            continue;
                        }
                    }
                    $salp1 = ($salp1a + $salp1b) / 2;
                    $calp1 = ($calp1a + $calp1b) / 2;
                    [$salp1, $calp1] = normalise($salp1, $calp1);
                    $tripn = false;
                    $tripb = (abs($salp1a - $salp1) + ($calp1a - $calp1) < PHP_FLOAT_EPSILON
                        || abs($salp1 - $salp1b) + ($calp1 - $calp1b) < PHP_FLOAT_EPSILON);
                }
                [$s12x] = $this->lengths(
                    $eps,

                    $sig12,
                    $ssig1,

                    $csig1,

                    $dn1,
                    $ssig2,

                    $csig2,

                    $dn2,
                );

                $s12x *= $this->b;
            }
        }

        return new Metre($s12x);
    }
}

function copySign(float $x, float $y): float
{
    if ($y >= 0) {
        return abs($x);
    } else {
        return -abs($x);
    }
}

/**
 * @return float[]
 */
function normalise(float $x, float $y): array
{
    $r = hypot($x, $y);

    return [$x / $r, $y / $r];
}

/**
 * @param float[] $p
 */
function polyVal(int $N, array $p, int $s, float $x): float
{
    $y = $p[$s];
    while ($N > 0) {
        --$N;
        ++$s;
        $y = $y * $x + $p[$s];
    }

    return $y;
}

function angleDiff(float $x, float $y): float
{
    $d = ($y - $x);
    if ($d < -180) {
        $d += 360;
    }

    return $d;
}

/**
 * @return float[]
 */
function sinCos(float $x): array
{
    $r = fmod($x, 360);
    $q = (int) round($r / 90);
    $r -= 90 * $q;
    $r = deg2rad($r);
    $s = sin($r);
    $c = cos($r);
    $q %= 4;
    if ($q < 0) {
        $q += 4;
    }
    if ($q == 1) {
        [$s, $c] = [$c, -$s];
    } elseif ($q == 2) {
        [$s, $c] = [-$s, -$c];
    } elseif ($q == 3) {
        [$s, $c] = [-$c, $s];
    }

    return [$s, $c];
}

/**
 * @param float[] $c
 */
function sinCosSeries(float $sinx, float $cosx, array $c): float
{
    $k = count($c);
    $n = $k - 1;
    $ar = 2 * ($cosx - $sinx) * ($cosx + $sinx);
    $y1 = 0;
    if ($n & 1) {
        --$k;
        $y0 = $c[$k];
    } else {
        $y0 = 0;
    }
    $n = intdiv($n, 2);
    while ($n--) {
        --$k;
        $y1 = $ar * $y0 - $y1 + $c[$k];
        --$k;
        $y0 = $ar * $y1 - $y0 + $c[$k];
    }

    return 2 * $sinx * $cosx * $y0;
}

function astroid(float $x, float $y): float
{
    $p = $x ** 2;
    $q = $y ** 2;
    $r = ($p + $q - 1) / 6;
    if (!($q == 0 && $r <= 0)) {
        $S = $p * $q / 4;
        $r2 = $r ** 2;
        $r3 = $r * $r2;
        $disc = $S * ($S + 2 * $r3);
        $u = $r;
        if ($disc >= 0) {
            $T3 = $S + $r3;
            $T3 += $T3 < 0 ? -sqrt($disc) : sqrt($disc);
            $T = ($T3 ** (1 / 3));
            $u += $T != 0 ? ($T + ($r2 / $T)) : 0;
        } else {
            $ang = atan2(sqrt(-$disc), -($S + $r3));
            $u += 2 * $r * cos($ang / 3);
        }
        $v = sqrt($u ** 2 + $q);
        $uv = $u < 0 ? ($q / ($v - $u)) : $u + $v;
        $w = ($uv - $q) / (2 * $v);
        $k = $uv / (sqrt($uv + $w ** 2) + $w);
    } else {
        $k = 0;
    }

    return $k;
}

function a1(float $eps): float
{
    $coeff = [
        1, 4, 64, 0, 256,
    ];
    $t = polyVal(3, $coeff, 0, $eps ** 2) / 256;

    return ($t + $eps) / (1 - $eps);
}

/**
 * @return float[]
 */
function c1(float $eps): array
{
    $c = [];
    $coeff = [
        -1, 6, -16, 32,
        -9, 64, -128, 2048,
        9, -16, 768,
        3, -5, 512,
        -7, 1280,
        -7, 2048,
    ];
    $d = $eps;
    $o = 0;
    foreach (range(1, 6) as $l) {
        $m = intdiv(6 - $l, 2);
        $c[$l] = $d * polyVal($m, $coeff, $o, $eps ** 2) / $coeff[$o + $m + 1];
        $o += $m + 2;
        $d *= $eps;
    }

    return $c;
}

function a2(float $eps): float
{
    $coeff = [
        -11, -28, -192, 0, 256,
    ];
    $t = polyVal(3, $coeff, 0, $eps ** 2) / $coeff[4];

    return ($t - $eps) / (1 + $eps);
}

/**
 * @return float[]
 */
function c2(float $eps): array
{
    $c = [];
    $coeff = [
        1, 2, 16, 32,
        35, 64, 384, 2048,
        15, 80, 768,
        7, 35, 512,
        63, 1280,
        77, 2048,
    ];
    $d = $eps;
    $o = 0;
    foreach (range(1, 6) as $l) {
        $m = intdiv(6 - $l, 2);
        $c[$l] = $d * polyVal($m, $coeff, $o, $eps ** 2) / $coeff[$o + $m + 1];
        $o += $m + 2;
        $d *= $eps;
    }

    return $c;
}


1		<?php
2
3		/**
4		* PHPCoord.
5		*
6		* @author Doug Wright
7		*/
8		declare(strict_types=1);
9
10		namespace PHPCoord\Geometry;
11
12		use PHPCoord\CoordinateOperation\GeographicValue;
13		use PHPCoord\Datum\Ellipsoid;
14		use PHPCoord\UnitOfMeasure\Length\Length;
15		use PHPCoord\UnitOfMeasure\Length\Metre;
16
17		use function abs;
18		use function atan2;
19		use function cos;
20		use function count;
21		use function deg2rad;
22		use function fmod;
23		use function hypot;
24		use function intdiv;
25		use function max;
26		use function min;
27		use function range;
28		use function round;
29		use function sin;
30		use function sqrt;
31
32		use const M_PI;
33		use const PHP_FLOAT_EPSILON;
34		use const PHP_FLOAT_MIN;
35
36		/**
37		* This is a cut-down and "PHP native" translation of the distance calculation code found inside GeographicLib,
38		* Original author/copyright Charles Karney (2011-2022) released under MIT license.
39		* https://geographiclib.sourceforge.io/.
40		*/
41		class Geodesic
42		{
43		private float $a;
44		private float $b;
45		private float $f;
46		private float $f1;
47		private float $e2;
48		private float $ep2;
49		private float $n;
50
51		/**
52		* @var float[]
53		*/
54		private array $a3Coefficients;
55
56		/**
57		* @var float[]
58		*/
59		private array $c3Coefficients;
60	207
61		public function __construct(Ellipsoid $ellipsoid)
62	207	{
63	207	$this->a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
64	207	$this->b = $ellipsoid->getSemiMinorAxis()->asMetres()->getValue();
65	207	$this->f = $ellipsoid->getFlattening();
66	207	$this->f1 = 1 - $this->f;
67	207	$this->e2 = $ellipsoid->getEccentricitySquared();
68	207	$this->ep2 = $this->e2 / $this->f1 ** 2;
69	207	$this->n = $this->f / (2 - $this->f);
70	207	$this->a3Coefficients = $this->a3Coefficients();
71		$this->c3Coefficients = $this->c3Coefficients();
72		}
73
74		/**
75		* @return float[]
76	207	*/
77		private function a3Coefficients(): array
78	207	{
79	207	$A3x = [];
80	207	$coeff = [
81	207	-3, 128,
82	207	-2, -3, 64,
83	207	-1, -3, -1, 16,
84	207	3, -1, -2, 8,
85	207	1, -1, 2,
86	207	1, 1,
87	207	];
88	207	$o = 0;
89	207	$k = 0;
90	207	foreach (range(5, 0, -1) as $j) {
91	207	$m = min(5 - $j, $j);
92	207	$A3x[$k] = polyVal($m, $coeff, $o, $this->n) / $coeff[$o + $m + 1];
93	207	++$k;
94		$o += $m + 2;
95		}
96	207
97		return $A3x;
98		}
99
100		/**
101		* @return float[]
102	207	*/
103		private function c3Coefficients(): array
104	207	{
105		$C3x = [];
106	207
107	207	$coeff = [
108	207	3, 128,
109	207	2, 5, 128,
110	207	-1, 3, 3, 64,
111	207	-1, 0, 1, 8,
112	207	-1, 1, 4,
113	207	5, 256,
114	207	1, 3, 128,
115	207	-3, -2, 3, 64,
116	207	1, -3, 2, 32,
117	207	7, 512,
118	207	-10, 9, 384,
119	207	5, -9, 5, 192,
120	207	7, 512,
121	207	-14, 7, 512,
122	207	21, 2560,
123	207	];
124	207	$o = 0;
125	207	$k = 0;
126	207	foreach (range(1, 5) as $l) {
127	207	foreach (range(5, $l, -1) as $j) {
128	207	$m = min(5 - $j, $j);
129	207	$C3x[$k] = polyVal($m, $coeff, $o, $this->n) / $coeff[$o + $m + 1]; // @phpstan-ignore binaryOp.invalid
130	207	++$k;
131		$o += $m + 2;
132		}
133		}
134	207
135		return $C3x;
136		}
137
138		/**
139		* @return float[]
140	180	*/
141		private function c3(float $eps): array
142	180	{
143	180	$c = [];
144	180	$mult = 1;
145	180	$o = 0;
146	180	foreach (range(1, 5) as $l) {
147	180	$m = 5 - $l;
148	180	$mult *= $eps;
149	180	$c[$l] = $mult * polyVal($m, $this->c3Coefficients, $o, $eps);
150		$o += $m + 1;
151		}
152	180
153		return $c;
154		}
155
156		/**
157		* @return float[]
158	198	*/
159		private function lengths(
160		float $eps,
161		float $sig12,
162		float $ssig1,
163		float $csig1,
164		float $dn1,
165		float $ssig2,
166		float $csig2,
167		float $dn2,
168	198	): array {
169	198	$C1a = c1($eps);
170	198	$C2a = c2($eps);
171	198	$A1 = a1($eps);
172	198	$A1 = 1 + $A1;
173	198	$A2 = a2($eps);
174	198	$A2 = 1 + $A2;
175	198	$m0x = $A1 - $A2;
176	198	$B1 = (sinCosSeries($ssig2, $csig2, $C1a) - sinCosSeries($ssig1, $csig1, $C1a));
177	198	$B2 = (sinCosSeries($ssig2, $csig2, $C2a) - sinCosSeries($ssig1, $csig1, $C2a));
178	198	$J12 = $m0x * $sig12 + ($A1 * $B1 - $A2 * $B2);
179	198	$s12b = $A1 * ($sig12 + $B1);
180		$m12b = ($dn2 * ($csig1 * $ssig2) - $dn1 * ($ssig1 * $csig2) - $csig1 * $csig2 * $J12);
181	198
182		return [$s12b, $m12b];
183		}
184
185		/**
186		* @return float[]
187	180	*/
188		private function inverseStart(
189		float $sbet1,
190		float $cbet1,
191		float $sbet2,
192		float $cbet2,
193		float $lam12,
194		float $slam12,
195		float $clam12,
196	180	): array {
197	180	$sig12 = -1;
198	180	$dnm = 1;
199	180	$sbet12 = $sbet2 * $cbet1 - $cbet2 * $sbet1;
200	180	$cbet12 = $cbet2 * $cbet1 + $sbet2 * $sbet1;
201	180	$sbet12a = $sbet2 * $cbet1;
202		$sbet12a += $cbet2 * $sbet1;
203	180
204	180	$shortline = ($cbet12 >= 0 && $sbet12 < 0.5 && $cbet2 * $lam12 < 0.5);
205	54	if ($shortline) {
206	54	$sbetm2 = ($sbet1 + $sbet2) ** 2;
207	54	$sbetm2 /= $sbetm2 + ($cbet1 + $cbet2) ** 2;
208	54	$dnm = sqrt(1 + $this->ep2 * $sbetm2);
209	54	$omg12 = $lam12 / ($this->f1 * $dnm);
210	54	$somg12 = sin($omg12);
211		$comg12 = cos($omg12);
212	126	} else {
213	126	$somg12 = $slam12;
214		$comg12 = $clam12;
215	180	}
216	180	$salp1 = $cbet2 * $somg12;
217	180	$calp1 = $comg12 >= 0 ?
218		$sbet12 + $cbet2 * $sbet1 * $somg12 ** 2 / (1 + $comg12) : $sbet12a - $cbet2 * $sbet1 * $somg12 ** 2 / (1 - $comg12);
219	180
220	180	$ssig12 = hypot($salp1, $calp1);
221		$csig12 = $sbet1 * $sbet2 + $cbet1 * $cbet2 * $comg12;
222	180
223		if ($shortline && $ssig12 < PHP_FLOAT_EPSILON) {
224	180	$sig12 = atan2($ssig12, $csig12);
225	72	} elseif (!(abs($this->n) >= 0.1 \|\| $csig12 >= 0 \|\| $ssig12 >= 6 * abs($this->n) * M_PI * $cbet1 ** 2)) {
226	72	$lam12x = atan2(-$slam12, -$clam12);
227	72	$k2 = $sbet1 ** 2 * $this->ep2;
228	72	$eps = $k2 / (2 * (1 + sqrt(1 + $k2)) + $k2);
229	72	$lamscale = $this->f * $cbet1 * polyVal(5, $this->a3Coefficients, 0, $eps) * M_PI;
230	72	$betscale = $lamscale * $cbet1;
231	72	$x = $lam12x / $lamscale;
232		$y = $sbet12a / $betscale;
233	72
234	18	if ($y > -PHP_FLOAT_EPSILON && $x > -1 - PHP_FLOAT_EPSILON) {
235	18	$salp1 = min(1.0, -$x);
236		$calp1 = -sqrt(1 - $salp1 ** 2);
237	54	} else {
238	54	$k = astroid($x, $y);
239	54	$omg12a = $lamscale * ($this->f >= 0 ? -$x * $k / (1 + $k) : -$y * (1 + $k) / $k);
240	54	$somg12 = sin($omg12a);
241	54	$comg12 = -cos($omg12a);
242	54	$salp1 = $cbet2 * $somg12;
243		$calp1 = $sbet12a - $cbet2 * $sbet1 * $somg12 ** 2 / (1 - $comg12);
244		}
245	180	}
246	180	if (!($salp1 <= 0)) {
247		[$salp1, $calp1] = normalise($salp1, $calp1);
248		} else {
249		$salp1 = 1;
250		$calp1 = 0;
251		}
252	180
253		return [$sig12, $salp1, $calp1, $dnm];
254		}
255
256		/**
257		* @return float[]
258	180	*/
259		private function lambda(
260		float $sbet1,
261		float $cbet1,
262		float $dn1,
263		float $sbet2,
264		float $cbet2,
265		float $dn2,
266		float $salp1,
267		float $calp1,
268		float $slam120,
269		float $clam120,
270	180	): array {
271		if ($sbet1 == 0 && $calp1 == 0) {
272		$calp1 = -PHP_FLOAT_MIN;
273		}
274	180
275	180	$salp0 = $salp1 * $cbet1;
276		$calp0 = hypot($calp1, $salp1 * $sbet1);
277	180
278	180	$ssig1 = $sbet1;
279	180	$somg1 = $salp0 * $sbet1;
280	180	$csig1 = $comg1 = $calp1 * $cbet1;
281		[$ssig1, $csig1] = normalise($ssig1, $csig1);
282	180
283	153	$calp2 = $cbet2 !== $cbet1 \|\| abs($sbet2) !== -$sbet1 ?
284	72	sqrt(($calp1 * $cbet1) ** 2 + ($cbet1 < -$sbet1 ?
285	153	($cbet2 - $cbet1) * ($cbet1 + $cbet2) :
286	180	($sbet1 - $sbet2) * ($sbet1 + $sbet2))) /
287		$cbet2 : abs($calp1);
288	180
289	180	$ssig2 = $sbet2;
290	180	$somg2 = $salp0 * $sbet2;
291	180	$csig2 = $comg2 = $calp2 * $cbet2;
292		[$ssig2, $csig2] = normalise($ssig2, $csig2);
293	180
294	180	$sig12 = atan2(
295	180	max(0.0, $csig1 * $ssig2 - $ssig1 * $csig2),
296	180	$csig1 * $csig2 + $ssig1 * $ssig2
297		);
298	180
299	180	$somg12 = max(0.0, $comg1 * $somg2 - $somg1 * $comg2);
300	180	$comg12 = $comg1 * $comg2 + $somg1 * $somg2;
301	180	$eta = atan2(
302	180	$somg12 * $clam120 - $comg12 * $slam120,
303	180	$comg12 * $clam120 + $somg12 * $slam120
304		);
305	180
306	180	$k2 = $calp0 ** 2 * $this->ep2;
307	180	$eps = $k2 / (2 * (1 + sqrt(1 + $k2)) + $k2);
308	180	$C3a = $this->c3($eps);
309	180	$B312 = (sinCosSeries($ssig2, $csig2, $C3a) - sinCosSeries($ssig1, $csig1, $C3a));
310	180	$domg12 = -$this->f * polyVal(5, $this->a3Coefficients, 0, $eps) * $salp0 * ($sig12 + $B312);
311		$lam12 = $eta + $domg12;
312	180
313		if ($calp2 == 0) {
314		$dlam12 = -2 * $this->f1 * $dn1 / $sbet1;
315	180	} else {
316	180	[, $dlam12] = $this->lengths(
317	180	$eps,
318	180	$sig12,
319	180	$ssig1,
320	180	$csig1,
321	180	$dn1,
322	180	$ssig2,
323	180	$csig2,
324	180	$dn2,
325	180	);
326		$dlam12 = $this->f1 / ($calp2 $cbet2);
327		}
328	180
329		return [$lam12, $sig12, $ssig1, $csig1, $ssig2, $csig2, $eps, $dlam12];
330		}
331	207
332		public function distance(GeographicValue $from, GeographicValue $to): Length
333	207	{
334	207	$lon12 = abs(angleDiff($from->getLongitude()->asDegrees()->getValue(), $to->getLongitude()->asDegrees()->getValue()));
335	207	$lam12 = deg2rad($lon12);
336	207	[$slam12, $clam12] = sinCos($lon12);
337		$lon12s = (180 - $lon12);
338	207
339	207	$lat1 = $from->getLatitude()->asDegrees()->getValue();
340	207	$lat2 = $to->getLatitude()->asDegrees()->getValue();
341	207	$swapp = abs($lat1) < abs($lat2) ? -1 : 1;
342	63	if ($swapp < 0) {
343		[$lat2, $lat1] = [$lat1, $lat2];
344	207	}
345	207	$latsign = copySign(1, -$lat1);
346	207	$lat1 *= $latsign;
347		$lat2 *= $latsign;
348	207
349	207	[$sbet1, $cbet1] = sinCos($lat1);
350	207	$sbet1 *= $this->f1;
351	207	[$sbet1, $cbet1] = normalise($sbet1, $cbet1);
352		$cbet1 = max(PHP_FLOAT_MIN, $cbet1);
353	207
354	207	[$sbet2, $cbet2] = sinCos($lat2);
355	207	$sbet2 *= $this->f1;
356	207	[$sbet2, $cbet2] = normalise($sbet2, $cbet2);
357		$cbet2 = max(PHP_FLOAT_MIN, $cbet2);
358	207
359	207	$dn1 = sqrt(1 + $this->ep2 * $sbet1 ** 2);
360		$dn2 = sqrt(1 + $this->ep2 * $sbet2 ** 2);
361	207
362		$meridian = $lat1 == -90 \|\| $slam12 == 0;
363	207
364	18	if ($meridian) {
365	18	$calp1 = $clam12;
366		$calp2 = 1.0;
367	18
368	18	$ssig1 = $sbet1;
369	18	$csig1 = $calp1 * $cbet1;
370	18	$ssig2 = $sbet2;
371		$csig2 = $calp2 * $cbet2;
372	18
373	18	$sig12 = atan2(
374	18	max(0.0, $csig1 * $ssig2 - $ssig1 * $csig2),
375	18	$csig1 * $csig2 + $ssig1 * $ssig2
376		);
377	18
378	18	[$s12x] = $this->lengths(
379	18	$this->n,
380	18	$sig12,
381	18	$ssig1,
382	18	$csig1,
383	18	$dn1,
384	18	$ssig2,
385	18	$csig2,
386	18	$dn2,
387	189	);
388	9
389		if ($sig12 < 1) {
390	180	$s12x *= $this->b;
391	180	}
392	180	} elseif ($sbet1 == 0 && $lon12s >= $this->f * 180) {
393	180	$s12x = $this->a * $lam12;
394	180	} else {
395	180	[$sig12, $salp1, $calp1, $dnm] = $this->inverseStart(
396	180	$sbet1,
397	180	$cbet1,
398	180	$sbet2,
399		$cbet2,
400	180	$lam12,
401		$slam12,
402		$clam12,
403	180	);
404	180
405	180	if ($sig12 >= 0) {
406	180	$s12x = $sig12 * $this->b * $dnm;
407	180	} else {
408	180	$numit = 0;
409		$tripn = $tripb = false;
410	180	$salp1a = PHP_FLOAT_MIN;
411	180	$calp1a = 1.0;
412	180	$salp1b = PHP_FLOAT_MIN;
413	180	$calp1b = -1.0;
414	180
415	180	while ($numit < 20) {
416	180	[$v, $sig12, $ssig1, $csig1, $ssig2, $csig2,
417	180	$eps, $dv] = $this->lambda(
418	180	$sbet1,
419	180	$cbet1,
420	180	$dn1,
421	180	$sbet2,
422	180	$cbet2,
423	180	$dn2,
424	180	$salp1,
425	180	$calp1,
426		$slam12,
427	171	$clam12,
428	117	);
429	117	if ($tripb \|\| !(abs($v) >= ($tripn ? 8 : 1) * PHP_FLOAT_EPSILON)) {
430	108	break;
431	108	}
432	108	if ($v > 0 && ($calp1 / $salp1 > $calp1b / $salp1b)) {
433		$salp1b = $salp1;
434		$calp1b = $calp1;
435	171	} elseif ($v < 0 && ($calp1 / $salp1 < $calp1a / $salp1a)) {
436	171	$salp1a = $salp1;
437	171	$calp1a = $calp1;
438	171	}
439	171
440	171	++$numit;
441	171	if ($dv > 0) {
442	171	$dalp1 = -$v / $dv;
443	171	$sdalp1 = sin($dalp1);
444	171	$cdalp1 = cos($dalp1);
445	171	$nsalp1 = $salp1 * $cdalp1 + $calp1 * $sdalp1;
446	171	if ($nsalp1 > 0 && abs($dalp1) < M_PI) {
447		$calp1 = $calp1 * $cdalp1 - $salp1 * $sdalp1;
448		$salp1 = $nsalp1;
449		[$salp1, $calp1] = normalise($salp1, $calp1);
450		$tripn = abs($v) <= 16 * PHP_FLOAT_EPSILON;
451		continue;
452		}
453		}
454		$salp1 = ($salp1a + $salp1b) / 2;
455		$calp1 = ($calp1a + $calp1b) / 2;
456	180	[$salp1, $calp1] = normalise($salp1, $calp1);
457	180	$tripn = false;
458	180	$tripb = (abs($salp1a - $salp1) + ($calp1a - $calp1) < PHP_FLOAT_EPSILON
459	180	\|\| abs($salp1 - $salp1b) + ($calp1 - $calp1b) < PHP_FLOAT_EPSILON);
460	180	}
461	180	[$s12x] = $this->lengths(
462	180	$eps,
		0 ignored issues – show Comprehensibility Best Practice introduced 2024-08-03 17:21 UTC by Report Bug Copy Issue Report The variable `$eps` does not seem to be defined for all execution paths leading up to this point. Loading history...
463	180	$sig12,
464	180	$ssig1,
		0 ignored issues – show Comprehensibility Best Practice introduced 2024-08-03 17:21 UTC by Report Bug Copy Issue Report The variable `$ssig1` does not seem to be defined for all execution paths leading up to this point. Loading history...
465	180	$csig1,
		0 ignored issues – show Comprehensibility Best Practice introduced 2024-08-03 17:21 UTC by Report Bug Copy Issue Report The variable `$csig1` does not seem to be defined for all execution paths leading up to this point. Loading history...
466		$dn1,
467	180	$ssig2,
		0 ignored issues – show Comprehensibility Best Practice introduced 2024-08-03 17:21 UTC by Report Bug Copy Issue Report The variable `$ssig2` does not seem to be defined for all execution paths leading up to this point. Loading history...
468		$csig2,
		0 ignored issues – show Comprehensibility Best Practice introduced 2024-08-03 17:21 UTC by Report Bug Copy Issue Report The variable `$csig2` does not seem to be defined for all execution paths leading up to this point. Loading history...
469		$dn2,
470		);
471	207
472		$s12x *= $this->b;
473		}
474		}
475
476		return new Metre($s12x);
477	207	}
478	36	}
479
480	171	function copySign(float $x, float $y): float
481		{
482		if ($y >= 0) {
483		return abs($x);
484		} else {
485		return -abs($x);
486		}
487		}
488
489	207	/**
490		* @return float[]
491	207	*/
492		function normalise(float $x, float $y): array
493		{
494		$r = hypot($x, $y);
495
496		return [$x / $r, $y / $r];
497		}
498
499	207	/**
500	207	* @param float[] $p
501	207	*/
502	207	function polyVal(int $N, array $p, int $s, float $x): float
503	207	{
504		$y = $p[$s];
505		while ($N > 0) {
506	207	--$N;
507		++$s;
508		$y = $y * $x + $p[$s];
509		}
510
511	207	return $y;
512	207	}
513	9
514		function angleDiff(float $x, float $y): float
515		{
516	207	$d = ($y - $x);
517		if ($d < -180) {
518		$d += 360;
519		}
520
521		return $d;
522		}
523
524	207	/**
525	207	* @return float[]
526	207	*/
527	207	function sinCos(float $x): array
528	207	{
529	207	$r = fmod($x, 360);
530	207	$q = (int) round($r / 90);
531	207	$r -= 90 * $q;
532	81	$r = deg2rad($r);
533		$s = sin($r);
534	207	$c = cos($r);
535	27	$q %= 4;
536	207	if ($q < 0) {
537	99	$q += 4;
538	207	}
539	81	if ($q == 1) {
540		[$s, $c] = [$c, -$s];
541		} elseif ($q == 2) {
542	207	[$s, $c] = [-$s, -$c];
543		} elseif ($q == 3) {
544		[$s, $c] = [-$c, $s];
545		}
546
547		return [$s, $c];
548		}
549
550	198	/**
551	198	* @param float[] $c
552	198	*/
553	198	function sinCosSeries(float $sinx, float $cosx, array $c): float
554	198	{
555	198	$k = count($c);
556	198	$n = $k - 1;
557		$ar = 2 * ($cosx - $sinx) * ($cosx + $sinx);
558	180	$y1 = 0;
559		if ($n & 1) {
560	198	--$k;
561	198	$y0 = $c[$k];
562	198	} else {
563	198	$y0 = 0;
564	198	}
565	198	$n = intdiv($n, 2);
566		while ($n--) {
567		--$k;
568	198	$y1 = $ar * $y0 - $y1 + $c[$k];
569		--$k;
570		$y0 = $ar * $y1 - $y0 + $c[$k];
571		}
572
573	54	return 2 * $sinx * $cosx * $y0;
574	54	}
575	54
576	54	function astroid(float $x, float $y): float
577	54	{
578	54	$p = $x ** 2;
579	54	$q = $y ** 2;
580	54	$r = ($p + $q - 1) / 6;
581	54	if (!($q == 0 && $r <= 0)) {
582	54	$S = $p * $q / 4;
583	36	$r2 = $r ** 2;
584	36	$r3 = $r * $r2;
585	36	$disc = $S * ($S + 2 * $r3);
586	36	$u = $r;
587		if ($disc >= 0) {
588	18	$T3 = $S + $r3;
589	18	$T3 += $T3 < 0 ? -sqrt($disc) : sqrt($disc);
590		$T = ($T3 ** (1 / 3));
591	54	$u += $T != 0 ? ($T + ($r2 / $T)) : 0;
592	54	} else {
593	54	$ang = atan2(sqrt(-$disc), -($S + $r3));
594	54	$u += 2 * $r * cos($ang / 3);
595		}
596		$v = sqrt($u ** 2 + $q);
597		$uv = $u < 0 ? ($q / ($v - $u)) : $u + $v;
598		$w = ($uv - $q) / (2 * $v);
599	54	$k = $uv / (sqrt($uv + $w ** 2) + $w);
600		} else {
601		$k = 0;
602		}
603
604	198	return $k;
605	198	}
606	198
607	198	function a1(float $eps): float
608		{
609	198	$coeff = [
610		1, 4, 64, 0, 256,
611		];
612		$t = polyVal(3, $coeff, 0, $eps ** 2) / 256;
613
614		return ($t + $eps) / (1 - $eps);
615		}
616
617	198	/**
618	198	* @return float[]
619	198	*/
620	198	function c1(float $eps): array
621	198	{
622	198	$c = [];
623	198	$coeff = [
624	198	-1, 6, -16, 32,
625	198	-9, 64, -128, 2048,
626	198	9, -16, 768,
627	198	3, -5, 512,
628	198	-7, 1280,
629	198	-7, 2048,
630	198	];
631	198	$d = $eps;
632	198	$o = 0;
633		foreach (range(1, 6) as $l) {
634		$m = intdiv(6 - $l, 2);
635	198	$c[$l] = $d * polyVal($m, $coeff, $o, $eps ** 2) / $coeff[$o + $m + 1];
636		$o += $m + 2;
637		$d *= $eps;
638		}
639
640	198	return $c;
641	198	}
642	198
643	198	function a2(float $eps): float
644		{
645	198	$coeff = [
646		-11, -28, -192, 0, 256,
647		];
648		$t = polyVal(3, $coeff, 0, $eps ** 2) / $coeff[4];
649
650		return ($t - $eps) / (1 + $eps);
651		}
652
653	198	/**
654	198	* @return float[]
655	198	*/
656	198	function c2(float $eps): array
657	198	{
658	198	$c = [];
659	198	$coeff = [
660	198	1, 2, 16, 32,
661	198	35, 64, 384, 2048,
662	198	15, 80, 768,
663	198	7, 35, 512,
664	198	63, 1280,
665	198	77, 2048,
666	198	];
667	198	$d = $eps;
668	198	$o = 0;
669		foreach (range(1, 6) as $l) {
670		$m = intdiv(6 - $l, 2);
671	198	$c[$l] = $d * polyVal($m, $coeff, $o, $eps ** 2) / $coeff[$o + $m + 1];
672		$o += $m + 2;
673		$d *= $eps;
674		}
675
676		return $c;
677		}
678

dvdoug / PHPCoord

Geodesic::inverseStart() C last analyzed 2025-01-31 22:38 UTC

Complexity

Size

Duplication

Code Coverage

Importance

How to fix Long Method Complexity

Long Method

Duplication Side-by-Side

Filter issues like

Geodesic::inverseStart() C
last analyzed 2025-01-31 22:38 UTC