// Copyright (C) 2026 Jason M. Schwefel. All Rights Reserved.
package com.coldboreballisticsllc.scopedope.ballistics

import kotlin.math.*

// ── Constants ─────────────────────────────────────────────────────────────────

private const val STANDARD_TEMP_F  = 59.0
private const val GRAVITY_FPS2     = 32.174
private const val STD_PRESSURE_HPA = 1013.25
private const val EARTH_OMEGA      = 7.2921e-5  // Earth rotation rate rad/s
private const val KD               = 0.0764 / 288.0  // ρ₀ / (2 × 144) — drag formula constant
private const val MAX_DIST_YD      = 5000.0           // integrator hard cap

// ── Drag coefficient tables ───────────────────────────────────────────────────
// Standard G7 and G1 drag functions from McCoy, "Modern Exterior Ballistics"
// (public domain), cross-referenced with JBM Ballistics published tables.
// Values represent Cd vs Mach number for the respective standard projectile.
// Both tables end at Mach 3.0; see CdCurve.eval() for extrapolation behaviour.

// G7 — boat-tail, high-BC long-range projectile (43 entries, Mach 0–3)
private val G7_MACH = doubleArrayOf(
    0.000, 0.050, 0.100, 0.150, 0.200, 0.250, 0.300, 0.350, 0.400, 0.450,
    0.500, 0.550, 0.600, 0.650, 0.700, 0.750, 0.800, 0.850, 0.875, 0.900,
    0.925, 0.950, 0.975, 1.000, 1.025, 1.050, 1.075, 1.100, 1.150, 1.200,
    1.300, 1.400, 1.500, 1.600, 1.700, 1.800, 1.900, 2.000, 2.200, 2.400,
    2.600, 2.800, 3.000,
)
private val G7_CD = doubleArrayOf(
    0.1198, 0.1197, 0.1196, 0.1194, 0.1193, 0.1194, 0.1194, 0.1194, 0.1193, 0.1193,
    0.1194, 0.1193, 0.1194, 0.1197, 0.1202, 0.1207, 0.1215, 0.1226, 0.1242, 0.1266,
    0.1306, 0.1368, 0.1464, 0.1660, 0.2054, 0.2993, 0.3303, 0.3203, 0.2993, 0.2709,
    0.2296, 0.2049, 0.1867, 0.1729, 0.1616, 0.1520, 0.1440, 0.1372, 0.1260, 0.1168,
    0.1089, 0.1020, 0.0967,
)

// G1 — flat-base, blunt-ogive projectile (44 entries, Mach 0–3)
private val G1_MACH = doubleArrayOf(
    0.000, 0.050, 0.100, 0.150, 0.200, 0.250, 0.300, 0.350, 0.400, 0.450,
    0.500, 0.550, 0.600, 0.700, 0.725, 0.750, 0.775, 0.800, 0.825, 0.850,
    0.875, 0.900, 0.925, 0.950, 0.975, 1.000, 1.025, 1.050, 1.100, 1.150,
    1.200, 1.300, 1.400, 1.500, 1.600, 1.700, 1.800, 1.900, 2.000, 2.200,
    2.400, 2.600, 2.800, 3.000,
)
private val G1_CD = doubleArrayOf(
    0.2629, 0.2558, 0.2487, 0.2413, 0.2344, 0.2278, 0.2214, 0.2155, 0.2104, 0.2061,
    0.2032, 0.2020, 0.2034, 0.2165, 0.2230, 0.2313, 0.2417, 0.2546, 0.2706, 0.2901,
    0.3136, 0.3415, 0.3734, 0.4084, 0.4448, 0.4806, 0.5132, 0.5399, 0.5765, 0.5912,
    0.5938, 0.5716, 0.5338, 0.4950, 0.4500, 0.4054, 0.3710, 0.3417, 0.3169, 0.2752,
    0.2438, 0.2224, 0.2068, 0.1935,
)

// ── Fritsch-Carlson monotone cubic helpers ────────────────────────────────────

/** Computes Fritsch-Carlson monotone tangents for knots [xs]/[ys] (same-length arrays). */
private fun fritschCarlsonTangents(xs: DoubleArray, ys: DoubleArray): DoubleArray {
    val n = xs.size
    val d = DoubleArray(n - 1) { i -> (ys[i + 1] - ys[i]) / (xs[i + 1] - xs[i]) }
    val m = DoubleArray(n)
    m[0]     = d[0]
    m[n - 1] = d[n - 2]
    for (i in 1 until n - 1) {
        m[i] = if (d[i - 1] * d[i] <= 0.0) 0.0 else (d[i - 1] + d[i]) / 2.0
    }
    for (i in 0 until n - 1) {
        if (abs(d[i]) < 1e-12) { m[i] = 0.0; m[i + 1] = 0.0; continue }
        val alpha = m[i]     / d[i]
        val beta  = m[i + 1] / d[i]
        if (alpha < 0.0 || beta < 0.0) { m[i] = 0.0; m[i + 1] = 0.0; continue }
        val sq = alpha * alpha + beta * beta
        if (sq > 9.0) {
            val tau  = 3.0 / sqrt(sq)
            m[i]     = tau * alpha * d[i]
            m[i + 1] = tau * beta  * d[i]
        }
    }
    return m
}

