code: public static double alongTrackDistanceRad(double lat1, double lon1, double lat2, double lon2, double latQ, double lonQ) {
      // TODO: inline and share some of the trigonometric computations!
      double dist1Q = haversineFormulaRad(lat1, lon1, latQ, lonQ);
      double ctd = crossTrackDistanceRad(lat1, lon1, lat2, lon2, latQ, lonQ, dist1Q);
      return alongTrackDistanceRad(lat1, lon1, lat2, lon2, latQ, lonQ, dist1Q, ctd);
    }
docs: The along track distance is the distance from S to Q along the track from
  S to E.
  <p>
  ATD=acos(cos(dist_1Q)/cos(XTD))
  <p>
  TODO: optimize.
  
  @param lat1 Latitude of starting point in radians.
  @param lon1 Longitude of starting point in radians.
  @param lat2 Latitude of destination point in radians.
  @param lon2 Longitude of destination point in radians.
  @param latQ Latitude of query point in radians.
  @param lonQ Longitude of query point in radians.
  @return Along-track distance in radians. May be negative - this gives the
  side.