#!/usr/bin/env python3
"""Independent standard G7/G1 point-mass solver.

Purpose
-------
1. Cross-validate ScopeDOPE's BallisticsEngine.kt physics with a *separately
   implemented* solver (different integration scheme + different Cd interpolation),
   corroborating the JBM cross-check in BALLISTICS_ENGINE_ACCURACY_ANALYSIS.md.
2. Generate G7 reference trajectories for the Sub-3 test cases so the ShootersCalculator
   reference (which uses a lower effective BC) can be replaced with values that use the
   same standard McCoy G7 drag table the engine implements (== what JBM produces).

Independence
------------
The drag *table* values are the published McCoy/JBM standard (the same numbers JBM uses;
not engine-specific). Everything else — the RK4 step loop, the piecewise-LINEAR Cd
interpolation (engine uses monotone cubic Hermite), the drop/mil geometry — is implemented
here from first principles. Agreement with the engine therefore tests the engine, and the
output is a legitimate stand-in for JBM's standard-table trajectory.

Standard point-mass drag, distance domain:
    dv/ds = -KD * Cd(v_rel / a) * v_rel^2 / (bcEff * v)
    KD = rho0 / (2 * 144) = 0.0764 / 288   (lb/ft^3 -> consistent units)
"""
from __future__ import annotations
import math

KD = 0.0764 / 288.0
STD_PRESSURE_HPA = 1013.25
STANDARD_TEMP_F = 59.0

# ── Standard McCoy/JBM G7 drag table (Cd vs Mach) — matches BallisticsEngine.kt nodes ──
G7_MACH = [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]
G7_CD =   [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_MACH = [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]
G1_CD =   [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]


def cd_linear(mach: float, table_mach, table_cd) -> float:
    """Piecewise-LINEAR interpolation (deliberately different from engine's cubic)."""
    if mach <= table_mach[0]:
        return table_cd[0]
    if mach >= table_mach[-1]:
        return table_cd[-1]
    lo, hi = 0, len(table_mach) - 1
    while lo + 1 < hi:
        mid = (lo + hi) // 2
        if table_mach[mid] <= mach:
            lo = mid
        else:
            hi = mid
    t = (mach - table_mach[lo]) / (table_mach[hi] - table_mach[lo])
    return table_cd[lo] + t * (table_cd[hi] - table_cd[lo])


def speed_of_sound_fps(temp_f: float) -> float:
    return 1116.45 * math.sqrt((temp_f + 459.67) / 518.67)


def air_density_ratio(alt_ft: float, temp_f: float, press_hpa: float = 0.0,
                      humidity_pct: float = 0.0) -> float:
    temp_k = (temp_f + 459.67) * 5.0 / 9.0
    std_temp_k = (STANDARD_TEMP_F + 459.67) * 5.0 / 9.0
    p = press_hpa if press_hpa > 0 else STD_PRESSURE_HPA * (1.0 - 6.87559e-6 * alt_ft) ** 5.2561
    temp_c = (temp_f - 32.0) * 5.0 / 9.0
    e_sat = 6.1078 * math.exp(17.27 * temp_c / (temp_c + 237.3))
    e_w = (humidity_pct / 100.0) * e_sat
    humidity_factor = 1.0 - 0.378 * e_w / max(p, 1.0)
    return (p / STD_PRESSURE_HPA) * (std_temp_k / temp_k) * humidity_factor


def bc_eff_for_velocity(v_rel, single_bc, zones, density_ratio):
    """Velocity-banded BC (mirrors the engine's new multi-zone path)."""
    valid = [z for z in zones if z[1] > 0 and z[0] >= 0]
    if len(valid) < 2:
        return single_bc / density_ratio
    valid = sorted(valid, key=lambda z: z[0], reverse=True)  # descending floor
    for floor, bc in valid:
        if v_rel >= floor:
            return bc / density_ratio
    return valid[-1][1] / density_ratio


def integrate(mv, single_bc, zones, table_mach, table_cd, dist_yd, sos,
              density_ratio, headwind=0.0, gravity=32.174, step_yd=0.5):
    """Independent RK4 spatial integrator. Returns (drop_ft, tof_s, term_vel)."""
    v = mv
    step_ft = step_yd * 3.0
    drop_ft = 0.0
    drop_vel = 0.0
    t_flight = 0.0
    total_yd = 0.0

    def dvds(vel):
        vr = max(vel + headwind, 1.0)
        bc_eff = bc_eff_for_velocity(vr, single_bc, zones, density_ratio)
        return -KD * cd_linear(vr / sos, table_mach, table_cd) * vr * vr / (bc_eff * max(vel, 1.0))

    while total_yd < dist_yd:
        k1 = dvds(v)
        k2 = dvds(max(v + 0.5 * step_ft * k1, 1.0))
        k3 = dvds(max(v + 0.5 * step_ft * k2, 1.0))
        k4 = dvds(max(v + step_ft * k3, 1.0))
        v_new = max(v + step_ft * (k1 + 2 * k2 + 2 * k3 + k4) / 6.0, 1.0)
        dt = step_ft / ((v + v_new) * 0.5)
        drop_ft += drop_vel * dt + 0.5 * gravity * dt * dt
        drop_vel += gravity * dt
        t_flight += dt
        v = v_new
        total_yd += step_yd
    return drop_ft, t_flight, v


def trajectory_mil(drop_d, drop_z, dist_yd, zero_yd, sight_height_in):
    if dist_yd <= 0 or zero_yd <= 0:
        return 0.0
    h_ft = sight_height_in / 12.0
    scale = dist_yd / zero_yd
    needed_ft = drop_d + h_ft - scale * (drop_z + h_ft)
    return (needed_ft * 12.0) / (dist_yd * 36.0 / 1000.0)


def trajectory(bc_model, bc, mv, zero_yd, sight_height_in, dists_yd,
               alt_ft=0.0, temp_f=59.0, press_hpa=0.0, humidity_pct=0.0,
               zones=None, step_yd=0.5):
    zones = zones or []
    table_mach, table_cd = (G7_MACH, G7_CD) if bc_model == "G7" else (G1_MACH, G1_CD)
    dr = air_density_ratio(alt_ft, temp_f, press_hpa, humidity_pct)
    sos = speed_of_sound_fps(temp_f)
    drop_z, _, _ = integrate(mv, bc, zones, table_mach, table_cd, zero_yd, sos, dr, step_yd=step_yd)
    out = {}
    for d in dists_yd:
        if d == 0:
            out[d] = 0.0
            continue
        drop_d, _, _ = integrate(mv, bc, zones, table_mach, table_cd, d, sos, dr, step_yd=step_yd)
        out[d] = trajectory_mil(drop_d, drop_z, d, zero_yd, sight_height_in)
    return out


if __name__ == "__main__":
    # Self-check: 6.5 CM Hornady 140gr ELD-M, G7 0.326, 2674 fps, 100yd zero, sea level/59F.
    # Engine (per analysis cross-check) gives 0.528 mil @200, 8.247 mil @1000.
    t = trajectory("G7", 0.326, 2674, 100, 1.5, [200, 1000])
    print(f"65cm Hornady ELD-M  200yd={t[200]:.3f} mil (engine 0.528)  1000yd={t[1000]:.3f} mil (engine 8.247)")