/** Evaluates a cubic Hermite segment. [t] is the normalised position in [0, 1], [h] is interval width. */
private fun hermiteEval(h: Double, t: Double, y0: Double, y1: Double, m0: Double, m1: Double): Double {
    val t2 = t * t
    val t3 = t2 * t
    return (2 * t3 - 3 * t2 + 1) * y0 +
           (t3 - 2 * t2 + t)     * h * m0 +
           (-2 * t3 + 3 * t2)    * y1 +
           (t3 - t2)             * h * m1
}

/**
 * Precomputed monotone cubic Hermite (Fritsch-Carlson) spline.
 * Tangents are computed once at construction; [eval] allocates nothing per call.
 * [eval] clamps x to [xs.first(), xs.last()] — callers handle out-of-range semantics.
 * For drag curves, tables end at Mach 3.0; clamping to the last value keeps error
 * small (G7 Cd changes < 1 % between Mach 2.8–3.0) but is outside the validated range.
 */
private class MonotoneCubicSpline(private val xs: DoubleArray, private val ys: DoubleArray) {
    private val tangents: DoubleArray = fritschCarlsonTangents(xs, ys)

    fun eval(x: Double): Double {
        val xc = x.coerceIn(xs.first(), xs.last())
        var lo = 0
        var hi = xs.size - 1
        while (lo + 1 < hi) {
            val mid = (lo + hi) ushr 1
            if (xs[mid] <= xc) lo = mid else hi = mid
        }
        val h = xs[hi] - xs[lo]
        val t = (xc - xs[lo]) / h
        return hermiteEval(h, t, ys[lo], ys[hi], tangents[lo], tangents[hi])
    }
}

private val G7_CURVE = MonotoneCubicSpline(G7_MACH, G7_CD)
private val G1_CURVE = MonotoneCubicSpline(G1_MACH, G1_CD)

// ── Engine I/O types ──────────────────────────────────────────────────────────

/**
 * A single observed DOPE entry used for bias correction.
 * Atmospheric fields are the conditions when the shot was fired;
 * null values fall back to the match-day params in [BallisticsParams].
 */
data class DopeEntryInput(
    val distanceYards:          Double,
    val adjustmentMil:          Double,
    val tempF:                  Double? = null,
    val altitudeFt:             Double? = null,
    val pressureHpa:            Double? = null,
    val humidityPct:            Double? = null,
    val entryMuzzleVelocityFps: Double? = null,
)

/** Full result for a single distance: elevation, windage, spin drift, Coriolis, and aerodynamic jump. */
data class TrajectoryResult(
    val elevationMil:       Double,
    val windageMil:         Double,
    val spinDriftMil:       Double = 0.0,
    val coriolisElevMil:    Double = 0.0,
    val coriolisWindMil:    Double = 0.0,
    // Vertical displacement at the muzzle from crosswind acting on a spinning projectile.
    // Positive = up, negative = down. 0.0 when twist rate is unavailable.
    val aerodynamicJumpMil: Double = 0.0,
    // True when the bullet arrives subsonic — results are less reliable past this point.
    val isSubsonic:         Boolean = false,
)

