sinCos() - Code Metrics - Inspection of "Switch from Vincenty to Karney for distance calcul..." - dvdoug/PHPCoord - Measure and Improve Code Quality continuously with Scrutinizer

Passed

Push — master ( 4157e1...9ef2b1 )

by Doug

created 2024-08-03 16:21 UTC

sinCos() A

↳ Parent: Project

Complexity

Conditions	5
Paths	8

Size

Total Lines	21
Code Lines	16

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
cc	5
eloc	16
nc	8
nop	1
dl	0
loc	21
rs	9.4222
c	0
b	0
f	0

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

dvdoug / PHPCoord

Push — master ( 4157e1...9ef2b1 )

sinCos() A

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like