spb/internetpressure

1/** Great-circle interpolation for Pressure Front arcs (spherical linear interpolation, N points). */2export function greatCircle(from: [number, number], to: [number, number], n = 64): [number, number][] {3  const toRad = Math.PI / 180;4  const toDeg = 180 / Math.PI;5  const [lon1, lat1] = [from[0] * toRad, from[1] * toRad];6  const [lon2, lat2] = [to[0] * toRad, to[1] * toRad];7  const d = 2 * Math.asin(Math.sqrt(Math.sin((lat2 - lat1) / 2) ** 2 + Math.cos(lat1) * Math.cos(lat2) * Math.sin((lon2 - lon1) / 2) ** 2));8  if (d === 0) return [from, to];9  const pts: [number, number][] = [];10  for (let i = 0; i <= n; i++) {11    const f = i / n;12    const A = Math.sin((1 - f) * d) / Math.sin(d);13    const B = Math.sin(f * d) / Math.sin(d);14    const x = A * Math.cos(lat1) * Math.cos(lon1) + B * Math.cos(lat2) * Math.cos(lon2);15    const y = A * Math.cos(lat1) * Math.sin(lon1) + B * Math.cos(lat2) * Math.sin(lon2);16    const z = A * Math.sin(lat1) + B * Math.sin(lat2);17    const lat = Math.atan2(z, Math.sqrt(x * x + y * y)) * toDeg;18    let lon = Math.atan2(y, x) * toDeg;19    // keep the line continuous when it would cross the antimeridian20    if (pts.length) {21      const prev = pts[pts.length - 1]![0];22      if (lon - prev > 180) lon -= 360;23      else if (prev - lon > 180) lon += 360;24    }25    pts.push([lon, lat]);26  }27  return pts;28}29