/**
 * @param  p1 Point 1
 * @param brng1 Bearing from Point 1
 * @param  p2 Point 2
 * @param brng2 bearing from Point 2
 * @returns Intersection point
 */
export const computeIntersection = (
  p1: NavFix,
  brng1: number,
  p2: NavFix,
  brng2: number
): NavFix | undefined | null => {
  if (isNaN(brng1)) throw new TypeError(`invalid brng1 ${brng1}`);
  if (isNaN(brng2)) throw new TypeError(`invalid brng2 ${brng2}`);

  const π = Math.PI;

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

  const φ1 = p1.latitude.toRadians(),
    λ1 = p1.longitude.toRadians();
  const φ2 = p2.latitude.toRadians(),
    λ2 = p2.longitude.toRadians();
  const θ13 = Number(brng1).toRadians(),
    θ23 = Number(brng2).toRadians();
  const Δφ = φ2 - φ1,
    Δλ = λ2 - λ1;

  // angular distance p1-p2
  const δ12 =
    2 *
    Math.asin(
      Math.sqrt(Math.sin(Δφ / 2) * Math.sin(Δφ / 2) + Math.cos(φ1) * Math.cos(φ2) * Math.sin(Δλ / 2) * Math.sin(Δλ / 2))
    );
  if (Math.abs(δ12) < Number.EPSILON) return p1; // coincident points

  // initial/final bearings between points
  const cosθa = (Math.sin(φ2) - Math.sin(φ1) * Math.cos(δ12)) / (Math.sin(δ12) * Math.cos(φ1));
  const cosθb = (Math.sin(φ1) - Math.sin(φ2) * Math.cos(δ12)) / (Math.sin(δ12) * Math.cos(φ2));
  const θa = Math.acos(Math.min(Math.max(cosθa, -1), 1)); // protect against rounding errors
  const θb = Math.acos(Math.min(Math.max(cosθb, -1), 1)); // protect against rounding errors

  const θ12 = Math.sin(λ2 - λ1) > 0 ? θa : 2 * π - θa;
  const θ21 = Math.sin(λ2 - λ1) > 0 ? 2 * π - θb : θb;

  const α1 = θ13 - θ12; // angle 2-1-3
  const α2 = θ21 - θ23; // angle 1-2-3

  if (Math.sin(α1) == 0 && Math.sin(α2) == 0) return undefined; // infinite intersections
  if (Math.sin(α1) * Math.sin(α2) < 0) return p2; // ambiguous intersection (antipodal/360°)

  const cosα3 = -Math.cos(α1) * Math.cos(α2) + Math.sin(α1) * Math.sin(α2) * Math.cos(δ12);

  const δ13 = Math.atan2(Math.sin(δ12) * Math.sin(α1) * Math.sin(α2), Math.cos(α2) + Math.cos(α1) * cosα3);

  const φ3 = Math.asin(
    Math.min(Math.max(Math.sin(φ1) * Math.cos(δ13) + Math.cos(φ1) * Math.sin(δ13) * Math.cos(θ13), -1), 1)
  );

  const Δλ13 = Math.atan2(Math.sin(θ13) * Math.sin(δ13) * Math.cos(φ1), Math.cos(δ13) - Math.sin(φ1) * Math.sin(φ3));
  const λ3 = λ1 + Δλ13;

  const lat = φ3.toDegrees();
  const lon = λ3.toDegrees();

  return {
    ...p1,
    latitude: lat,
    longitude: lon,
    name: 'INTC',
    isIntersection: true,
  };
};