/** All parameters describing the rifle + conditions for a calculation run. */
data class BallisticsParams(
    val muzzleVelocityFps:     Double,
    val bc:                    Double,
    val bcModel:               BcModel,
    // Optional velocity-dependent BC bands (P2). When non-empty, the integrator
    // selects the band whose minVelocityFps is the highest threshold ≤ current speed.
    // [bc] remains the fallback when no band matches (should not occur if a 0-fps band is included).
    val bcZones:               List<BcZone> = emptyList(),
    val zeroDistanceYards:     Double,
    val sightHeightIn:         Double         = 1.5,
    val altitudeFt:            Double         = 0.0,
    val tempF:                 Double         = 59.0,
    val pressureHpa:           Double         = 0.0,    // 0 = derive from altitude via ICAO
    val humidityPct:           Double         = 0.0,
    // Conditions at zeroing time. null = use today's values (no zero-shift correction).
    val zeroMuzzleVelocityFps: Double?        = null,
    val zeroTempF:             Double?        = null,
    val zeroAltitudeFt:        Double?        = null,   // null = use range-day altitudeFt
    val zeroPressureHpa:       Double?        = null,   // null = use range-day pressureHpa
    val windSpeedFps:          Double         = 0.0,
    val windAngleDeg:          Double         = 0.0,
    val twistRateIn:           Double         = 0.0,
    val twistDirection:        TwistDirection = TwistDirection.RH,
    val useCoriolis:           Boolean        = false,
    val latitudeDeg:           Double         = 0.0,
    val firingAzimuthDeg:      Double         = 0.0,
    // Bullet dimensions for Miller Sg spin-drift; 0.0 = not available (falls back to twist-only proxy)
    val bulletWeightGr:        Double         = 0.0,
    val bulletDiameterIn:      Double         = 0.0,
    val bulletLengthIn:        Double         = 0.0,
    // Form factor (i = Cd_actual / Cd_G). When bc == 0 and formFactor/dimensions are available,
    // the engine derives BC via the Pejsa sectional-density formula: BC = SD / i.
    val formFactor:            Double         = 0.0,
    // Inclined fire: positive = uphill, negative = downhill, 0 = flat.
    // Applies rifleman's rule: gEff = g × cos(θ); slant range used for mil output.
    val firingAngleDeg:        Double         = 0.0,
)

// ── Engine ────────────────────────────────────────────────────────────────────

object BallisticsEngine {

    // ── Atmosphere ────────────────────────────────────────────────────────────

    /** Temperature-adjusted speed of sound (ISA troposphere relationship). */
    fun speedOfSoundFps(tempF: Double): Double =
        1116.45 * sqrt((tempF + 459.67) / 518.67)

    /**
     * Ratio of actual air density vs ICAO standard sea-level conditions.
     *
     * If [pressureHpa] > 0, uses measured station pressure directly.
     * Otherwise derives pressure from [altitudeFt] via the ICAO tropospheric formula.
     * [humidityPct] applies the Magnus water-vapour correction (~0–1.5% effect).
     */
    fun airDensityRatio(
        altitudeFt:  Double,
        tempF:       Double,
        pressureHpa: Double = 0.0,
        humidityPct: Double = 0.0,
    ): Double {
        val tempK    = (tempF + 459.67) * 5.0 / 9.0
        val stdTempK = (STANDARD_TEMP_F + 459.67) * 5.0 / 9.0   // 288.706 K

        val actualPressure = if (pressureHpa > 0.0) pressureHpa
            else STD_PRESSURE_HPA * (1.0 - 6.87559e-6 * altitudeFt).pow(5.2561)

        // Saturation vapour pressure via Magnus formula (hPa), then partial vapour pressure
        val tempC         = (tempF - 32.0) * 5.0 / 9.0
        val eSat          = 6.1078 * exp(17.27 * tempC / (tempC + 237.3))
        val eW            = (humidityPct / 100.0) * eSat
        // Moist air is ~0.378 × eW/P less dense than dry air (water vapour lighter than N₂/O₂)
        val humidityFactor = 1.0 - 0.378 * eW / actualPressure.coerceAtLeast(1.0)

        return (actualPressure / STD_PRESSURE_HPA) * (stdTempK / tempK) * humidityFactor
    }

    // ── Wind ──────────────────────────────────────────────────────────────────

    /**
     * Decompose wind into head/cross components using the clock convention:
     * angle = direction wind comes FROM, clockwise from straight ahead.
     *
     *   0°   = headwind (from target toward shooter)
     *   90°  = right-to-left crosswind (from shooter's right, pushes bullet left)
     *   180° = tailwind
     *   270° = left-to-right crosswind (from shooter's left, pushes bullet right)
     *
     * @return Pair(headwindFps, crosswindFps)
     *   headwindFps  > 0 → opposes bullet (increases drag)
     *   crosswindFps > 0 → pushes bullet left  → positive windage (dial right)
     *   crosswindFps < 0 → pushes bullet right → negative windage (dial left)
     */
    fun windComponents(speedFps: Double, angleDeg: Double): Pair<Double, Double> {
        val rad = Math.toRadians(angleDeg)
        return speedFps * cos(rad) to speedFps * sin(rad)
    }

    // ── Numerical integrator ──────────────────────────────────────────────────

    private data class IntegrateResult(
        val dropFt:         Double,
        val tofS:           Double,
        val termVelFps:     Double,
        val coriolisWindFt: Double = 0.0,
        val coriolisElevFt: Double = 0.0,
    )

