intersectionWith method

Geographic? intersectionWith({ required double bearing, required Geographic other, required double otherBearing, })

Returns the point of intersection of two paths both defined by a position and a bearing.

The two paths are defined by:

• the current position and the bearing parameter
• other position and otherBearing both as parameters

Both bearings are measured in degrees from north (0°..360°).

The destination point is returned as a geographic position (or null if no unique intersection can be calculated).

Examples:

  const p1 = Geographic(lat: 51.8853, lon: 0.2545);
  const brng1 = 108.547;
  const p2 = Geographic(lat: 49.0034, lon: 2.5735);
  const brng2 = 32.435;

  // intersection point (lat: 50.9078°N, lon: 004.5084°E)
  final pInt = p1.spherical.intersectionWith(bearing: brng1, other: p2,
      otherBearing: brng2);

Implementation

Geographic? intersectionWith({
  required double bearing,
  required Geographic other,
  required double otherBearing,
}) {
  if (position == other) return position;

  // see www.edwilliams.org/avform.htm#Intersection

  final lat1 = position.lat.toRadians();
  final lon1 = position.lon.toRadians();
  final lat2 = other.lat.toRadians();
  final lon2 = other.lon.toRadians();
  final dlat = lat2 - lat1;
  final dlon = lon2 - lon1;

  final brng13 = bearing.toRadians();
  final brng23 = otherBearing.toRadians();

  // angular distance p1-p2
  final dst12 = 2.0 *
      asin(
        sqrt(
          sin(dlat / 2.0) * sin(dlat / 2.0) +
              cos(lat1) * cos(lat2) * sin(dlon / 2.0) * sin(dlon / 2.0),
        ),
      );

  // check for coincident points
  if (dst12.abs() < doublePrecisionEpsilon) return position;

  // initial/final bearings between points
  final cosBrngA =
      (sin(lat2) - sin(lat1) * cos(dst12)) / (sin(dst12) * cos(lat1));
  final cosBrngB =
      (sin(lat1) - sin(lat2) * cos(dst12)) / (sin(dst12) * cos(lat2));
  final brngA =
      acos(cosBrngA.clamp(-1.0, 1.0)); // protect against rounding errors
  final brngB =
      acos(cosBrngB.clamp(-1.0, 1.0)); // protect against rounding errors

  final brng12 = sin(lon2 - lon1) > 0 ? brngA : 2.0 * pi - brngA;
  final brng21 = sin(lon2 - lon1) > 0 ? 2.0 * pi - brngB : brngB;

  final ang1 = brng13 - brng12; // angle 2-1-3
  final ang2 = brng21 - brng23; // angle 1-2-3

  // check for infinite intersections
  if (sin(ang1) == 0 && sin(ang2) == 0) return null;

  // check for ambiguous intersection (antipodal/360°)
  if (sin(ang1) * sin(ang2) < 0) return null;

  final cosAng3 = -cos(ang1) * cos(ang2) + sin(ang1) * sin(ang2) * cos(dst12);

  final dst13 = atan2(
    sin(dst12) * sin(ang1) * sin(ang2),
    cos(ang2) + cos(ang1) * cosAng3,
  );

  final lat3 = asin(
    (sin(lat1) * cos(dst13) + cos(lat1) * sin(dst13) * cos(brng13))
        .clamp(-1.0, 1.0),
  );

  final dlon13 = atan2(
    sin(brng13) * sin(dst13) * cos(lat1),
    cos(dst13) - sin(lat1) * sin(lat3),
  );
  final lon3 = lon1 + dlon13;

  return Geographic(lat: lat3.toDegrees(), lon: lon3.toDegrees());
}