Вычисление расстояния между двумя точками на планете Земля

<?php

    /**
     * Class Earth
     * Земляшка
     *
     * @copyright https://tech.yandex.ru/maps/
     */
    abstract class Earth {
        /**
         * Длина одного градуса на конкретной широте
         *
         * @param float $gps_y Широта
         *
         * @return float
         */
        static function get_x_length($gps_y) {
            return 0.00008516452991 * pow($gps_y, 3) - 0.02385660256 * pow($gps_y, 2) + 0.2589428868 * $gps_y + 107.75;
        }

        /**
         * Длина половины градуса на конкретной широте
         *
         * @param float $gps_x Широта
         *
         * @return float
         */
        static function get_y_length_05($gps_x) {
            return -0.000002053898482 * pow($gps_x, 3) + 0.0002813579278 * pow($gps_x, 2) + 0.007 * $gps_x + 55.135;
        }

        static function P($t) {
            return max(min($t, 89.999), -89.999);
        }

        static function h($t) {
            return 0 > $t ? -1 : $t > 0 ? 1 : 0;
        }

        const M = 6378137;
        const i = 0.00335281066474682000653026208241;// 1 - sqrt(1 - 0.00669437999014);
        const o = 6356752.3142451837116840090937925;// (1 - i) * M;
        const c = 1e-10;

        static function r($t) {
            return abs($t - self::h($t) * M_PI / 2) < self::c ? self::h($t) * M_PI / 2 : atan((1 - self::i) * tan($t));
        }

        /**
         * Вычисление расстояния между двумя точками на планете Земля
         *
         * @param float[] $t
         * @param float[] $s
         *
         * @return float
         */
        static function get_distance($t, $s) {
            if ($t == $s || $t[0] == $s[0] && $t[1] == $s[1]) {
                return 0;
            }
            $t[0] = self::P($t[0]);
            $s[0] = self::P($s[0]);
            if (($t[0] == 0) && ($s[0] == 0)) {
                $O = min(abs($s[1] - $t[1]), 180);

                return 2 * M_PI * self::M / 360 * abs($O);
            }

            $f = $t[0] * M_PI / 180;
            $g = $s[0] * M_PI / 180;
            $m = ($s[1] - $t[1]) * M_PI / 180;
            $D = self::r($f);
            $b = self::r($g);
            $q = $m;
            $p = $q;
            $y = 0;
            do {
                $R = sqrt(pow(cos($b) * sin($q), 2) + pow(cos($D) * sin($b) - sin($D) * cos($b) * cos($q), 2));
                $S = sin($D) * sin($b) + cos($D) * cos($b) * cos($q);
                $z = atan2($R, $S);
                $j = cos($D) * cos($b) * sin($q) / sin($z);
                $k = 1 - pow($j, 2);
                $A = cos($z) - 2 * sin($D) * sin($b) / $k;
                $B = self::i / 16 * $k * (4 + self::i * (4 - 3 * $k));
                $q = min(
                    max($m + (1 - $B) * self::i * $j * ($z + $B * sin($z) * ($A + $B * cos($z) * (-1 + 2 * $A * $A))),
                        -2 * M_PI), 2 * M_PI);
                $v = abs($p - $q);
                $p = $q;
                $y++;
            } while ($v > self::c && 10 > $y);
            $C = $k * (self::M * self::M - self::o * self::o) / (self::o * self::o);
            $E = 1 + $C * (4096 + $C * (-768 + $C * (320 - 175 * $C))) / 16384;
            $G = $C * (256 + $C * (-128 + $C * (74 - 47 * $C))) / 1024;
            $H = $G * sin($z) *
                 ($A + 0.25 * $G *
                       (cos($z) * (-1 + 2 * $A * $A) -
                        1 / 6 * $G * $A * (-3 + 4 * sin($z) * sin($z)) * (-3 + 4 * $A * $A)));

            return self::o * $E * ($z - $H);
        }
    }

    $collection = [];
    for ($i = 0; $i < 20; $i++) {
        $a = [mt_rand(20, 169), mt_rand(40, 77)];
        $b = [mt_rand(20, 169), mt_rand(40, 77)];

        print_r(a, $b, Earth::get_distance($a, $b));
    }

?>