Vincenty inverse formula for ellipsoids

vif.php

/**
 * Calculates the geodesic distance between two points specified by radian latitude/longitude using
 * Vincenty inverse formula for ellipsoids (vif)
 */
function vif($lng1, $lat1, $lng2, $lat2)
{
    $lng1 = deg2rad($lng1);
    $lat1 = deg2rad($lat1);
    $lng2 = deg2rad($lng2);
    $lat2 = deg2rad($lat2);

    // WGS-84 ellipsoid parameters

    // length of major axis of the ellipsoid (radius at equator)
    $a = 6378137;

    // earth of minor axis of the ellipsoid (radius at the poles)
    $b = 6356752.314245;

    // flattening of the ellipsoid
    $f = 1 / 298.257223563;

    # difference in longitude
    $L = $lng2 - $lng1;

    # reduced latitude
    $U1 = atan((1 - $f) * tan($lat1));

    # reduced latitude
    $U2 = atan((1 - $f) * tan($lat2));

    $sinU1 = sin($U1);
    $cosU1 = cos($U1);
    $sinU2 = sin($U2);
    $cosU2 = cos($U2);

    $cosSqAlpha = null;
    $sinSigma = null;
    $cosSigma = null;
    $cos2SigmaM = null;
    $sigma = null;

    $lambda = $L;
    $lambdaP = 0;
    $iterLimit = 100;

    while (abs($lambda - $lambdaP) > 1e-12 & $iterLimit > 0) {
        $sinLambda = sin($lambda);
        $cosLambda = cos($lambda);
        $sinSigma = sqrt(($cosU2 * $sinLambda) * ($cosU2 * $sinLambda) +
            ($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosLambda) * ($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosLambda));
        if ($sinSigma == 0) return (0);  # Co-incident points
        $cosSigma = $sinU1 * $sinU2 + $cosU1 * $cosU2 * $cosLambda;
        $sigma = atan2($sinSigma, $cosSigma);
        $sinAlpha = $cosU1 * $cosU2 * $sinLambda / $sinSigma;
        $cosSqAlpha = 1 - $sinAlpha * $sinAlpha;
        $cos2SigmaM = $cosSigma - 2 * $sinU1 * $sinU2 / $cosSqAlpha;
        if (is_nan($cos2SigmaM)) $cos2SigmaM = 0;  # Equatorial line: cosSqAlpha=0
        $C = $f / 16 * $cosSqAlpha * (4 + $f * (4 - 3 * $cosSqAlpha));
        $lambdaP = $lambda;
        $lambda = $L + (1 - $C) * $f * $sinAlpha * ($sigma + $C * $sinSigma * ($cos2SigmaM + $C * $cosSigma * (-1 + 2 * $cos2SigmaM * $cos2SigmaM)));
        $iterLimit = $iterLimit - 1;
    }

    if ($iterLimit === 0) { # formula failed to converge
        return NAN;
    }

    $uSq = $cosSqAlpha * ($a * $a - $b * $b) / ($b * $b);
    $A = 1 + $uSq / 16384 * (4096 + $uSq * (-768 + $uSq * (320 - 175 * $uSq)));
    $B = $uSq / 1024 * (256 + $uSq * (-128 + $uSq * (74 - 47 * $uSq)));

    $deltaSigma = $B * $sinSigma * ($cos2SigmaM + $B / 4 * ($cosSigma * (-1 + 2 * $cos2SigmaM ^ 2) - $B / 6 * $cos2SigmaM * (-3 + 4 * $sinSigma ^ 2) * (-3 + 4 * $cos2SigmaM ^ 2)));

    $s = $b * $A * ($sigma - $deltaSigma);

    return $s; # Distance in km
}