Along Track distance

[John Coffey]

Dr Karney,
I had asked you to help me figure out the cross track distance and that worked well, but I forgot to ask you about how to calculate the Along Track distance. See '@TODO not used yet.' comment in my code snippet below.

/**
 * Calculates Cross-Track and Along-Track distances for a point
 * <code>rLocD</code> from Waypoint's <code>rLocA</code> and
 * <code>rLocB</code>
 *
 * @param rLocA  [in] point 'A' lat/lon.
 * @param rLocB  [in] point 'B' lat/lon.
 * @param rLocD  [in] point 'D' lat/lon.
 *
 * @return cross track error and along track distance (both
 *         specified in Nautical Miles.
 */
/*static*/
std::pair<double, double>
CoPilotUtils::getTrackInfoNM(
    const GeoLocation& rLocA,
    const GeoLocation& rLocB,
    const GeoLocation& rLocD) {

    // GeographicLib prefers degrees.
    GeoCoords geoA = {
        rLocA.getLat(),
        rLocA.getLon()
    };
    GeoCoords geoB = {
        rLocB.getLat(),
        rLocB.getLon()
    };
    GeoCoords geoD = {
        rLocD.getLat(),
        rLocD.getLon()
    };

    // Solve the Geodesic inverse problem from A->D
    double s12, azi1, azi2, m12, M12, M21, S12;
    auto arcAD = gGeodesic.Inverse(
        geoA.Latitude(), geoA.Longitude(),
        geoD.Latitude(), geoD.Longitude(),
        s12, azi1, azi2, m12, M12, M21, S12);
    const auto azi1AD = azi1;
    const auto m12AD = m12;
    auto arcAB = gGeodesic.Inverse(
        geoA.Latitude(), geoA.Longitude(),
        geoB.Latitude(), geoB.Longitude(),
        s12, azi1, azi2, m12, M12, M21, S12);
    const auto azi1AB = azi1;
    const auto xTrackKarneyNm = ((azi1AD - azi1AB) *
        Constants::degree() * m12AD) / Constants::nauticalmile();
    // @TODO not used yet.
    const auto aTrackKarneyNm = 0.0;
    std::pair<double, double> result = {
        xTrackKarneyNm, aTrackKarneyNm
    };
    return result;
#if 0
    // Assuming D is close to the geodesic AB, then the
    // ellipsoidal formula for cross track distance is
    // [azi1(AD) - azi1(AB)] * m12(AD). (m12 = reduced length).
    const auto dist_AD = havDistanceRads(rLocA, rLocD);
    const auto crs_AD = course(rLocA, rLocD, dist_AD);
    const auto dist_AB = havDistanceRads(rLocA, rLocB);
    const auto crs_AB = course(rLocA, rLocB, dist_AB);
    const auto xTrackRads = xTrack(dist_AD, crs_AD, crs_AB);
    const auto aTrackRads = aTrack(dist_AD, xTrackRads);
    const auto xTrackWilliamsNM = (Constants::WGS84_a() * xTrackRads / Constants::nauticalmile());
    //result.first = xTrackWilliamsNM;
    //result = {
    //    xTrackWilliamsNM,
    //    Constants::WGS84_a() * aTrackRads / Constants::nauticalmile()
    //};
    return result;
#endif
}

[Charles Karney]