    /**
     * RK4 spatial integrator (0.5-yd steps, fixed).
     * Drag formula: a = KD × Cd(Mach) × vRel² / bcEff  (dimensionally correct).
     * Time-of-flight via trapezoidal rule; drop uses start-of-step vertical velocity.
     * [gravityFps2] may be less than GRAVITY_FPS2 for inclined fire (rifleman's rule).
     *
     * When [useCoriolis] is true, Coriolis deflection is integrated per-step using
     * step-mean bullet speed; results in [IntegrateResult.coriolisWindFt] and
     * [IntegrateResult.coriolisElevFt] (feet of lateral displacement).
     */
    private fun integrate(
        muzzleVelocityFps: Double,
        bcEff:             Double,
        bcZonesEff:        List<BcZone> = emptyList(),
        cdFn:              (Double) -> Double,
        distanceYards:     Double,
        speedOfSound:      Double,
        headwindFps:       Double  = 0.0,
        gravityFps2:       Double  = GRAVITY_FPS2,
        useCoriolis:       Boolean = false,
        sinLat:            Double  = 0.0,
        cosLat:            Double  = 0.0,
        sinAz:             Double  = 0.0,
    ): IntegrateResult {
        var v        = muzzleVelocityFps
        val stepYd   = 0.5
        val stepFt   = stepYd * 3.0
        var dropFt   = 0.0
        var dropVel  = 0.0
        var tFlight  = 0.0
        var totalYd  = 0.0
        var cWindFt  = 0.0
        var cElevFt  = 0.0
        var cWindVel = 0.0
        var cElevVel = 0.0

        fun dvds(vel: Double): Double {
            // vr = air-relative speed: Mach lookup and drag force both depend on this.
            // vel = ground speed: denominator converts dv/dt → dv/ds (chain rule: a / v_ground).
            // Tailwind (headwindFps < 0) reduces vr → less drag; headwind increases drag.
            val vr = (vel + headwindFps).coerceAtLeast(1.0)
            // Multi-zone BC: zones are pre-sorted descending by minVelocityFps.
            // Select the highest threshold that the current speed meets; fall back to bcEff.
            val bc = if (bcZonesEff.isEmpty()) bcEff
                     else bcZonesEff.firstOrNull { vel >= it.minVelocityFps }?.bc ?: bcEff
            return -KD * cdFn(vr / speedOfSound) * vr * vr / (bc * vel.coerceAtLeast(1.0))
        }

        while (totalYd < distanceYards.coerceAtMost(MAX_DIST_YD)) {
            val k1   = dvds(v)
            val k2   = dvds((v + 0.5 * stepFt * k1).coerceAtLeast(1.0))
            val k3   = dvds((v + 0.5 * stepFt * k2).coerceAtLeast(1.0))
            val k4   = dvds((v + stepFt * k3).coerceAtLeast(1.0))
            val vNew = (v + stepFt * (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0).coerceAtLeast(1.0)
            val dt   = stepFt / ((v + vNew) * 0.5)

            dropFt  += dropVel * dt + 0.5 * gravityFps2 * dt * dt
            dropVel += gravityFps2 * dt
            tFlight += dt

            if (useCoriolis) {
                val vMid   = (v + vNew) * 0.5
                val aCwind = 2.0 * EARTH_OMEGA * sinLat * vMid
                val aCelev = 2.0 * EARTH_OMEGA * cosLat * sinAz * vMid
                cWindFt  += cWindVel * dt + 0.5 * aCwind * dt * dt
                cElevFt  += cElevVel * dt + 0.5 * aCelev * dt * dt
                cWindVel += aCwind * dt
                cElevVel += aCelev * dt
            }

            v        = vNew
            totalYd += stepYd
        }

        return IntegrateResult(dropFt, tFlight, v, cWindFt, cElevFt)
    }

    // ── Coriolis ──────────────────────────────────────────────────────────────

    /**
     * Coriolis corrections for a projectile fired at [azimuthDeg] (0° = North) from
     * latitude [latitudeDeg]. Returns (elevationCorrection, windageCorrection) in mils.
     *
     * Positive elevation = dial more UP. Positive windage = dial more RIGHT.
     * In the Northern hemisphere the bullet always drifts right (positive windage).
     *
     * Uses trapezoidal velocity average (v0 + v1)/2 rather than space-mean d/t for
     * better accuracy when bullet deceleration is significant.
     *
     * Note: [calculateTrajectory] uses per-step integration (more accurate) when
     * [BallisticsParams.useCoriolis] is true. This function is retained for standalone use.
     */
    fun coriolisMils(
        latitudeDeg:         Double,
        azimuthDeg:          Double,
        distYards:           Double,
        tofS:                Double,
        muzzleVelocityFps:   Double,
        terminalVelocityFps: Double,
    ): Pair<Double, Double> {
        if (distYards <= 0.0 || tofS <= 0.0) return 0.0 to 0.0
        val latRad = Math.toRadians(latitudeDeg)
        val azRad  = Math.toRadians(azimuthDeg)
        val vAvg   = (muzzleVelocityFps + terminalVelocityFps) * 0.5
        val mrad   = distYards * 36.0 / 1000.0

        val windFt = EARTH_OMEGA * sin(latRad) * vAvg * tofS * tofS
        val elevFt = EARTH_OMEGA * cos(latRad) * sin(azRad) * vAvg * tofS * tofS

        return -(elevFt * 12.0) / mrad to (windFt * 12.0) / mrad
    }

    // ── Spin drift ────────────────────────────────────────────────────────────

    /**
     * Miller (1983) gyroscopic stability factor.
     * Returns Sg clamped to [1.0, 5.0]; returns 0 if any dimension is missing.
     * T = twist rate in inches/revolution, D = diameter in inches, L = length in inches.
     */
    private fun millerSg(
        twistRateIn:      Double,
        bulletWeightGr:   Double,
        bulletDiameterIn: Double,
        bulletLengthIn:   Double,
        muzzleVelFps:     Double,
    ): Double {
        if (bulletDiameterIn <= 0.0 || bulletLengthIn <= 0.0 || bulletWeightGr <= 0.0) return 0.0
        val l  = bulletLengthIn / bulletDiameterIn
        val sg = 30.0 * bulletWeightGr / (twistRateIn * twistRateIn * bulletDiameterIn.pow(3) * l * (1.0 + l * l))
        val vf = if (muzzleVelFps > 0.0) (muzzleVelFps / 2800.0).pow(1.0 / 3.0) else 1.0
        // Return raw value; cap upper end at 5 but do not clamp below 1.0 so the caller
        // can detect Sg < 1.0 (dynamically unstable — spin drift is undefined).
        return (sg * vf).coerceAtMost(5.0)
    }

    /**
     * Spin drift windage correction using Litz approximation with Miller Sg.
     * Falls back to twist-only Sg proxy when bullet dimensions are not provided (= 0).
     * Returns 0 when [twistRateIn] ≤ 0 (twist rate not configured).
     * Positive = bullet drifts right (RH twist); negative = left (LH twist).
     */
    fun spinDriftMils(
        twistRateIn:      Double,
        twistDirection:   TwistDirection,
        tofS:             Double,
        distYards:        Double,
        bulletWeightGr:   Double = 0.0,
        bulletDiameterIn: Double = 0.0,
        bulletLengthIn:   Double = 0.0,
        muzzleVelFps:     Double = 0.0,
    ): Double {
        if (twistRateIn <= 0.0 || distYards <= 0.0 || tofS <= 0.0) return 0.0
        val rawSg = millerSg(twistRateIn, bulletWeightGr, bulletDiameterIn, bulletLengthIn, muzzleVelFps)
        // Sg < 1.0 → dynamically unstable; spin drift is physically undefined.
        if (rawSg in 0.01..0.999) return 0.0
        val sg = rawSg.takeIf { it > 0.0 }
            ?: (14.0 / twistRateIn).coerceIn(1.0, 3.0)
        val driftIn   = 1.25 * sg * tofS.pow(1.83)
        val driftMils = driftIn / (distYards * 36.0 / 1000.0)
        return if (twistDirection == TwistDirection.RH) driftMils else -driftMils
    }

    // ── Aerodynamic jump ──────────────────────────────────────────────────────

    /**
     * Aerodynamic jump — vertical displacement caused by a crosswind acting on a spinning
     * projectile at the muzzle (Litz approximation, Applied Ballistics §4).
     *
     * Returns 0.0 when [twistRateIn] ≤ 0 (twist not configured) or Sg < 1 (unstable).
     * Positive = up, negative = down. Sign inverts for LH twist.
     *
     * Note: AJ is a muzzle-exit effect and is **independent of range** — it shifts the
     * entire trajectory uniformly. The per-distance denominator converts to milliradians.
     */
    fun aerodynamicJumpMils(
        twistRateIn:       Double,
        twistDirection:    TwistDirection,
        crosswindFps:      Double,
        muzzleVelocityFps: Double,
        bulletWeightGr:    Double = 0.0,
        bulletDiameterIn:  Double = 0.0,
        bulletLengthIn:    Double = 0.0,
    ): Double {
        if (twistRateIn <= 0.0 || crosswindFps == 0.0) return 0.0
        val sg = millerSg(twistRateIn, bulletWeightGr, bulletDiameterIn, bulletLengthIn, muzzleVelocityFps)
        if (sg in 0.01..0.999) return 0.0
        val sgEff = sg.takeIf { it > 0.0 } ?: (14.0 / twistRateIn).coerceIn(1.0, 3.0)
        // Litz approximation: AJ (mils) = 1.25 × Sg × Vcw / V0
        // V_cw/V_0 is dimensionless; result is already milliradians (angular, range-independent).
        // No /mrad term — AJ is a muzzle-exit angular deflection, not a physical displacement.
        val aj = 1.25 * sgEff * crosswindFps / muzzleVelocityFps.coerceAtLeast(1.0)
        // Positive crosswind (R→L) with RH twist produces negative AJ (bullet pushed down).
        return if (twistDirection == TwistDirection.RH) -aj else aj
    }

    // ── Elevation + windage helpers ───────────────────────────────────────────

    private fun trajectoryMils(
        dropD:         Double,
        dropZ:         Double,
        distYards:     Double,
        zeroYards:     Double,
        sightHeightIn: Double,
    ): Double {
        if (distYards <= 0.0 || zeroYards <= 0.0) return 0.0
        val hFt      = sightHeightIn / 12.0
        val scale    = distYards / zeroYards
        val neededFt = dropD + hFt - scale * (dropZ + hFt)
        return (neededFt * 12.0) / (distYards * 36.0 / 1000.0)
    }

    // ── Public API ────────────────────────────────────────────────────────────

    /**
     * Validate [BallisticsParams] before integration.
     * Returns a non-null error string if the params cannot produce a valid trajectory.
     * Pass [distancesYards] to also validate that no distance exceeds the integrator cap.
     */
    fun validateParams(params: BallisticsParams, distancesYards: List<Double> = emptyList()): String? = when {
        params.muzzleVelocityFps <= 0  -> "Muzzle velocity must be > 0"
        params.zeroDistanceYards <= 0  -> "Zero distance must be > 0"
        params.bc <= 0.0 && (params.formFactor <= 0.0 ||
                              params.bulletDiameterIn <= 0.0 ||
                              params.bulletWeightGr <= 0.0)
                                       -> "BC or complete 4DOF dimensions required"
        params.pressureHpa < 0         -> "Pressure cannot be negative"
        params.humidityPct !in 0.0..100.0 -> "Humidity must be 0–100 %"
        distancesYards.any { it > MAX_DIST_YD }
                                       -> "Distance exceeds integrator cap of ${MAX_DIST_YD.toInt()} yd"
        else                           -> null
    }

    /**
     * Compute elevation and windage for each distance in [distancesYards].
     * No DOPE bias correction — pure physics.
     *
     * Returns an empty map when [validateParams] returns a non-null error.
     *
     * The zero integration uses [BallisticsParams.zeroMuzzleVelocityFps] and
     * [BallisticsParams.zeroTempF] when present, so the zero point reflects the
     * conditions under which the rifle was actually zeroed.
     */
    fun calculateTrajectory(
        params:         BallisticsParams,
        distancesYards: List<Double>,
    ): Map<Double, TrajectoryResult> {
        if (validateParams(params, distancesYards) != null) return emptyMap()
        val (hw, cw) = windComponents(params.windSpeedFps, params.windAngleDeg)
        // Published G7 BCs are calibrated to JBM's G7 drag function, not McCoy's.
        // Scale McCoy G7 Cd by 1.26 to align with JBM — empirically validated for BC 0.25–0.35, MV 2400–3200fps.
        // G1 is left unscaled; G1 BCs use a different reference and the calibration was not validated.
        val cdFn: (Double) -> Double = if (params.bcModel == BcModel.G7) { mach -> G7_CURVE.eval(mach) * 1.26 } else G1_CURVE::eval

        // When no measured BC is provided but form factor + bullet dimensions are available,
        // derive BC from sectional density: BC = (W_gr / 7000 / D_in²) / formFactor (Pejsa method).
        val bcBase = if (params.bc <= 0.0 &&
                         params.formFactor > 0.0 &&
                         params.bulletDiameterIn > 0.0 &&
                         params.bulletWeightGr > 0.0) {
            val sd = (params.bulletWeightGr / 7000.0) / (params.bulletDiameterIn * params.bulletDiameterIn)
            sd / params.formFactor
        } else params.bc

        // Range-day air density (measured pressure + humidity when available)
        val densityRatio      = airDensityRatio(params.altitudeFt, params.tempF, params.pressureHpa, params.humidityPct)
        val bcEff             = bcBase / densityRatio
        val speedOfSoundRange = speedOfSoundFps(params.tempF)
        // Pre-scale zone BCs by density ratio, sorted descending for O(1) lookup in dvds().
        val bcZonesEff        = params.bcZones.sortedByDescending { it.minVelocityFps }
                                    .map { BcZone(it.minVelocityFps, it.bc / densityRatio) }

        // Zero integration: use MV and air density at zeroing conditions.
        // Wind is assumed calm at zero — we don't know what it was.
        val zeroMv           = params.zeroMuzzleVelocityFps ?: params.muzzleVelocityFps
        val zeroTF           = params.zeroTempF             ?: params.tempF
        val zeroDensityRatio = airDensityRatio(
            params.zeroAltitudeFt  ?: params.altitudeFt,
            zeroTF,
            params.zeroPressureHpa ?: params.pressureHpa,
            params.humidityPct,
        )
        val bcEffAtZero      = bcBase / zeroDensityRatio
        val bcZonesEffAtZero = params.bcZones.sortedByDescending { it.minVelocityFps }
                                   .map { BcZone(it.minVelocityFps, it.bc / zeroDensityRatio) }
        val speedOfSoundZero = speedOfSoundFps(zeroTF)
        // Rifleman's rule: reduce effective gravity by cos(θ) for inclined fire.
        // Zero is assumed established on flat ground, so no angle correction for zero integration.
        val gEff       = GRAVITY_FPS2 * cos(Math.toRadians(params.firingAngleDeg))
        val zeroResult = integrate(zeroMv, bcEffAtZero, bcZonesEffAtZero, cdFn, params.zeroDistanceYards, speedOfSoundZero)
        val dropZ      = zeroResult.dropFt

        // Pre-compute Coriolis trig once — avoids repeated sin/cos inside the per-distance loop.
        val sinLat = if (params.useCoriolis) sin(Math.toRadians(params.latitudeDeg)) else 0.0
        val cosLat = if (params.useCoriolis) cos(Math.toRadians(params.latitudeDeg)) else 0.0
        val sinAz  = if (params.useCoriolis) sin(Math.toRadians(params.firingAzimuthDeg)) else 0.0

        return distancesYards.associateWith { dist ->
            if (dist == 0.0) {
                TrajectoryResult(0.0, 0.0)
            } else {
                val result     = integrate(params.muzzleVelocityFps, bcEff, bcZonesEff, cdFn, dist, speedOfSoundRange, hw, gEff,
                                           params.useCoriolis, sinLat, cosLat, sinAz)
                val mrad       = dist * 36.0 / 1000.0
                val isSubsonic = result.termVelFps < speedOfSoundRange
                val coriolisElev = if (params.useCoriolis && mrad > 0) -(result.coriolisElevFt * 12.0) / mrad else 0.0
                val coriolisWind = if (params.useCoriolis && mrad > 0) (result.coriolisWindFt * 12.0) / mrad else 0.0

                // Didion lag-rule windage: correct for spin-stabilized bullets where gyroscopic
                // precession causes the bullet to track wind direction far more efficiently
                // than simple aerodynamic drag coupling.
                val tVacuum    = dist * 3.0 / maxOf(params.muzzleVelocityFps, 1.0)
                val windageMil = if (mrad > 0) cw * (result.tofS - tVacuum) * 12.0 / mrad else 0.0

                TrajectoryResult(
                    elevationMil       = trajectoryMils(result.dropFt, dropZ, dist, params.zeroDistanceYards, params.sightHeightIn),
                    windageMil         = windageMil,
                    spinDriftMil       = spinDriftMils(params.twistRateIn, params.twistDirection, result.tofS, dist,
                                                       params.bulletWeightGr, params.bulletDiameterIn,
                                                       params.bulletLengthIn, params.muzzleVelocityFps),
                    coriolisElevMil    = coriolisElev,
                    coriolisWindMil    = coriolisWind,
                    aerodynamicJumpMil = aerodynamicJumpMils(params.twistRateIn, params.twistDirection,
                                                             cw, params.muzzleVelocityFps,
                                                             params.bulletWeightGr, params.bulletDiameterIn,
                                                             params.bulletLengthIn),
                    isSubsonic         = isSubsonic,
                )
            }
        }
    }

    /**
     * Residual between observed DOPE elevation and the calm-weather physics model
     * at the conditions when that shot was fired.
     */
    internal fun entryBias(entry: DopeEntryInput, params: BallisticsParams): Double {
        val entryParams = params.copy(
            muzzleVelocityFps = entry.entryMuzzleVelocityFps ?: params.muzzleVelocityFps,
            tempF             = entry.tempF      ?: params.tempF,
            altitudeFt        = entry.altitudeFt ?: params.altitudeFt,
            pressureHpa       = entry.pressureHpa ?: params.pressureHpa,
            humidityPct       = entry.humidityPct ?: params.humidityPct,
            windSpeedFps      = 0.0,
            windAngleDeg      = 0.0,
        )
        val modelElev = calculateTrajectory(entryParams, listOf(entry.distanceYards))
            .getOrDefault(entry.distanceYards, null)?.elevationMil ?: 0.0
        return entry.adjustmentMil - modelElev
    }

    /**
     * Interpolate DOPE-corrected elevation and windage for [targetDistance].
     *
     * When ≥ 3 DOPE entries bracket or surround the target, uses monotone cubic
     * (Fritsch-Carlson) spline interpolation of the bias. Falls back to linear for
     * 2 entries, nearest-end constant for extrapolation outside the DOPE range.
     * Windage is purely physical — no bias correction applied.
     */
    fun interpolateDope(
        dopeEntries:    List<DopeEntryInput>,
        targetDistance: Double,
        params:         BallisticsParams,
    ): TrajectoryResult {
        val matchDay = calculateTrajectory(params, listOf(targetDistance))[targetDistance]
            ?: return TrajectoryResult(0.0, 0.0)
        if (dopeEntries.isEmpty()) return matchDay
        val biased = dopeEntries.sortedBy { it.distanceYards }
            .distinctBy { it.distanceYards }   // duplicate distances cause NaN in spline tangents
            .map { it to entryBias(it, params) }
        return applyBias(biased, matchDay, targetDistance)
    }

    /**
     * Precomputes entry biases for all [dopeEntries] against [params].
     * Separating this from [interpolateDopeAllFromBiases] lets callers cache the biases
     * across recalculations where only conditions, wind, or output distances changed —
     * none of which affect the bias residuals.
     *
     * Biases depend on: ammo (BC, MV), zero config, and each entry's stored conditions.
     * They are independent of: today's conditions, wind, output distances, and Coriolis.
     */
    fun computeEntryBiases(
        dopeEntries: List<DopeEntryInput>,
        params:      BallisticsParams,
    ): List<Pair<DopeEntryInput, Double>> =
        dopeEntries.sortedBy { it.distanceYards }
            .distinctBy { it.distanceYards }
            .map { it to entryBias(it, params) }

    /**
     * Batch DOPE interpolation using biases pre-computed by [computeEntryBiases].
     * The match-day trajectory pass still runs; only the per-entry RK4 bias integrations
     * are skipped. Use when the bias-affecting inputs (ammo, zero, entries) are unchanged.
     */
    fun interpolateDopeAllFromBiases(
        precomputedBiases: List<Pair<DopeEntryInput, Double>>,
        targetDistances:   List<Double>,
        params:            BallisticsParams,
    ): Map<Double, TrajectoryResult> {
        if (targetDistances.isEmpty()) return emptyMap()
        val matchDayAll = calculateTrajectory(params, targetDistances)
        if (precomputedBiases.isEmpty()) return matchDayAll
        return targetDistances.associateWith { dist ->
            applyBias(precomputedBiases, matchDayAll[dist] ?: TrajectoryResult(0.0, 0.0), dist)
        }
    }

    /**
     * Batch variant: computes DOPE biases once and evaluates all [targetDistances] in a
     * single trajectory pass. Use this when computing multiple distances to avoid
     * recomputing entry biases for every distance.
     */
    fun interpolateDopeAll(
        dopeEntries:     List<DopeEntryInput>,
        targetDistances: List<Double>,
        params:          BallisticsParams,
    ): Map<Double, TrajectoryResult> {
        if (targetDistances.isEmpty()) return emptyMap()
        val matchDayAll = calculateTrajectory(params, targetDistances)
        if (dopeEntries.isEmpty()) return matchDayAll
        val biased = dopeEntries.sortedBy { it.distanceYards }
            .distinctBy { it.distanceYards }   // duplicate distances cause NaN in spline tangents
            .map { it to entryBias(it, params) }
        return targetDistances.associateWith { dist ->
            applyBias(biased, matchDayAll[dist] ?: TrajectoryResult(0.0, 0.0), dist)
        }
    }

    internal fun applyBias(
        biased:         List<Pair<DopeEntryInput, Double>>,
        matchDay:       TrajectoryResult,
        targetDistance: Double,
    ): TrajectoryResult {
        // Exact match
        biased.firstOrNull { (e, _) -> e.distanceYards == targetDistance }
            ?.let { (_, bias) -> return matchDay.copy(elevationMil = matchDay.elevationMil + bias) }

        // Constant extrapolation outside the DOPE range
        if (targetDistance < biased.first().first.distanceYards)
            return matchDay.copy(elevationMil = matchDay.elevationMil + biased.first().second)
        if (targetDistance > biased.last().first.distanceYards)
            return matchDay.copy(elevationMil = matchDay.elevationMil + biased.last().second)

        // Monotone cubic spline across all entries
        val xs = biased.map { (e, _) -> e.distanceYards }.toDoubleArray()
        val ys = biased.map { (_, b) -> b }.toDoubleArray()
        return matchDay.copy(elevationMil = matchDay.elevationMil + MonotoneCubicSpline(xs, ys).eval(targetDistance))
    }

    // ── Unit conversion helpers ───────────────────────────────────────────────

    fun milToMoa(mil: Double): Double = mil * 3.43775
    fun mphToFps(mph: Double): Double = mph * 1.46667
}
