const TWO_PI = Math.PI * 2;
const WGS84_SEMI_MAJOR_METERS = 6_378_137;
const WGS84_FLATTENING = 1 / 298.257_223_563;
const WGS84_SEMI_MINOR_METERS = WGS84_SEMI_MAJOR_METERS * (1 - WGS84_FLATTENING);

// This survey mode is deliberately regional. Radially projecting an
// Archimedean spiral onto WGS84 stays within a negligible speed/tangent error
// at this radius; a future continental mode needs an ellipsoid-arc
// reparameterization instead of silently pretending the flat formula is exact.
export const MAX_SPIRAL_RADIUS_METERS = 250_000;

export type GeodeticRadians = {
  longitude: number;
  latitude: number;
};

export type SpiralSurfaceFrame = GeodeticRadians & {
  angle: number;
  radiusMeters: number;
  radialBearing: number;
  tangentHeading: number;
};

export function normalizeRadians(value: number) {
  const normalized = ((value + Math.PI) % TWO_PI + TWO_PI) % TWO_PI - Math.PI;
  return normalized === -Math.PI ? Math.PI : normalized;
}

export function spiralRadiusAtAngle(angle: number, pitchMetersPerTurn: number) {
  const safeAngle = Math.max(0, finiteNumber(angle, "spiral_angle_invalid"));
  const pitch = positiveNumber(pitchMetersPerTurn, "spiral_pitch_invalid");
  return (pitch / TWO_PI) * safeAngle;
}

export function spiralArcLengthAtAngle(angle: number, pitchMetersPerTurn: number) {
  const safeAngle = Math.max(0, finiteNumber(angle, "spiral_angle_invalid"));
  const b = positiveNumber(pitchMetersPerTurn, "spiral_pitch_invalid") / TWO_PI;
  return (b / 2) * (
    safeAngle * Math.sqrt(1 + safeAngle * safeAngle)
    + Math.asinh(safeAngle)
  );
}

/**
 * Inverts the Archimedean spiral arc-length function. This keeps the camera's
 * configured surface speed constant instead of accelerating as the radius
 * grows. Newton is monotone here; the bounded fallback protects the first
 * samples near the origin from numerical overshoot.
 */
export function spiralAngleAtArcLength(surfaceDistanceMeters: number, pitchMetersPerTurn: number) {
  const distance = Math.max(0, finiteNumber(surfaceDistanceMeters, "spiral_distance_invalid"));
  const pitch = positiveNumber(pitchMetersPerTurn, "spiral_pitch_invalid");
  if (distance === 0) return 0;
  const b = pitch / TWO_PI;
  let angle = distance <= b ? distance / b : Math.sqrt((2 * distance) / b);
  for (let iteration = 0; iteration < 10; iteration += 1) {
    const error = spiralArcLengthAtAngle(angle, pitch) - distance;
    if (Math.abs(error) <= Math.max(1e-7, distance * 1e-12)) break;
    const derivative = b * Math.sqrt(1 + angle * angle);
    const next = angle - error / derivative;
    angle = Number.isFinite(next) && next >= 0 ? next : angle / 2;
  }
  return angle;
}

/**
 * Vincenty's WGS84 direct solution: move from a geodetic point along an
 * initial bearing by an exact ellipsoid surface distance. Longitude/latitude
 * are radians. The animation therefore remains stable at datelines and high
 * latitudes and never treats degrees as a flat Cartesian plane.
 */