Consistent with the approximation that you used for the cross track distance (i.e., that the point is close to the geodesic), you can approximate the distance along the track as the distance from A to D. (It's easy to refine these results to give the cross and along track distances with no approximation. See our earlier discussion.)

[Anonymous]

I'm sorry I did not get back to you earlier on this post as I have been using the Williams formulary for cross track and along track distance calculations. I changed to use your more accurate non spherical model to calculate the cross track distance and got that working fine today. However, I am having some trouble with the Along Track distance in the ABD triangle. In the case of the Cross Track distance the simplification I was able to successfully use was: Assuming D is close to the geodesic AB, then the ellipsoidal formula for cross track distance is [azi1(AD) - azi1(AB)] * m12(AD). (m12 = reduced length).

The problem is that the assumption that the point D is close to the AD geodesic is not true. For flight legs, it is typically within 99 nautical miles of a flight leg so the course difference between the AB geodesic and the AD geodesic can be significant. From the code snippet below you can see I am using a formula for aTrackKarneyNm.

With the fix I have - I get a very similar cross track error of 8.87NM vs -8.96NM (for your formula vs williams') - I am not worried about the sign (that only means which side of the AD geodesic, however in the same ABD triangle, I get 53NM for williams' vs -0.7254 for aTrackKarneyNm.

I wonder if you could check this formula and spot anything obvious. Maybe its radians or out by the radius of the earth of something else I might have missed.

Thanks in advance

John

// TODO: Follow up with Karney regarding the 'along track' divergence.
double s12, azi1, azi2, m12, M12, M21, S12;
/*auto arcAD = */gGeodesic.Inverse(
    geoA.Latitude(), geoA.Longitude(),
    geoD.Latitude(), geoD.Longitude(),
    s12, azi1, azi2, m12, M12, M21, S12);
const auto azi1AD = azi1;
const auto m12AD = m12;
/*auto arcAB = */gGeodesic.Inverse(
    geoA.Latitude(), geoA.Longitude(),
    geoB.Latitude(), geoB.Longitude(),
    s12, azi1, azi2, m12, M12, M21, S12);
const auto azi1AB = azi1;
const auto xTrackKarneyNm = (
    Math::AngDiff(azi1AB, azi1AD) * Constants::degree() * m12AD) /
    Constants::nauticalmile();
const auto aTrackKarneyNm = std::asin(
    std::sqrt(std::pow(std::sin(m12AD), 2) -
        std::pow(std::sin(xTrackKarneyNm), 2)) /
    std::cos(xTrackKarneyNm));
std::pair<double, double> result = {xTrackKarneyNm, aTrackKarneyNm};

[Charles Karney]

I'm sorry but I have no idea what you intend with your formula for aTrackKarneyNm. Given that m12AD is a distance, sin(m12AD) makes no sense. Similarly for sin(xTrackKarneyNm).

[Charles Karney]

Probably your best bet is to solve the interception problem: find the point C on geodesic AB which is closest to D. This is discussed here. Then the cross track distance is CD and the along track distance is AC.

[Anonymous]

Thanks, I followed your suggestion along with the example above and I have the values more or less matching the Williams' formulary values after using 4 iterations (instead of 10) as the results converge pretty quickly (I believe your method of the gnomic projection is more accurate as the assumption that the point D is near the AB geodesic is not valid. This is why I am pursuing this approach the first place. The only issue I encountered is that the s12 length for cross track distance is always positive - with the Williams' formula I am getting a sign (so I know which side of the flight path I am on). Am I doing something incorrect - do you know a way I can easily tell the side of the AB geodesic the point D is located on?

/*static*/
std::pair<double, double>
CoPilotUtils::getTrackInfoNM(
    const GeoLocation& rLocA,
    const GeoLocation& rLocB,
    const GeoLocation& rLocD) {

    // GeographicLib prefers degrees.
    GeoCoords geoA = {rLocA.getLat(), rLocA.getLon()};
    GeoCoords geoB = {rLocB.getLat(), rLocB.getLon()};
    GeoCoords geoD = {rLocD.getLat(), rLocD.getLon()};

#if 1
    static const GeographicLib::Gnomonic gn(gGeodesic);
    // intersection lat0/lon0. Initialize to a decent guess
    auto lat0 = (geoA.Latitude() + geoB.Latitude()) / 2;
    // Possibly need to deal with longitudes wrapping around
    auto lon0 = (geoA.Longitude() + geoB.Longitude()) / 2;

    // should only need a few iterations
    for (auto i = 0; i < 4; ++i) {
        double xa, ya, xb, yb;
        double xd, yd;
        gn.Forward(lat0, lon0, geoA.Latitude(), geoA.Longitude(), xa, ya);
        gn.Forward(lat0, lon0, geoB.Latitude(), geoB.Longitude(), xb, yb);
        gn.Forward(lat0, lon0, geoD.Latitude(), geoD.Longitude(), xd, yd);
        //
        vector3 va(xa, ya);
        vector3 vb(xb, yb);
        // la is homogeneous representation of line A,B
        vector3 la = va.cross(vb);
        vector3 vd(xd, yd);
        // lb is homogeneous representation of line thru D perpendicular to la
        vector3 lb(la._y, -la._x, la._x * yd - la._y * xd);
        // p0 is homogeneous representation of intersection of la and lb
        vector3 p0 = la.cross(lb);
        p0.norm();
        double lat1, lon1;
        gn.Reverse(lat0, lon0, p0._x, p0._y, lat1, lon1);
        lat0 = lat1;
        lon0 = lon1;
    }
    // distance (s12) from B to intersection lat0, lon0
    double aTrackNm;
    gGeodesic.Inverse(
        geoA.Latitude(), geoA.Longitude(),
        lat0, lon0, aTrackNm);
    aTrackNm /= Constants::nauticalmile();
    double xTrackNm;
    gGeodesic.Inverse(
        geoD.Latitude(),
        geoD.Longitude(),
        lat0, lon0, xTrackNm);
    xTrackNm /= Constants::nauticalmile();
    return {xTrackNm, aTrackNm};
#else
    // Assuming D is close to the geodesic AB, then the
    // ellipsoidal formula for cross track distance is
    // [azi1(AD) - azi1(AB)] * m12(AD). (m12 = reduced length).
    // https://sourceforge.net/p/geographiclib/discussion/1026621/thread/299518a3e4/
    const auto dist_AD = havDistanceRads(rLocA, rLocD);
    const auto crs_AD = course(rLocA, rLocD, dist_AD);
    const auto dist_AB = havDistanceRads(rLocA, rLocB);
    const auto crs_AB = course(rLocA, rLocB, dist_AB);
    const auto xTrackRads = xTrack(dist_AD, crs_AD, crs_AB);
    const auto aTrackRads = aTrack(dist_AD, xTrackRads);
    const auto aTrackWilliamsNm = (Constants::WGS84_a() *
        aTrackRads / Constants::nauticalmile());
    const auto xTrackWilliamsNm = (Constants::WGS84_a() *
        xTrackRads / Constants::nauticalmile());
    return { xTrackWilliamsNm, aTrackWilliamsNm };
#endif
}

The xTrack code is shown below (this gives signed results indicating the side of the flight line I am on.

/**
 * Calculates the 'Cross Track' distance.
 *
 * @param distRads [in] Distance from 'A to 'B' (per Ed Williams
 *                 Formulary).
 * @param crs1     [in] Course from 'A to 'B' (per Ed Williams
 *                 Formulary).
 * @param crs2     [in] Course from 'A to 'D' (per Ed Williams
 *                 Formulary).
 *
 * @return Cross Track distance (radians).
 */
/*static*/
double
CoPilotUtils::xTrack(
    const double distRads,
    const double crs1,
    const double crs2)
{
    return std::asin(std::sin(distRads) * std::sin(crs1 - crs2));
}

[Charles Karney]

Let me describe what needs to be done. (I don't want to fix your code; you should be doing this!) Having computed the geodesics from A to C (the along track path) and from C to D (the across track path), take the difference of the azimuths at C. This will be ±90° with the sign indicating which side of AB the point C lies. Draw a picture to make this clearer for you.