Расстояние между двумя точками, эллипсоид. Не новость конечно, но народ интересуется часто

Очередной раз столкнулся с задачей поточнее посчитать расстояние между двумя точками на местности. Некоторые скажут "боян"... да и хрен с ними. Короче, кто до сих пор не знает как посчитать расстояние на эллипсоиде - прошу:

Координаты в десятичных градусах, типа долгота lon: 37.53 широта lat: 55.40. Результат в метрах а заодно и азимут направления движения в градусах от 0 до 360.
Тем кто с GPS работает - нужно часто. Правда еще нужно уметь перевести координаты из NMEA записи в десятичные градусы, это для пятого класса задача... ладно, сначала переведем NMEA запись координат в градусы:

function nmea2DecDegree($nmeaCoord) { // decimal(9,7), convert nmea ddmm.mmmm (dddmm.mmmm) to decimal degree
    $res = 0.0;
    $pp = strpos($nmeaCoord, '.');
    if ($pp === false || $pp <>

    $d = substr($nmeaCoord, 0, $pp-2);
    $m = substr($nmeaCoord, $pp-2);
    $res = (double)$d + ((double)$m / 60.0);
    return $res;
}

ну а теперь посчитаем расстояние в метрах и азимут в градусах:

function calcDistance ($StartLat, $StartLon, $EndLat, $EndLon) { // array, 'dist' => 'x' meters, 'bearing' => 'y' degree; from decimal degree coords
//~ 10 сантиметров точность, когда
//~ abs(lon1-lon2) <= 0.0000014
//~ abs(lat1-lat2) <= 0.0000008
    $res = array('dist'=>0.0, 'bearing'=>0.0);
    $StartLat = (double)str_replace(',', '.', $StartLat); $StartLon = (double)str_replace(',', '.', $StartLon);
    $EndLat = (double)str_replace(',', '.', $EndLat); $EndLon = (double)str_replace(',', '.', $EndLon);

    if ($StartLat == 0.0 || $StartLon == 0.0 || $EndLat == 0.0 || $EndLon == 0.0) return $res;
    if ( abs($StartLon - $EndLon) <= 0.0000014 && abs($StartLat - $EndLat) <= 0.0000008 ) return $res;

    $D2R = 0.01745329251994330; // Pi/180
    $R2D = 57.29577951308230000; // 180/Pi
/*
// http://www.pcigeomatics.com/cgi-bin/pcihlp/PROJ%7CEARTH+MODELS%7CELLIPSOIDS
// ELLIPS  Descriptor                 Semi-Major Axis (A) (metres)      Semi-Minor Axis (B) (metres)
// 12      WGS               1984  6378137.000000  6356752.314245
// 15      Krassovsky        1940  6378245.000000  6356863.018800
*/
    $a = 6378137.0; // Semi-major axis of ellipsoid in meters
//    $b = 6356752.314245;  // Semi-minor axis of ellipsoid
    $e2 = 0.00673949674233346; // 2nd eccentricity squared  Geocent_ep2 = (Geocent_a2 - Geocent_b2) / Geocent_b2

    $fdLambda = ($StartLon - $EndLon) * $D2R;
    $fdPhi = ($StartLat - $EndLat) * $D2R;
    $fPhimean = ($StartLat + $EndLat) / 2.0 * $D2R;
    $fTemp = 1 - $e2 * pow( sin($fPhimean), 2);
    $fRho = $a * (1 - $e2) / pow($fTemp, 1.5);
    $fNu = $a / sqrt(1 - $e2 * sin($fPhimean) * sin($fPhimean) );
    $fz = 2 * asin(sqrt (
        pow(sin( $fdPhi / 2.0), 2 ) + cos( $EndLat * $D2R )
        * cos( $StartLat * $D2R ) * pow(sin($fdLambda / 2.0), 2)
    ) );

    $fAlpha = asin(cos( $EndLat * $D2R ) * sin($fdLambda) / sin($fz) );
    $fR = $fRho * $fNu / ($fRho * pow(sin($fAlpha), 2) + $fNu * pow(cos($fAlpha), 2));
    $res['dist'] = $fz * $fR;

    if ($res['dist'] <= 0.999) return $res;

    $Bearing = abs($fAlpha * $R2D);
    if (($StartLat <= $EndLat) and ($StartLon > $EndLon)) $Bearing = 360 - $Bearing;
    elseif (($StartLat > $EndLat) and ($StartLon >= $EndLon)) $Bearing = 180 + $Bearing;
    elseif (($StartLat > $EndLat) and ($StartLon < $EndLon)) $Bearing = 180 - $Bearing;

    $res['bearing'] = $Bearing;
    return $res;
}