Inverting the longitude integral

[Anonymous]

Hi everyone, and thank you Charles for these awesome papers and robust library.
I've used your code to assemble together C++ tools for solving the common geodesic problems, as well as various intersections between different types of curves (currently geodesic segments and circles... while i'm also planning for offset curves, rhumb lines and some others used in some aeronautical geographic definitions). Each curve type is currently able to intersect with others, as well as parallels and meridians. For parallels, solving for a point on a geodesic line with a given latitude was straightforward, however when solving against given longitude, I'm hitting that hairy "I3" integral.

My temporary solution for intersecting stuff with a meridian is to treat it as a geodesic (works as a 180° arclength geodesic segment at least for oblate ellipsoids) and use the geodesic segment intersection routines instead. However it's itching that I could do better than this... I see that you give the series expansion coefficients for reverting the distance integral, which is used in the direct problem. Would it be possible to also invert the lambda integral in a similar manner ? Quickly browsing for Lagrange inversion give me formulas which are quite beyond my own grasp on mathematics, and I'm not versed in Maxima either. In case this kind of exercise is easy for you, could you maybe show me the way to get the coefficients of the inverted expansion ?

Thanks,
Guillaume Mirey

[FaleneLogics]

By the way, I improved my current implementation. now trying to solve for it iteratively.
Using your naming scheme, I would begin from a first guess sig12 == lam12, recomputing each time an omg12 and a new lam12 from this sigma, comparing the resulting lam12 delta to desired (dErr), and correcting back (in the majority of cases where salp0 isn't too small) omg12 with

omg12 -= dErr * (1 - _A3c);

and recovering the next iteration sig12 from this updated omg12...
(I'm thus simply ignoring the I3 term in the correction, as I have no clue how to handle it...)

This seems already quite precise and way faster than using an intersection routine with a geodesic segment in place of the meridian, although I'm still wondering whether inverting I3, in a similar manner to inverting I1, is doable ^^

[Charles Karney]

Yes the I3 series is surely invertible. But rather than messing up your
code with another compilicated series expansion, why not invert it with
Newton's method? Compute the error with the series exansion; however
use the integrand in Eq. (8) to give the derivative needed to compute
the Newton correction. This method will converge quadratically (double
the number of correct digits with each iteration). This will do a lot
better that the method given in your second message (which seems to be a
rather ad hoc approach).

[FaleneLogics]

Thanks a lot for pointing me the way. My lack of mathematical background makes those "obvious" methods hard to spot for me.
so...

omg(sig) = sqrt(1-e²cos²(bet(sig)))
         = sqrt(1-e²(1 - sin²(bet(sig))))
         = sqrt(1-e²(1 - cos²(alp0) * sin²(sig)))

omg'(sig) = e² * cos²(alp0) * sin(sig)*cos(sig) / omg(sig)

lam(sig) = omg(sig) + f * sin(alp0) * I3(sig)

lam'(sig) = omg'(sig) + f * sin(alp0) * I3'(sig)
          = omg'(sig) + f * sin(alp0) * 
                       (2-f) / (1 + (1-f) * sqrt(1 + k² + sin²(sig)))

solving for (lam(sig) - lam0) == 0,
given sig.n estimation for sig,
Newton's method next iteration :
    sig.n+1 = sig.n - (lam(sig.n) - lam0) / lam'(sig.n)

Am I on the right track ?
Thanks !
Guillaume Mirey

[Charles Karney]

A quick check says:
Need a - instead of + in the second equation for lam'(sig).
Oh and I see you fixed the iteration to use lam(sig.n) - lam0 instead of lam(sig.n).

However, a really good test is to print out the errors as the method progresses
and you need to see the error going like
1e-2
1e-4
1e-8
1e-16
1e-16 (at round off limit)

I will also repeat the check with MPFR because then you can extend the
test arbitrarily far.

[Anonymous]

You, sir, rox.
Many thanks !

[Anonymous]

In case anyone is interested, my attempt at derivating lambda above was quite wrong. (I miserad a w for omega in Charles Karney's paper). the omega derivative part is much more like :

[edit] Omg'(sig) = sin(alp0) / (cos²(sig) + sin²(alp0) * sin²(sig)))

Lambda as a function of sigma however is far less well-behaved than was the distance function : in case of geodesics which are quite close to meridians themselves, there are two cases where applying the Newton's method naively would fail :
First, when getting really close to the apex point of those almost-meridian geodesics, that derivative explodes. dividing the error by it thus would result in adding zero and no progress being made, even if the error is still significant.
Second, on the contrary, the "flat" parts of the function are really flat for almost meridians, thus the derivative is almost zero and dividing the error by it would result in correction of sigma by huge amounts (I cut those when they'd get above a quarter pi).
In those cases, I "did" revert the newton iteration to that "ad-hoc" method of mine in the second post.

Now works like a charm for all my test cases so far ;)

Note that if the geodesic line is really, really close to a meridian (sin(alp0) ~= 0), then this time I shortcut that whole iterative part altogether and force sigma2 to either the northern or southern apex of the line (whichever closer), ie an almost-polar point.
[Edit] Also, first guess for sigma is now computed from a first guess omega12 == target lambda12

Regards,
Guillaume Mirey

[Charles Karney]

Indeed... β, ω, σ, and α are all quantities defined on the auxiliary sphere. So the formulas connecting them, Eqs. (9) to (18) in Geodesics on an ellipsoid of revolution, do not depend on the ellipsoid parameters, f or e. Sorry I didn't catch this earlier.

[Anonymous]

As discussed in this thread, I need to solve for the point on a geodesic line with a given longitude. Has this method been refined? Is there code somewhere that I can use? Thanks!

[FaleneLogics]

For what it's worth, I can post the code I ended up with (it's been a while, I can't really tell the whys of everything, now).

[Note: judging by all this, in retrospect, you'd probably be better checking for intersect between two geodesics. On Earth at least, the line with constant longitude is a geodesic. Simply initialize it at equator for example, and with your desired longitude, towards azimuth 0.]

Anyway. This is the function doing the main computation work:

    bool Ellipsoid::initPointByLonDiffFromRef(const LineEqRefParams& lineEqRef, const PointOnLineParams& refPointOnLine, double dLon12, PointOnLineParams& resultingPointOnLine) const
    {
        double dLambda12Target = Maths::degToRad(dLon12);
        if (std::abs(dLambda12Target) < tol0_)
        {
            resultingPointOnLine = refPointOnLine;
            return true;
        }

        double dSinAlpha0 = lineEqRef._dSinEqAlpha;
        double dCosAlpha0 = lineEqRef._dCosEqAlpha;

        double dSigma1 = Maths::degToRad(refPointOnLine._dArcLength);
        double dSinSigma1 = refPointOnLine._dSinSigma;
        double dCosSigma1 = refPointOnLine._dCosSigma;
        double dSinOmega1 = dSinAlpha0 * dSinSigma1;
        double dCosOmega1 = dCosSigma1;
        double dOmega1 = std::atan2(dSinOmega1, dCosOmega1);
        double dLambda1 = Maths::degToRad(refPointOnLine._dLonDiff);
        double B31 = refPointOnLine._B3;

        double dK2 = lineEqRef._k2;
        double dCos2Alp0TimesE2 = _e2 * Maths::sq(dCosAlpha0);
        double dSinAlp0TimesFTimes2MinusF = _f * dSinAlpha0 * (2 - _f);

        int iIter = maxit1_;
        double dSigma12;
        double dSinArcMore, dCosArcMore;
        if (std::abs(dSinAlpha0) <= tol0_)
        {
            // Almost-meridian :
            // Target arclength becomes northern or southern apex (ie pole), whichever is closer, and get done with it
            double dArcLengthToNorthPole = Maths::angleDiff(refPointOnLine._dArcLength, 90.0);
            double dArcLengthToSouthPole = Maths::angleDiff(refPointOnLine._dArcLength, -90.0);
            double dArcLengthMore = std::abs(dArcLengthToSouthPole) < std::abs(dArcLengthToNorthPole) ?
                dArcLengthToSouthPole : dArcLengthToNorthPole;
            dSigma12 = Maths::degToRad(dArcLengthMore);
            Maths::sincosd(dArcLengthMore, dSinArcMore, dCosArcMore);
            iIter = 0;
        }
        else
        {
            // First estim : omega12 == lambda12
            double dOmega12 = dLambda12Target;
            double dOmega2 = dOmega1 + dOmega12;
            double dSinOmega2, dCosOmega2;
            Maths::sincosd(Maths::radToDeg(dOmega2), dSinOmega2, dCosOmega2);
            double dCosSig2 = dCosOmega2;
            double dSinSig2 = dSinOmega2 / dSinAlpha0;
            double dArcLength2 = Maths::atan2d(dSinSig2, dCosSig2);
            double dArcLengthMore = Maths::angleDiff(refPointOnLine._dArcLength, dArcLength2);
            dSigma12 = Maths::degToRad(dArcLengthMore);
        }
        Maths::sincosd(Maths::radToDeg(dSigma12), dSinArcMore, dCosArcMore);
        double dPreviousErr = Maths::maxValue();
        for (; iIter; iIter--)
        {
            double dSinSig2 = dSinSigma1 * dCosArcMore + dCosSigma1 * dSinArcMore;
            double dCosSig2 = dCosSigma1 * dCosArcMore - dSinSigma1 * dSinArcMore;
            double dSinOmega2 = dSinAlpha0 * dSinSig2;
            double dCosOmega2 = dCosSig2;
            //Maths::norm(dCosOmega2, dSinOmega2);
            double dOmega12 = std::atan2(dSinOmega2 * dCosOmega1 - dCosOmega2 * dSinOmega1, dCosOmega2 * dCosOmega1 + dSinOmega2 * dSinOmega1);
            double B32 = _SinCosSeries(true, dSinSig2, dCosSig2, lineEqRef._C3a, nC3_-1);
            double dLambda12 = dOmega12 + lineEqRef._A3c * (dSigma12 + (B32 - B31));

            // Error from target
            double dErr = dLambda12 - dLambda12Target;
            if (std::abs(dErr - dPreviousErr) < tol0_ || std::abs(dErr) < tol0_)
                break; // Hurray
            dPreviousErr = dErr;

            // Newton's method
            // dSigma12 -= dErr / Lam12'(sig12)
            // Where Lam12'(sig12) = Lam'(sig2)
            //                     = Omg'(sig2) - f * sin(alp0) * I3'(sig2)
            // Where I3'(sig) = (2-f) / (1 + (1-f) * sqrt(1 + k² + sin²(sig)))
            // Where Omg(sig) = atan(salp0*tan(sig))
            // => Omg'(sig) = salp0*sec²(sig)/(1 + salp0²tan²(sig)) == salp0/(cos²(sig) * (1 + salp0²tan²(sig))) == salp0/(cos²(sig) + salp0²sin²(sig)))
            //    
            double dOmg2p = dSinAlpha0 / (Maths::sq(dCosSig2) + Maths::sq(dSinAlpha0 * dSinSig2));

            // I3'(sig) = (2-f) / (1 + (1-f) * sqrt(1 + k² + sin²(sig)))
            double dI3pWithFactor = dSinAlp0TimesFTimes2MinusF / (1.0 + _f1 * std::sqrt(1.0 + dK2 * Maths::sq(dSinSig2)));

            double dLam12p = (dOmg2p - dI3pWithFactor);

            // would applying newton's method as is :
            // - switch sigma by quarter pi or more => we seem to have diverged
            if (std::abs(dLam12p) <= std::abs(dErr) * GEODESYCORE_PI * 0.25 || 
            // - or result in no correction
                std::abs(dErr / dLam12p) < tol0_)
            {
                // => we revert to correcting omega directly
                dOmega12 -= dErr * (1 - lineEqRef._A3c);
                double dOmega2 = dOmega1 + dOmega12;
                Maths::sincosd(Maths::radToDeg(dOmega2), dSinOmega2, dCosOmega2);
                dCosSig2 = dCosOmega2;
                dSinSig2 = dSinOmega2 / lineEqRef._dSinEqAlpha;
                double dArcLength2 = Maths::atan2d(dSinSig2, dCosSig2);
                double dArcLengthMore = Maths::angleDiff(refPointOnLine._dArcLength, dArcLength2);
                dSigma12 = Maths::degToRad(dArcLengthMore);
            }
            else dSigma12 -= dErr / dLam12p;

            // getting sigma12 sin/cos ready for next iteration (or even loop exit)
            Maths::sincosd(Maths::radToDeg(dSigma12), dSinArcMore, dCosArcMore);
        }   

        _onInitPointBySigmaFromRef(lineEqRef, refPointOnLine, dSigma12, dSinArcMore, dCosArcMore, resultingPointOnLine);

        double dArcLength = refPointOnLine._dArcLength + Maths::radToDeg(dSigma12);
        double dSigma2 = dSigma1 + dSigma12;

        // computes B1 and distance from eq ref
        double B12 = _SinCosSeries(true, resultingPointOnLine._dSinSigma, resultingPointOnLine._dCosSigma, lineEqRef._C1a, nC1_);
        double dDistance = (1.0 + lineEqRef._A1m1) * (dSigma2 + B12) * _b;

        resultingPointOnLine._dArcLength = 0 + dArcLength;
        resultingPointOnLine._B1 = B12;
        resultingPointOnLine._dDistance_m = 0 + dDistance;

        return true;
    }

This is the calling function, where we query intersection between geodesic line and meridian:

    eIntersectResult LineOps::computeSegmentMeridianIntersection(LineEqRefParams& lineRef, const PointOnLineParams& refPointOnLine, double dSegArcLengthBefore, double dSegArcLengthAfter,
        double dLonOfMeridian, double dEpsilonDistance, PointOnLineParams& outPointOnLine) const
    {
        PointOnLineParams arcMidPoint = refPointOnLine;
        double dDiffBeforeAfter = dSegArcLengthAfter - dSegArcLengthBefore;
        if (std::abs(dDiffBeforeAfter) > Maths::radToDeg(_pEllipsoid->tol1_))
        {
            _pEllipsoid->initPointByArcLengthFromRef(lineRef, refPointOnLine, dDiffBeforeAfter * 0.5, arcMidPoint);
        }
        if (_pEllipsoid->almostEqualsCoords(arcMidPoint._coords, PointCoords(arcMidPoint._coords.m_dLatitude, dLonOfMeridian), dEpsilonDistance))
        {
            PointOnLineParams segStart = refPointOnLine;
            if (dSegArcLengthBefore > Maths::radToDeg(_pEllipsoid->tol1_))
                _pEllipsoid->initPointByArcLengthFromRef(lineRef, refPointOnLine, -dSegArcLengthBefore, segStart);
            if (_pEllipsoid->almostEqualsCoords(segStart._coords, PointCoords(segStart._coords.m_dLatitude, dLonOfMeridian), dEpsilonDistance))
            {
                PointOnLineParams segEnd = refPointOnLine;
                if (dSegArcLengthAfter > Maths::radToDeg(_pEllipsoid->tol1_))
                    _pEllipsoid->initPointByArcLengthFromRef(lineRef, refPointOnLine, dSegArcLengthAfter, segEnd);
                if (_pEllipsoid->almostEqualsCoords(segEnd._coords, PointCoords(segEnd._coords.m_dLatitude, dLonOfMeridian), dEpsilonDistance))
                {
                    int iContainedApex = _containsApex(segStart, segEnd);
                    if (iContainedApex)
                    {
                        double dLineApexLat = _pEllipsoid->getLineApexLat(lineRef) * (double)iContainedApex;
                        PointOnLineParams segApex;
                        // call "compute intersect on tangent line" only to find apex point
                        _computeSegmentParallelIntersectOnTangentLine(lineRef, refPointOnLine, dSegArcLengthBefore, dSegArcLengthAfter, dLineApexLat,
                            dEpsilonDistance * _avgDistToArc, segApex);
                        if (_pEllipsoid->almostEqualsCoords(segApex._coords, PointCoords(segApex._coords.m_dLatitude, dLonOfMeridian), dEpsilonDistance))
                            return eIntersectResult::eIntersect_Infinity;
                        else
                        {
                            outPointOnLine = segStart;
                            return eIntersectResult::eIntersect_Two;
                        }
                    }
                    else
                        return eIntersectResult::eIntersect_Infinity;
                }
                else
                {
                    outPointOnLine = segStart;
                    return eIntersectResult::eIntersect_One;
                }
            }
        }

        double dMidPointLon = arcMidPoint._coords.m_dLongitude;
        double dLonDiff = Maths::angleDiff(dMidPointLon, Maths::longitudeNormalize(dLonOfMeridian));
        if (!_pEllipsoid->initPointByLonDiffFromRef(lineRef, arcMidPoint, dLonDiff, outPointOnLine))
            return eIntersectResult::eIntersect_Error;

        double dEpsilonArc = dEpsilonDistance * _avgDistToArc;
        double dStartArcLengthInclEps = refPointOnLine._dArcLength - dSegArcLengthBefore - dEpsilonArc;
        double dEndArcLengthInclEps = refPointOnLine._dArcLength + dSegArcLengthAfter + dEpsilonArc;

        if (outPointOnLine._dArcLength >= dStartArcLengthInclEps && outPointOnLine._dArcLength <= dEndArcLengthInclEps)
            return eIntersectResult::eIntersect_One;
        else
            return eIntersectResult::eIntersect_None;
    }

And the structures that I use, appearing as parameters in the functions above.
Hopefully self-explanatory if you understand Charles code (I have simply divided the concerns of a point on a line, from the point-agnostic params of the line, and my whole codebase is organized around this)

    // Stores every parameter of a longitude agnostic geodesic line
    // When correctly initialized, should hold coherent values for all properties but :
    // - half values (ie parameters betwwen equatorial crossings) : _dEqCrossingLonDiff _dEqCrossingDistance
    //      optional, computed-on-demand. Not-yet computed check == !Maths::isFinite(_dEqCrossingLonDiff)
    // - A2/C2 (ie differential integrals properties) : _A2m1, _C2a
    //      optional, computed-on-demand. Not-yet computed check == !Maths::isFinite(_A2m1)
    // - A4/C4 (ie area integrals properties) : _A4, _C4a
    //      optional, computed-on-demand. Not-yet computed check == !Maths::isFinite(_A4)
    struct LineEqRefParams
    {
        LineEqRefParams():
            _dEqCrossingLonDiff(Maths::NaN()), _A2m1(Maths::NaN()), _A4(Maths::NaN()) {}

        bool areHalfValuesKnown() const { return Maths::isFinite(_dEqCrossingLonDiff); }
        bool areDifferentialIntegralValuesKnown() const { return Maths::isFinite(_A2m1); }
        bool areAreaIntegralValuesKnown() const { return Maths::isFinite(_A4); }

        double getEqAzimuth() const { return Maths::azimuthNormalize(Maths::atan2d(_dSinEqAlpha, _dCosEqAlpha)); }

        double getEqCrossingLonDiff() const { return _dEqCrossingLonDiff; }
        double getEqCrossingDistance() const { return _dEqCrossingDistance; }

    private:
        double _dSinEqAlpha, _dCosEqAlpha;
        double _dEqCrossingLonDiff, _dEqCrossingDistance;
        double _k2, _eps; // values used by most integral computations
        double _A1m1, _A2m1, _A3c, _A4; // pre-computed factors for integrals (1:distance, 2:differentials, 3:longitude, 4:area)
        double _C1a[nC1_ + 1], _C1pa[nC1p_ + 1], _C2a[nC2_ + 1], _C3a[nC3_], _C4a[nC4_]; // pre-computed factors for integral sin cos series
    };

    // Stores every parameter of a point along with an azimuth uniquely qualifying a geodesic line at that point
    // When correctly initialized, should hold coherent values for all properties
    struct PointOnLineParams
    {
        PointOnLineParams() {}

        const PointCoords& getCoords() const { return _coords; }
        double getAzimuth() const { return _dAzimuth; }

        double getDistanceFromEqRef() const { return _dDistance_m; }
        double getLonDiffFromEqRef() const { return _dLonDiff; }
        double getArcLengthFromEqRef() const { return _dArcLength; }

    private:
        PointCoords _coords;
        double _dAzimuth;

        double _dSinSigma, _dCosSigma;
        double _dSinAlpha, _dCosAlpha;
        double _dDistance_m, _dLonDiff, _dArcLength;
        double _B1, _B3;    // Distance integral, Longitude integral
    };

Cheers,
Guillaume.

[Anonymous]

Guillaume,

Thank you very much for posting. I was thinking that this approach would have better performance than the intersect between two geodesics using the gnomonic projection iteratively. I'll compare and see; I need to do this with thousands of lines, and some of the lines change dynamically as the user edits them.

Best regards,
Sean

[Charles Karney]

There is certainly some overhead in using solution to the intersection problem via the gnomonic projection and coding the problem directly will give you a faster result. However, it may be that the gnomonic method is "fast enough". The gnomonic method also has the drawback that everything needs to be within one hemisphere. But, you can get around this restriction by finding an approximate solution using spherical geometry and moving the starting point for the geodesic to be close to the intersection point.