Vincenty Formula Jacked From Navigator PHP

vincenty.php

// Inputs in Radians so degtorad degrees first.
function vincenty($lat1, $lon1, $lat2, $lon2) {
        // Equitorial Radius
        $a = 6378137.0;
        // Polar Radius
        $b = 6356752.31424518;
        //Flattening of the ellipsoid
        $f = 0.00335281066;
        // Difference in longitude
        $L = $lon2 - $lon1;  
        $U1 = atan((1 - $f) * tan($lat1));  //U is 'reduced latitude'
        $U2 = atan((1 - $f) * tan($lat2));
        $sinU1 = sin($U1);
        $sinU2 = sin($U2);
        $cosU1 = cos($U1);
        $cosU2 = cos($U2);

        $lambda = $L;
        $lambdaP = 2 * pi();
        $i = 20;

        while (abs($lambda - $lambdaP) > 1e-12 && --$i > 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 (6)
            $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)));
        }

        if ($i == 0)
            return false;  //formula failed to converge

        $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 * $cos2SigmaM) - $B / 6 * $cos2SigmaM * (-3 + 4 * $sinSigma * $sinSigma) * (-3 + 4 * $cos2SigmaM * $cos2SigmaM)));
        $d = $b * $A * ($sigma - $deltaSigma);
        return $d;
    }