mjaschen / phpgeo

        $iterations = 0;

        do {

            $sinλ   = sin($λ);

            $cosλ   = cos($λ);

            $sinSqσ = ($cosU2 * $sinλ) * ($cosU2 * $sinλ)

                      + ($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosλ) * ($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosλ);

            $sinσ   = sqrt($sinSqσ);

            if ($sinσ == 0) {

                return 0;

            }

            $cosσ   = $sinU1 * $sinU2 + $cosU1 * $cosU2 * $cosλ;

            $σ      = atan2($sinσ, $cosσ);

            $sinα   = $cosU1 * $cosU2 * $sinλ / $sinσ;

            $cosSqα = 1 - $sinα * $sinα;

            if ($cosSqα == 0.0) {

                $cos2σM = 0;

            } else {

                $cos2σM = $cosσ - 2 * $sinU1 * $sinU2 / $cosSqα;

            }

            $C  = $f / 16 * $cosSqα * (4 + $f * (4 - 3 * $cosSqα));

            $λp = $λ;

            $λ  = $L + (1 - $C) * $f * $sinα * ($σ + $C * $sinσ * ($cos2σM + $C * $cosσ * (- 1 + 2 * $cos2σM * $cos2σM)));

        } while (abs($λ - $λp) > 1e-12 && ++ $iterations < 200);

        if ($iterations >= 200) {

            throw new NotConvergingException('Inverse Vincenty Formula did not converge');

        $cosU1 = cos($U1);

        $cosU2 = cos($U2);

        do {

            $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.0;

            }

            $cosSigma = $sinU1 * $sinU2 + $cosU1 * $cosU2 * $cosLambda;

            $sigma = atan2($sinSigma, $cosSigma);

            $sinAlpha = $cosU1 * $cosU2 * $sinLambda / $sinSigma;

            $cosSqAlpha = 1 - $sinAlpha * $sinAlpha;

            if ($cosSqAlpha == 0) {

                $cos2SigmaM = 0;

            } else {

                $cos2SigmaM = $cosSigma - 2 * $sinU1 * $sinU2 / $cosSqAlpha;

            }

            $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)));

        } while (abs($lambda - $lambdaP) > 1e-12 && -- $iterationLimit > 0);

        if ($iterationLimit == 0) {

            throw new NotConvergingException();