export function directGeodesicDestination(
  origin: GeodeticRadians,
  initialBearing: number,
  distanceMeters: number,
): GeodeticRadians {
  const latitude1 = finiteNumber(origin.latitude, "spiral_origin_invalid");
  const longitude1 = finiteNumber(origin.longitude, "spiral_origin_invalid");
  const bearing = finiteNumber(initialBearing, "spiral_bearing_invalid");
  const distance = Math.max(0, finiteNumber(distanceMeters, "spiral_distance_invalid"));
  if (Math.abs(latitude1) > Math.PI / 2 + 1e-12) throw new Error("spiral_origin_invalid");
  if (distance === 0) return { longitude: normalizeRadians(longitude1), latitude: latitude1 };

  const sinBearing = Math.sin(bearing);
  const cosBearing = Math.cos(bearing);
  const tanReducedLatitude1 = (1 - WGS84_FLATTENING) * Math.tan(latitude1);
  const cosReducedLatitude1 = 1 / Math.sqrt(1 + tanReducedLatitude1 * tanReducedLatitude1);
  const sinReducedLatitude1 = tanReducedLatitude1 * cosReducedLatitude1;
  const sigma1 = Math.atan2(tanReducedLatitude1, cosBearing);
  const sinAlpha = cosReducedLatitude1 * sinBearing;
  const cosSquaredAlpha = 1 - sinAlpha * sinAlpha;
  const uSquared = cosSquaredAlpha
    * (WGS84_SEMI_MAJOR_METERS ** 2 - WGS84_SEMI_MINOR_METERS ** 2)
    / (WGS84_SEMI_MINOR_METERS ** 2);
  const coefficientA = 1 + (uSquared / 16_384)
    * (4096 + uSquared * (-768 + uSquared * (320 - 175 * uSquared)));
  const coefficientB = (uSquared / 1024)
    * (256 + uSquared * (-128 + uSquared * (74 - 47 * uSquared)));

  let sigma = distance / (WGS84_SEMI_MINOR_METERS * coefficientA);
  let previousSigma = Number.POSITIVE_INFINITY;
  let sinSigma = 0;
  let cosSigma = 1;
  let cosTwoSigmaMiddle = 0;
  for (let iteration = 0; iteration < 24 && Math.abs(sigma - previousSigma) > 1e-12; iteration += 1) {
    cosTwoSigmaMiddle = Math.cos(2 * sigma1 + sigma);
    sinSigma = Math.sin(sigma);
    cosSigma = Math.cos(sigma);
    const deltaSigma = coefficientB * sinSigma * (
      cosTwoSigmaMiddle
      + (coefficientB / 4) * (
        cosSigma * (-1 + 2 * cosTwoSigmaMiddle ** 2)
        - (coefficientB / 6) * cosTwoSigmaMiddle
          * (-3 + 4 * sinSigma ** 2)
          * (-3 + 4 * cosTwoSigmaMiddle ** 2)
      )
    );
    previousSigma = sigma;
    sigma = distance / (WGS84_SEMI_MINOR_METERS * coefficientA) + deltaSigma;
  }

  sinSigma = Math.sin(sigma);
  cosSigma = Math.cos(sigma);
  cosTwoSigmaMiddle = Math.cos(2 * sigma1 + sigma);
  const temporary = sinReducedLatitude1 * sinSigma
    - cosReducedLatitude1 * cosSigma * cosBearing;
  const latitude2 = Math.atan2(
    sinReducedLatitude1 * cosSigma + cosReducedLatitude1 * sinSigma * cosBearing,
    (1 - WGS84_FLATTENING) * Math.sqrt(sinAlpha ** 2 + temporary ** 2),
  );
  const lambda = Math.atan2(
    sinSigma * sinBearing,
    cosReducedLatitude1 * cosSigma - sinReducedLatitude1 * sinSigma * cosBearing,
  );
  const coefficientC = (WGS84_FLATTENING / 16) * cosSquaredAlpha
    * (4 + WGS84_FLATTENING * (4 - 3 * cosSquaredAlpha));
  const longitudeDelta = lambda - (1 - coefficientC) * WGS84_FLATTENING * sinAlpha * (
    sigma + coefficientC * sinSigma * (
      cosTwoSigmaMiddle + coefficientC * cosSigma * (-1 + 2 * cosTwoSigmaMiddle ** 2)
    )
  );
  return {
    longitude: normalizeRadians(longitude1 + longitudeDelta),
    latitude: latitude2,
  };
}

export function spiralSurfaceFrame(
  origin: GeodeticRadians,
  initialHeading: number,
  surfaceDistanceMeters: number,
  pitchMetersPerTurn: number,
): SpiralSurfaceFrame {
  const angle = spiralAngleAtArcLength(surfaceDistanceMeters, pitchMetersPerTurn);
  const radiusMeters = spiralRadiusAtAngle(angle, pitchMetersPerTurn);
  if (radiusMeters > MAX_SPIRAL_RADIUS_METERS) throw new Error("spiral_extent_limit");
  const radialBearing = normalizeRadians(initialHeading + angle);
  const destination = directGeodesicDestination(origin, radialBearing, radiusMeters);
  // In polar coordinates dr/dθ=b and r=bθ. The path tangent is therefore
  // rotated atan2(r, dr/dθ)=atan(θ) from the outward radial direction.
  const tangentHeading = normalizeRadians(radialBearing + Math.atan(angle));
  return {
    ...destination,
    angle,
    radiusMeters,
    radialBearing,
    tangentHeading,
  };
}

function finiteNumber(value: number, errorCode: string) {
  if (!Number.isFinite(value)) throw new Error(errorCode);
  return value;
}

function positiveNumber(value: number, errorCode: string) {
  const finite = finiteNumber(value, errorCode);
  if (finite <= 0) throw new Error(errorCode);
  return finite;
}