OpenBallistics: MCP Debugger Infrastructure for Parity Debugging

debug_py_step.py

#!/usr/bin/env python3
from py_ballisticcalc import (
    Ammo, Angular, Atmo, Calculator, Distance, DragModel,
    Pressure, Shot, TableG7, Temperature, Velocity, Weapon, Weight, Wind,
)

bc = 0.209
dm = DragModel(bc=bc, drag_table=TableG7,
               weight=Weight.Grain(77.0),
               diameter=Distance.Millimeter(5.56),
               length=Distance.Millimeter(22.1))
ammo = Ammo(dm=dm, mv=Velocity.MPS(752.1))
weapon = Weapon(sight_height=Distance.Millimeter(90.0),
                twist=Distance.Inch(9.0))

zero_shot = Shot(weapon=weapon, ammo=ammo, atmo=Atmo.icao())
calc = Calculator()
calc.set_weapon_zero(zero_shot, Distance.Meter(100.0))

current_atmo = Atmo(
    altitude=Distance.Meter(1363),
    pressure=Pressure.hPa(934.3),
    temperature=Temperature.Celsius(-7.3),
    humidity=71.7,
)

wind_dir = Angular.Degree(float(((4 - 6) * 30) % 360))
sustained_wind = Wind(velocity=Velocity.MPS(0.3), direction_from=wind_dir)

shot = Shot(
    weapon=weapon, ammo=ammo, atmo=current_atmo,
    winds=[sustained_wind],
    look_angle=Angular.Degree(3.6),
    cant_angle=Angular.Degree(0.8),
    latitude=-47.8, azimuth=83.5,
)

result = calc.fire(
    shot, trajectory_range=Distance.Meter(800),
    trajectory_step=Distance.Meter(50), raise_range_error=False,
)
print("DONE")

debugger-findings.md

Debugger-based parity investigation

Debugger setup

Infrastructure

Two mcp-debugger instances run simultaneously:

• Port 3001 -- Docker container (docker.io/debugmcp/mcp-debugger:latest), --network host, workspace mounted at /workspace. Provides Python, JavaScript, and Java/JDI adapters.
• Port 3002 -- Host npm package (@debugmcp/mcp-debugger). Provides Python adapter using the host's Python environment (where py_ballisticcalc is installed from .ai/sources/).

Java adapter requires the Docker image because the npm package ships without a compiled JDI bridge (JdiDapServer.class). Building from source (pnpm --filter @debugmcp/adapter-java run build:adapter) also works but Docker is simpler.

Both use Streamable HTTP transport. Communication is via curl:

# Helper script: .ai/scripts/mcp.sh <port> <session_id> <json_rpc_body>
curl -s -X POST "http://localhost:${PORT}/mcp" \
  -H "Content-Type: application/json" \
  -H "Accept: application/json, text/event-stream" \
  -H "Mcp-Session-Id: $SESSION_ID" \
  -d "$BODY" | grep '^data:' | sed 's/^data: //'

Python side

Script: .ai/scripts/debug_py_step.py -- calls calc.fire() with hardcoded fixture inputs (seed=42, fixture 0: .223 Rem, G7, BC=0.209). This exercises the full production path: set_weapon_zero (zero-finding) then engine.integrate (trajectory).

Breakpoints use conditions to distinguish zero-finding from trajectory:

• filter_flags == 0 -- zero-finding calls
• filter_flags != 0 -- trajectory calls
• integration_step_count == N -- specific step number

Kotlin/JVM side

Class: core/src/jvmTest/kotlin/.../DebugOneStep.kt -- constructs BallisticInput with the same fixture parameters, calls solver.findZeroAngle() then solver.computeTrajectory().

Launched with JDWP:

noglob java \
  "-agentlib:jdwp=transport=dt_socket,server=y,suspend=y,address=*:5005" \
  -cp "$CP" org.openballistics.engine.DebugOneStep

mcp-debugger attaches via create_debug_session with host=localhost, port=5005.

Breakpoints use FQCN org.openballistics.engine.TrajectorySolver and line numbers. To separate zero-finding from trajectory, breakpoints are placed at lines unique to each code path:

• Line 271 (simulateTrajectory) -- zero-finding
• Line 142 (integrate) -- trajectory

Experiment 1: initialization comparison (trajectory phase)

Breakpoints: Python at rk4.py:185 (after initial state setup, trajectory phase), Kotlin at TrajectorySolver.kt:136 (same logical point).

Results

Value	Python (metric)	Kotlin (metric)	Delta
pos.x	0	0	0
pos.y	-0.08999122716083738	-0.08999122716083738	0
pos.z	-0.00125659623052308	-0.00125659623052307	~1e-17
vel.x	750.528517415018	750.528513793092	3.6e-06
vel.y	48.593784592316	48.593822839031	-3.8e-05
vel.z	0.019115455799	0.019115990798	-5.3e-07
barrel_elev	0.064655840669	0.064655891727	-5.1e-08
barrel_az	0.000025469326	0.000025470039	-7.1e-10
cant_cos	0.999902524009	0.999902524009	0
cant_sin	0.013962180339	0.013962180339	~1e-17
look/slope_rad	0.062831853072	0.062831853072	~1e-16
sight_height	0.09	0.09	0
bc	0.209	0.209	0

Conclusion: Velocity diverges because barrel elevation differs by 5.1e-08 rad. barrel_elevation = slope + cos(cant) * zeroAngle. Since slope and cant match, the zero angles differ by 5.1e-08 rad.

Follow-up: using py's exact zero angle

Changed Kotlin test to use fixture.zero_angle_rad (with full precision, not rounded) instead of solver.findZeroAngle().

Results: errors dropped from ~1e-2 cm to ~1e-6 cm. Only drop still failed (velocity and windage passed). Remaining drop error at 2200m was 3.2e-05 cm (31.8x over 1e-6 tolerance).

Conclusion: zero angle difference is the primary error source (~99.9%). There is a residual per-step drop error that compounds over distance.

Important: zero_angle_rad is an output, not an input. The Kotlin engine must compute its own zero angle and match py's value. The experiment was diagnostic only to separate zero-angle error from RK4-step error.

Experiment 2: zero-finding RK4 step comparison

To eliminate the zero-angle difference as a confound, compared one RK4 step during zero-finding (both engines start with barrelElevation=0, no wind, no coriolis, ICAO atmosphere).

Breakpoints: Python at rk4.py:235 (before velocity update, step 2), Kotlin at TrajectorySolver.kt:314 (return statement, step 2).

Pre-step state (t=0.0025, after one completed step)

Value	Python (metric)	Kotlin (metric)	Delta
x	1.87859722387	1.87859722102	2.9e-09
y	-0.09003062783	-0.09003062783	~3e-17
vx	750.778551573	750.778551433	1.1e-06
vy	-0.024495088099	-0.024495088099	~3e-14
density_ratio	1.0004255875331	1.0004255875331	0
km	~equivalent	~equivalent	~1e-15

Key observation: vx diverges by 1.1e-06 but vy matches to 1e-14. Position x diverges by 2.9e-09 (accumulated from vx), position y matches perfectly. This pattern -- x diverges, y does not -- is unexpected because both components use the same formula.

Acceleration comparison (k1-k4)

Value	Python (metric)	Kotlin (metric)	Delta
a1.vx	-528.0651954	-528.0651938	-1.6e-06
a1.vy	-9.78942165338	-9.78942165340	2.5e-11
a2.vx	-527.1370583	-527.1370567	-1.6e-06
a2.vy	-9.78083754764	-9.78083754768	3.8e-11
a3.vx	-527.1386889	-527.1386873	-1.6e-06
a3.vy	-9.78084504821	-9.78084504824	3.7e-11
a4.vx	-526.2129903	-526.2129888	-1.6e-06
a4.vy	-9.77228352435	-9.77228352440	4.9e-11

The x-acceleration diverges by ~1.6e-06 consistently across k1-k4. The y-acceleration matches to ~1e-11.

Ruled-out culprits

1. Unit system (imperial vs metric) -- RULED OUT

Evidence: scripted analysis showed max floating-point difference from unit conversion is ~1e-14. Imperial RK4 rewrite proved the unit system is not the cause (made things 1e+8 worse). The Mach number computed via Rankine vs Kelvin constants matches to 1.5e-14.

2. Drag constant (DRAG_CONSTANT) -- RULED OUT

Evidence: Kotlin DRAG_CONSTANT = 0.3048 / 2.08551e-04 = 1461.513.... Python 2.08551e-04. Both use the same underlying value. The km values match to ~1e-15 when unit-converted.

3. Gravity constant -- RULED OUT

Evidence: Python -32.17405 ft/s^2. Kotlin 32.17405 * 0.3048 = 9.80665 m/s^2. Both standard g. Acceleration in y matches to ~1e-11, confirming gravity is identical.

4. PCHIP drag table interpolation -- RULED OUT

Evidence: cd_from_table values match to 1e-13 (from earlier instrumentation experiment, dimensionless values).

5. Coriolis constants -- RULED OUT

Evidence: Both use 7.2921159e-5 rad/s. sinLat, cosLat, sinAz, cosAz all match. Coriolis is disabled during zero-finding anyway.

6. Speed of sound formula (Rankine vs Kelvin constants) -- RULED OUT

Evidence: 49.0223 * sqrt(R) (imperial) and 20.0467 * sqrt(K) (metric) differ by ~4.6 ppm. However, this affects BOTH Python and Kotlin equally (both have the same base/altitude inconsistency). Mach number computed via either path matches to 1.5e-14 when properly unit-converted. Changing Kotlin to use 20.0467 * sqrt(K) made things WORSE (1.6e-03 m/s error) because Python uses the Rankine formula for base speed of sound.

7. Wind decomposition -- RULED OUT

Evidence: During zero-finding, wind is zero. For trajectory, the decomposition was verified: clock 4 -> headwind=-0.15 m/s, crosswind=0.2598 m/s in both engines.

8. Spin drift formula -- RULED OUT

Evidence: py /12 (inches to feet), kt *0.0254 (inches to meters). Both convert from inches, producing equivalent physical results.

9. Air density (CIPM-2007) -- RULED OUT

Evidence: density_ratio matches exactly (0.0 delta) between engines, confirming both compute identical air density.

Key finding: py_ballisticcalc uses 3.2808399 as ft/m constant

Python's Velocity unit uses 3.2808399 for fps-to-mps conversion, NOT the exact 1/0.3048 = 3.2808398950131.... Delta: 4.99e-09.

Evidence:

>>> Velocity.MPS(1.0) >> Velocity.FPS
3.2808399
>>> 1/0.3048
3.2808398950131

This means v0_fps = 752.1 * 3.2808399 = 2467.519688790 in Python, while 752.1 / 0.3048 = 2467.519685039. Delta: 3.75e-06 fps.

This 3.75e-06 fps offset in v0 propagates through:

• The drag computation k_m * v * |v| → 1.6e-06 acceleration offset
• Zero-finding convergence → 5.1e-08 zero angle difference
• Every RK4 step → cumulative position/velocity drift

Impact on DRAG_CONSTANT

The DRAG_CONSTANT converts py's imperial drag constant to metric. Original: 0.3048 / 2.08551e-04 = 1461.5130112059 (uses exact 0.3048). Should be: 1 / (2.08551e-04 * 3.2808399) = 1461.5130089843922 (uses py's constant).

Fix applied but test results unchanged (1047 failures). The DRAG_CONSTANT affects km by 2.2e-06 relative, but the velocity offset through v0_fps affects the entire initial state setup.

Why ax diverges but not ay

With barrelElevation=0: vx = v0 (large), vy = 0. The drag force is km * vAir * |vAir| where vAir ≈ (v0, 0, 0). So ax ≈ km*v0^2 (sensitive to v0 errors) while ay ≈ 0 - G (insensitive, dominated by gravity). Oracle consultation confirmed this analysis.

Root cause chain

1. py_ballisticcalc uses 3.2808399 (not 1/0.3048) as its ft/m
2. ALL imperial values in py carry a ~5e-9 relative offset
3. Kotlin must either: a. Use the same constant everywhere, OR b. Work entirely in imperial like py does
4. Currently: Atmosphere.kt uses 3.2808399 but TrajectorySolver.kt uses 3.2808398950131 -- inconsistent

What was attempted

1. Changed DRAG_CONSTANT from 1461.5130112059 (0.3048-based) to 1461.5130089843922 (3.2808399-based). Result: test failures unchanged (1047). The constant alone is not enough.

2. Changed M_TO_FT from 3.2808398950131 to 3.2808399. Result: made things worse (same failure count but different delta profile). The M_TO_FT constant is used for altitude/pressure computations where the exact value matters.

3. Changed gravity from 32.17405 * 0.3048 to 32.17405 / 3.2808399. Result: didn't compile as const val; when hardcoded, negligible impact since gravity affects y-component which already matches.

All changes reverted. Working tree clean.

Why individual constant fixes don't work

The problem is systemic: Python does ALL arithmetic in imperial using 3.2808399 for conversions. The intermediate results (v_fps, mach_fps, k_m, acceleration_fps2) all carry the ~5e-9 relative bias from this constant. These biases are self-consistent within the imperial system.

Kotlin works in metric. It converts to feet only for atmospheric lookups. The metric intermediate values are computed with exact constants (0.3048, 1/0.3048). There's no single constant to fix -- the issue is the unit system boundary itself.

Next steps

The clearest path to parity:

1. Step 0 verification: Capture the complete state at t=0 (before ANY RK4 step) from both debuggers and verify they're identical. If they match exactly, the per-step error is ~1.6e-06/528 ≈ 3e-9 relative, which over 1000+ steps compounds to the observed 1e-6 to 1e-2 cm errors.

2. Candidate approach: Make the Kotlin engine's drag computation replicate Python's arithmetic exactly. Compute drag in imperial (using py's constants) then convert the result to metric. This isolates the conversion to one place rather than having it spread across multiple constants.

   Concretely: inside rk4Step, compute km as densityRatio * cd * 2.08551e-04 / bc (same as py), use it with velocities converted to fps (v * 3.2808399), compute acceleration in fps^2, then convert back to m/s^2 (/ 3.2808399).

3. Alternative: Accept that metric vs imperial will always produce ~1e-9 relative per-step error, and tighten the zero-finding algorithm to converge to the same value as py's (which absorbs the systematic drag error into the zero angle). This might achieve 1e-6 parity without rewriting the drag computation.

Fixes applied (committed)

Fix 1: PCHIP interpolation in simulateTrajectory

Kotlin's zero-finding used linear interpolation at the target distance. Python's _integrate uses PCHIP (3-point cubic Hermite) interpolation. This caused the zero angle to converge to a different value by ~5.1e-08 radians.

Fix: Added prevPrevState tracking in simulateTrajectory and call pchipInterpolateAtDistance when 3 points are available.

Result: barrel elevation difference dropped from 5.1e-08 to 1.5e-12 radians. Test failures: 1027 -> 253.

Fix 2: Imperial drag constants

Kotlin's DRAG_CONSTANT used 0.3048 / 2.08551e-04. Python uses 2.08551e-04 directly with velocities in fps (using 3.2808399 as ft/m, not exact 1/0.3048). Gravity used 32.17405 * 0.3048 but should be 32.17405 / 3.2808399 to match py's arithmetic.

Fix: Compute km = density * cd * 2.08551e-04 * 3.2808399 / bc (matching py's conversion path). Gravity: G = 32.17405 / 3.2808399.

Result: Marginal improvement (~1.5 ppb per step reduction).

Current state: 253/1800 failures

All failures are drop-only, scaling with distance:

• Max delta: 4.6e-5 cm at 2200m (.338 Lapua Mag)
• Typical: 1-5e-6 cm at 600-800m

Conversion factor analysis

py_ballisticcalc uses different conversion constants for different unit types:

Unit type	py constant	Kotlin constant	Match?
Velocity (fps/mps)	3.2808399	N/A (metric)	N/A
Distance (ft/m)	39.37007874 in/m (exact)	3.2808398950131	exact
SOS (fps)	49.0223 * sqrt(R)	49.0223 * sqrt(R) / 3.2808399	yes
SOS (m/s) altitude	20.0467 * sqrt(K)	20.0467 * sqrt(K)	yes
Drag magic	2.08551e-04	2.08551e-04 * 3.2808399	yes
Gravity	32.17405 fps^2	32.17405 / 3.2808399 m/s^2	yes
Coriolis	omega * v_fps	omega * v_mps	equivalent

Outstanding: per-step y-position error ~1.3e-10

Debugger comparison at step 2 of trajectory (after PCHIP fix):

Value	py (metric)	kt (metric)	delta
x	1.874704473715650e+00	1.874704473715473e+00	1.8e-13
y	3.135813158212049e-02	3.135813144808181e-02	1.3e-10
z	-1.209663359953115e-03	-1.209663361824947e-03	1.9e-12
vx	7.492358061240687e+02	7.492358061239977e+02	7.1e-11
vy	4.848575743791208e+01	4.848575743901033e+01	-1.1e-09
vz	1.843123267766896e-02	1.843123269293303e-02	-1.5e-11

The y-position has the largest delta (1.3e-10), which compounds over ~1000 steps to the observed 1e-5 to 5e-5 cm errors.

Root cause: TBD. Need to compare k1-k4 within this step to find which substep introduces the y error.

Debugger setup preserved

Both mcp-debugger instances are still running:

• Port 3001: Docker (Java/JDI)
• Port 3002: Host npm (Python/debugpy)

Scripts in .ai/scripts/:

• mcp.sh -- generic MCP HTTP call helper
• debug_py_step.py -- Python debug target (calls calc.fire())
• dual_debug.py -- dual-debugger comparison driver

JVM debug target: DebugOneStep.kt (untracked, in jvmTest)

dual_debug.py

#!/usr/bin/env python3
import json
import subprocess
import sys
import os

FT2M = 0.3048
FPS2MPS = 0.3048

class McpDebugger:
    def __init__(self, port, mcp_session_id, dbg_session_id):
        self.port = port
        self.mcp_sid = mcp_session_id
        self.dbg_sid = dbg_session_id
        self._id = 1000

    def _call(self, method, params):
        self._id += 1
        body = json.dumps({
            "jsonrpc": "2.0", "id": self._id,
            "method": method, "params": params
        })
        script = os.path.join(os.path.dirname(__file__), "mcp.sh")
        r = subprocess.run(
            [script, str(self.port), self.mcp_sid, body],
            capture_output=True, text=True, timeout=30
        )
        if not r.stdout.strip():
            return None
        return json.loads(r.stdout.strip())

    def tool(self, name, args):
        result = self._call("tools/call", {"name": name, "arguments": args})
        if result and "result" in result:
            return json.loads(result["result"]["content"][0]["text"])
        return None

    def eval(self, expr):
        r = self.tool("evaluate_expression", {
            "sessionId": self.dbg_sid, "expression": expr
        })
        if r and r.get("success"):
            return r["result"]
        return f"ERR:{r}"

    def eval_float(self, expr):
        v = self.eval(expr)
        try:
            return float(v)
        except (ValueError, TypeError):
            return None

    def step_over(self):
        self.tool("step_over", {"sessionId": self.dbg_sid})

    def step_into(self):
        self.tool("step_into", {"sessionId": self.dbg_sid})

    def continue_exec(self):
        self.tool("continue_execution", {"sessionId": self.dbg_sid})

    def get_line(self):
        r = self.tool("get_stack_trace", {"sessionId": self.dbg_sid})
        if r and r.get("stackFrames"):
            f = r["stackFrames"][0]
            return f["name"], f["line"]
        return "?", -1

    def get_locals(self):
        r = self.tool("get_local_variables", {"sessionId": self.dbg_sid})
        if r and r.get("variables"):
            return {v["name"]: v["value"] for v in r["variables"]}
        return {}


def compare(label, py_val, kt_val, convert_py=True):
    if py_val is None or kt_val is None:
        print(f"  {label:25s}  py={py_val}  kt={kt_val}  MISSING")
        return
    if convert_py:
        py_metric = py_val * FT2M
    else:
        py_metric = py_val
    delta = abs(py_metric - kt_val)
    ok = "OK" if delta < 1e-14 else f"DELTA={delta:.6e}"
    if delta >= 1e-14:
        ok = f"*** DELTA={delta:.6e} ***"
    print(f"  {label:25s}  py={py_metric:+.15e}  kt={kt_val:+.15e}  {ok}")


def compare_raw(label, py_val, kt_val):
    if py_val is None or kt_val is None:
        print(f"  {label:25s}  py={py_val}  kt={kt_val}  MISSING")
        return
    delta = abs(py_val - kt_val)
    ok = "OK" if delta < 1e-14 else f"*** DELTA={delta:.6e} ***"
    print(f"  {label:25s}  py={py_val:+.15e}  kt={kt_val:+.15e}  {ok}")


if __name__ == "__main__":
    py_port = int(sys.argv[1])
    py_mcp = sys.argv[2]
    py_dbg = sys.argv[3]
    kt_port = int(sys.argv[4])
    kt_mcp = sys.argv[5]
    kt_dbg = sys.argv[6]

    py = McpDebugger(py_port, py_mcp, py_dbg)
    kt = McpDebugger(kt_port, kt_mcp, kt_dbg)

    print("=== Verifying breakpoint positions ===")
    print(f"  PY: {py.get_line()}")
    print(f"  KT: {kt.get_line()}")

    import math

    print("\n=== INITIALIZATION ===")
    print("--- Barrel elevation & azimuth (radians) ---")
    py_be = py.eval_float("props.barrel_elevation_rad")
    kt_be = kt.eval_float("barrelElevation")
    compare_raw("barrel_elevation_rad", py_be, kt_be)

    py_ba = py.eval_float("props.barrel_azimuth_rad")
    kt_ba = kt.eval_float("barrelAzimuth")
    compare_raw("barrel_azimuth_rad", py_ba, kt_ba)

    print("--- Muzzle velocity ---")
    py_v0 = py.eval_float("props.muzzle_velocity_fps")
    kt_v0_mps = kt.eval_float("state.getVx() / Math.cos(barrelElevation) / Math.cos(barrelAzimuth)")
    if py_v0 and kt_v0_mps:
        compare("v0", py_v0, kt_v0_mps)

    print("--- Initial position (ft->m) ---")
    py_rx = py.eval_float("range_vector.x")
    py_ry = py.eval_float("range_vector.y")
    py_rz = py.eval_float("range_vector.z")
    kt_rx = kt.eval_float("state.getX()")
    kt_ry = kt.eval_float("state.getY()")
    kt_rz = kt.eval_float("state.getZ()")
    compare("pos.x", py_rx, kt_rx)
    compare("pos.y", py_ry, kt_ry)
    compare("pos.z", py_rz, kt_rz)

    print("--- Initial velocity (fps->m/s) ---")
    py_vx = py.eval_float("velocity_vector.x")
    py_vy = py.eval_float("velocity_vector.y")
    py_vz = py.eval_float("velocity_vector.z")
    kt_vx = kt.eval_float("state.getVx()")
    kt_vy = kt.eval_float("state.getVy()")
    kt_vz = kt.eval_float("state.getVz()")
    compare("vel.x", py_vx, kt_vx)
    compare("vel.y", py_vy, kt_vy)
    compare("vel.z", py_vz, kt_vz)

    print("--- Dimensionless params ---")
    py_bc = py.eval_float("props.bc")
    compare_raw("bc", py_bc, 0.209)

    py_sg = py.eval_float("props.stability_coefficient")
    print(f"  {'stability':25s}  py={py_sg}")

    py_cant_cos = py.eval_float("props.cant_cosine")
    py_cant_sin = py.eval_float("props.cant_sine")
    kt_cant = kt.eval_float("this.input.getCant().getRadians()")
    if kt_cant is not None:
        compare_raw("cant_cos", py_cant_cos, math.cos(kt_cant))
        compare_raw("cant_sin", py_cant_sin, math.sin(kt_cant))

    py_look = py.eval_float("props.look_angle_rad")
    kt_slope = kt.eval_float("this.input.getSlope().getRadians()")
    compare_raw("look/slope_rad", py_look, kt_slope)

    py_sh = py.eval_float("props.sight_height_ft")
    kt_sh = kt.eval_float("this.input.getSightHeight().getMeters()")
    if py_sh and kt_sh:
        compare("sight_height", py_sh, kt_sh)

    print("--- Wind vector (fps->m/s) ---")
    py_wx = py.eval_float("wind_vector.x")
    py_wy = py.eval_float("wind_vector.y")
    py_wz = py.eval_float("wind_vector.z")
    print(f"  {'wind.x (py fps)':25s}  {py_wx}")
    print(f"  {'wind.y (py fps)':25s}  {py_wy}")
    print(f"  {'wind.z (py fps)':25s}  {py_wz}")

    print("\n=== Done initialization check ===")

dump_py_step.py

#!/usr/bin/env python3
"""Dump all intermediate values from ONE RK4 step in py_ballisticcalc."""
import math
import json

from py_ballisticcalc import (
    Ammo, Angular, Atmo, Calculator, Distance, DragModel,
    Pressure, Shot, TableG1, TableG7, Temperature, Velocity, Weapon, Weight, Wind,
)
from py_ballisticcalc.shot import ShotProps
from py_ballisticcalc.conditions import Coriolis
from py_ballisticcalc.engines.rk4 import RK4IntegrationEngine
from py_ballisticcalc.vector import Vector

# --- Fixture inputs (seed=42, fixture 0) ---
bc = 0.209
drag_model_str = "G7"
weight_grains = 77.0
diameter_mm = 5.56
length_mm = 22.1
v0_mps = 752.1
zero_distance_m = 100.0
sight_height_mm = 90.0
twist_inches = 9.0
temperature_c = -7.3
pressure_hpa = 934.3
humidity_pct = 71.7
altitude_m = 1363
wind_clock = 4
wind_speed_mps = 0.3
slope_deg = 3.6
cant_deg = 0.8
latitude_deg = -47.8
azimuth_deg = 83.5

# --- Build objects ---
drag_table = TableG7
dm = DragModel(bc=bc, drag_table=drag_table,
               weight=Weight.Grain(weight_grains),
               diameter=Distance.Millimeter(diameter_mm),
               length=Distance.Millimeter(length_mm))
ammo = Ammo(dm=dm, mv=Velocity.MPS(v0_mps))
weapon = Weapon(sight_height=Distance.Millimeter(sight_height_mm),
                twist=Distance.Inch(twist_inches))

# Zero the weapon
zero_shot = Shot(weapon=weapon, ammo=ammo, atmo=Atmo.icao())
calc = Calculator()
calc.set_weapon_zero(zero_shot, Distance.Meter(zero_distance_m))
zero_angle_rad = float(weapon.zero_elevation >> Angular.Radian)

# Current atmosphere
current_atmo = Atmo(
    altitude=Distance.Meter(altitude_m),
    pressure=Pressure.hPa(pressure_hpa),
    temperature=Temperature.Celsius(temperature_c),
    humidity=humidity_pct,
)

def clock_to_degrees(clock):
    return float(((clock - 6) * 30) % 360)

wind_dir = Angular.Degree(clock_to_degrees(wind_clock))
sustained_wind = Wind(velocity=Velocity.MPS(wind_speed_mps), direction_from=wind_dir)

shot = Shot(
    weapon=weapon, ammo=ammo, atmo=current_atmo,
    winds=[sustained_wind],
    look_angle=Angular.Degree(slope_deg),
    cant_angle=Angular.Degree(cant_deg),
    latitude=latitude_deg, azimuth=azimuth_deg,
)

# --- Build ShotProps (what the engine uses internally) ---
props = ShotProps.from_shot(shot)

print("=== SHOT PROPS (internal imperial units) ===")
print(f"muzzle_velocity_fps: {props.muzzle_velocity_fps:.15e}")
print(f"sight_height_ft: {props.sight_height_ft:.15e}")
print(f"barrel_elevation_rad: {props.barrel_elevation_rad:.15e}")
print(f"barrel_azimuth_rad: {props.barrel_azimuth_rad:.15e}")
print(f"look_angle_rad: {props.look_angle_rad:.15e}")
print(f"cant_cosine: {props.cant_cosine:.15e}")
print(f"cant_sine: {props.cant_sine:.15e}")
print(f"bc: {props.bc:.15e}")
print(f"alt0_ft: {props.alt0_ft:.15e}")
print(f"stability_coefficient: {props.stability_coefficient:.15e}")
engine = RK4IntegrationEngine(None)
props.calc_step = engine.get_calc_step()
print(f"calc_step: {props.calc_step:.15e}")
print(f"twist_inch: {props.twist_inch:.15e}")

# --- Initial state ---
velocity = props.muzzle_velocity_fps
range_vector = Vector(0.0, -props.cant_cosine * props.sight_height_ft,
                       -props.cant_sine * props.sight_height_ft)
velocity_vector = Vector(
    math.cos(props.barrel_elevation_rad) * math.cos(props.barrel_azimuth_rad),
    math.sin(props.barrel_elevation_rad),
    math.cos(props.barrel_elevation_rad) * math.sin(props.barrel_azimuth_rad),
) * velocity

print()
print("=== INITIAL STATE (imperial) ===")
print(f"range_vector: x={range_vector.x:.15e} y={range_vector.y:.15e} z={range_vector.z:.15e}")
print(f"velocity_vector: x={velocity_vector.x:.15e} y={velocity_vector.y:.15e} z={velocity_vector.z:.15e}")

# --- Wind ---
wind_sock_vector = sustained_wind.vector
print()
print("=== WIND VECTOR (imperial) ===")
print(f"wind: x={wind_sock_vector.x:.15e} y={wind_sock_vector.y:.15e} z={wind_sock_vector.z:.15e}")

# --- Step 1 intermediates ---
print()
print("=== RK4 STEP 1 INTERMEDIATES ===")

# Density and mach at initial altitude
density_ratio, mach_fps = props.get_density_and_mach_for_altitude(range_vector.y)
print(f"density_ratio: {density_ratio:.15e}")
print(f"mach_fps: {mach_fps:.15e}")

# Relative velocity
relative_velocity = velocity_vector - wind_sock_vector
relative_speed = relative_velocity.magnitude()
print(f"relative_velocity: x={relative_velocity.x:.15e} y={relative_velocity.y:.15e} z={relative_velocity.z:.15e}")
print(f"relative_speed: {relative_speed:.15e}")

# Mach number for drag lookup
mach_number = relative_speed / mach_fps
print(f"mach_number (for drag): {mach_number:.15e}")

# Drag
from py_ballisticcalc.interpolation import pchip_eval
cd = pchip_eval(props.drag_curve, mach_number)
drag_by_mach = cd * 2.08551e-04 / props.bc
print(f"cd_from_table: {cd:.15e}")
print(f"drag_by_mach (SDF): {drag_by_mach:.15e}")

k_m = density_ratio * drag_by_mach
print(f"k_m: {k_m:.15e}")

delta_time = props.calc_step
print(f"delta_time: {delta_time:.15e}")

# Gravity vector
gravity = Vector(0.0, -32.17405, 0.0)
print(f"gravity: x={gravity.x:.15e} y={gravity.y:.15e} z={gravity.z:.15e}")

# Coriolis
coriolis = props.coriolis
print(f"coriolis.sin_lat: {coriolis.sin_lat:.15e}")
print(f"coriolis.cos_lat: {coriolis.cos_lat:.15e}")
print(f"coriolis.sin_az: {coriolis.sin_az:.15e}")
print(f"coriolis.cos_az: {coriolis.cos_az:.15e}")
print(f"coriolis.range_east: {coriolis.range_east:.15e}")
print(f"coriolis.range_north: {coriolis.range_north:.15e}")
print(f"coriolis.cross_east: {coriolis.cross_east:.15e}")
print(f"coriolis.cross_north: {coriolis.cross_north:.15e}")

# --- k1 ---
v1 = velocity_vector
rel1 = v1 - wind_sock_vector

coriolis_a1 = coriolis.coriolis_acceleration_local(v1)
drag_a1 = -k_m * rel1 * rel1.magnitude()
a1 = gravity + coriolis_a1 + drag_a1

print()
print("--- k1 ---")
print(f"rel1: x={rel1.x:.15e} y={rel1.y:.15e} z={rel1.z:.15e}")
print(f"rel1.magnitude: {rel1.magnitude():.15e}")
print(f"coriolis_a1: x={coriolis_a1.x:.15e} y={coriolis_a1.y:.15e} z={coriolis_a1.z:.15e}")
print(f"drag_a1: x={drag_a1.x:.15e} y={drag_a1.y:.15e} z={drag_a1.z:.15e}")
print(f"a1: x={a1.x:.15e} y={a1.y:.15e} z={a1.z:.15e}")

# --- k2 ---
v2 = velocity_vector + 0.5 * delta_time * a1
rel2 = v2 - wind_sock_vector

coriolis_a2 = coriolis.coriolis_acceleration_local(v2)
drag_a2 = -k_m * rel2 * rel2.magnitude()
a2 = gravity + coriolis_a2 + drag_a2

print()
print("--- k2 ---")
print(f"v2: x={v2.x:.15e} y={v2.y:.15e} z={v2.z:.15e}")
print(f"rel2: x={rel2.x:.15e} y={rel2.y:.15e} z={rel2.z:.15e}")
print(f"coriolis_a2: x={coriolis_a2.x:.15e} y={coriolis_a2.y:.15e} z={coriolis_a2.z:.15e}")
print(f"drag_a2: x={drag_a2.x:.15e} y={drag_a2.y:.15e} z={drag_a2.z:.15e}")
print(f"a2: x={a2.x:.15e} y={a2.y:.15e} z={a2.z:.15e}")

# --- k3 ---
v3 = velocity_vector + 0.5 * delta_time * a2
rel3 = v3 - wind_sock_vector

coriolis_a3 = coriolis.coriolis_acceleration_local(v3)
drag_a3 = -k_m * rel3 * rel3.magnitude()
a3 = gravity + coriolis_a3 + drag_a3

print()
print("--- k3 ---")
print(f"v3: x={v3.x:.15e} y={v3.y:.15e} z={v3.z:.15e}")
print(f"rel3: x={rel3.x:.15e} y={rel3.y:.15e} z={rel3.z:.15e}")
print(f"coriolis_a3: x={coriolis_a3.x:.15e} y={coriolis_a3.y:.15e} z={coriolis_a3.z:.15e}")
print(f"drag_a3: x={drag_a3.x:.15e} y={drag_a3.y:.15e} z={drag_a3.z:.15e}")
print(f"a3: x={a3.x:.15e} y={a3.y:.15e} z={a3.z:.15e}")

# --- k4 ---
v4 = velocity_vector + delta_time * a3
rel4 = v4 - wind_sock_vector

coriolis_a4 = coriolis.coriolis_acceleration_local(v4)
drag_a4 = -k_m * rel4 * rel4.magnitude()
a4 = gravity + coriolis_a4 + drag_a4

print()
print("--- k4 ---")
print(f"v4: x={v4.x:.15e} y={v4.y:.15e} z={v4.z:.15e}")
print(f"rel4: x={rel4.x:.15e} y={rel4.y:.15e} z={rel4.z:.15e}")
print(f"coriolis_a4: x={coriolis_a4.x:.15e} y={coriolis_a4.y:.15e} z={coriolis_a4.z:.15e}")
print(f"drag_a4: x={drag_a4.x:.15e} y={drag_a4.y:.15e} z={drag_a4.z:.15e}")
print(f"a4: x={a4.x:.15e} y={a4.y:.15e} z={a4.z:.15e}")

# --- Final update ---
new_velocity = velocity_vector + (a1 + 2*a2 + 2*a3 + a4) * (delta_time / 6.0)
new_range = range_vector + (v1 + 2*v2 + 2*v3 + v4) * (delta_time / 6.0)

print()
print("=== AFTER STEP 1 ===")
print(f"new_range: x={new_range.x:.15e} y={new_range.y:.15e} z={new_range.z:.15e}")
print(f"new_velocity: x={new_velocity.x:.15e} y={new_velocity.y:.15e} z={new_velocity.z:.15e}")
print(f"new_speed: {new_velocity.magnitude():.15e}")

dump_zero_step.py

#!/usr/bin/env python3
"""Dump all intermediate values from ONE rk4Step in the ZERO-FINDING scenario.

Zero-finding conditions:
- ICAO standard atmosphere (15C, 1013.25 hPa, 0% humidity, sea level)
- No wind
- No Coriolis (lat=0, az=0)
- barrelElevation = 0 (first iteration)
- v0 = 752.1 m/s, sightHeight = 90mm, zeroDistance = 100m
"""
import math
import sys
import os

sys.path.insert(0, os.path.join(os.path.dirname(__file__), '../sources/py_ballisticcalc'))

from py_ballisticcalc import (
    Ammo, Angular, Atmo, Distance, DragModel,
    Shot, TableG7, Velocity, Weapon, Weight,
)
from py_ballisticcalc.shot import ShotProps
from py_ballisticcalc.engines.rk4 import RK4IntegrationEngine
from py_ballisticcalc.vector import Vector, ZERO_VECTOR
from py_ballisticcalc.interpolation import pchip_eval

# --- Inputs (zero-finding scenario) ---
bc = 0.209
v0_mps = 752.1
zero_distance_m = 100.0
sight_height_mm = 90.0
twist_inches = 9.0

# Build objects
dm = DragModel(bc=bc, drag_table=TableG7,
               weight=Weight.Grain(77.0),
               diameter=Distance.Millimeter(5.56),
               length=Distance.Millimeter(22.1))
ammo = Ammo(dm=dm, mv=Velocity.MPS(v0_mps))
weapon = Weapon(sight_height=Distance.Millimeter(sight_height_mm),
                twist=Distance.Inch(twist_inches))

# Zero shot: ICAO atmosphere, no wind, no coriolis
zero_shot = Shot(weapon=weapon, ammo=ammo, atmo=Atmo.icao())

# Build ShotProps exactly as set_weapon_zero does
props = ShotProps.from_shot(zero_shot)

# Force barrel_elevation = 0 (first iteration of zero-finding)
props.barrel_elevation_rad = 0.0
props.barrel_azimuth_rad = 0.0

engine = RK4IntegrationEngine(None)
props.calc_step = engine.get_calc_step()

FT_PER_M = 3.2808399  # Python's own constant

print("=== ZERO-FINDING SHOT PROPS ===")
print(f"muzzle_velocity_fps:    {props.muzzle_velocity_fps:.15e}")
print(f"sight_height_ft:        {props.sight_height_ft:.15e}")
print(f"barrel_elevation_rad:   {props.barrel_elevation_rad:.15e}")
print(f"barrel_azimuth_rad:     {props.barrel_azimuth_rad:.15e}")
print(f"bc:                     {props.bc:.15e}")
print(f"calc_step:              {props.calc_step:.15e}")
print(f"cant_cosine:            {props.cant_cosine:.15e}")
print(f"cant_sine:              {props.cant_sine:.15e}")

# Check coriolis
print(f"coriolis:               {props.coriolis}")

# Initial state (zero-finding uses cant=0, slope=0, barrelElev=0)
velocity = props.muzzle_velocity_fps
range_vector = Vector(0.0, -props.cant_cosine * props.sight_height_ft,
                      -props.cant_sine * props.sight_height_ft)
velocity_vector = Vector(
    math.cos(props.barrel_elevation_rad) * math.cos(props.barrel_azimuth_rad),
    math.sin(props.barrel_elevation_rad),
    math.cos(props.barrel_elevation_rad) * math.sin(props.barrel_azimuth_rad),
) * velocity

print()
print("=== INITIAL STATE (imperial) ===")
print(f"range_vector.x:         {range_vector.x:.15e}  ft")
print(f"range_vector.y:         {range_vector.y:.15e}  ft")
print(f"range_vector.z:         {range_vector.z:.15e}  ft")
print(f"velocity_vector.x:      {velocity_vector.x:.15e}  fps")
print(f"velocity_vector.y:      {velocity_vector.y:.15e}  fps")
print(f"velocity_vector.z:      {velocity_vector.z:.15e}  fps")

print()
print("=== INITIAL STATE (metric, /3.2808399) ===")
print(f"range_vector.x:         {range_vector.x/FT_PER_M:.15e}  m")
print(f"range_vector.y:         {range_vector.y/FT_PER_M:.15e}  m")
print(f"range_vector.z:         {range_vector.z/FT_PER_M:.15e}  m")
print(f"velocity_vector.x:      {velocity_vector.x/FT_PER_M:.15e}  m/s")
print(f"velocity_vector.y:      {velocity_vector.y/FT_PER_M:.15e}  m/s")
print(f"velocity_vector.z:      {velocity_vector.z/FT_PER_M:.15e}  m/s")

# Wind = zero (no wind in zero-finding)
wind_vector = Vector(0.0, 0.0, 0.0)
print()
print("=== WIND (imperial) ===")
print(f"wind_vector.x:          {wind_vector.x:.15e}  fps")
print(f"wind_vector.y:          {wind_vector.y:.15e}  fps")
print(f"wind_vector.z:          {wind_vector.z:.15e}  fps")

# Step 1 intermediates
print()
print("=== RK4 STEP 1 INTERMEDIATES ===")

# Density and mach at initial y (range_vector.y)
density_ratio, mach_fps = props.get_density_and_mach_for_altitude(range_vector.y)
print(f"density_ratio:          {density_ratio:.15e}")
print(f"mach_fps (speed of sound): {mach_fps:.15e}  fps")
print(f"mach_fps/FT_PER_M:      {mach_fps/FT_PER_M:.15e}  m/s")

# Relative velocity (no wind => same as ground velocity)
relative_velocity = velocity_vector - wind_vector
relative_speed = relative_velocity.magnitude()
print(f"relative_speed:         {relative_speed:.15e}  fps")
print(f"relative_speed/FT_PER_M:{relative_speed/FT_PER_M:.15e}  m/s")

# Mach number for drag
mach_number = relative_speed / mach_fps
print(f"mach_number:            {mach_number:.15e}")

# Drag
cd = pchip_eval(props.drag_curve, mach_number)
drag_by_mach = cd * 2.08551e-04 / props.bc
k_m = density_ratio * drag_by_mach
print(f"cd:                     {cd:.15e}")
print(f"drag_by_mach (SDF):     {drag_by_mach:.15e}  fps^-1")
print(f"k_m:                    {k_m:.15e}  ft^-1")
print(f"k_m * FT_PER_M:         {k_m * FT_PER_M:.15e}  m^-1  (= Kotlin km?)")

delta_time = props.calc_step
gravity = Vector(0.0, -32.17405, 0.0)  # Python's gravity in fps²
print(f"delta_time:             {delta_time:.15e}  s")
print(f"gravity.y:              {gravity.y:.15e}  fps²")
print(f"gravity.y/FT_PER_M:     {gravity.y/FT_PER_M:.15e}  m/s²")

# No coriolis in zero-finding
coriolis_fn = None

def acceleration(rel_vel: Vector, ground_vel: Vector) -> Vector:
    """Acceleration from gravity + drag (no coriolis in zero-finding)."""
    drag = -k_m * rel_vel * rel_vel.magnitude()
    return gravity + drag

# k1
v1 = velocity_vector
rel1 = v1 - wind_vector
a1 = acceleration(rel1, v1)

print()
print("--- k1 ---")
print(f"v1.x:                   {v1.x:.15e}  fps")
print(f"v1.y:                   {v1.y:.15e}  fps")
print(f"v1.z:                   {v1.z:.15e}  fps")
print(f"rel1.x:                 {rel1.x:.15e}  fps")
print(f"rel1.y:                 {rel1.y:.15e}  fps")
print(f"rel1.z:                 {rel1.z:.15e}  fps")
print(f"rel1.magnitude:         {rel1.magnitude():.15e}  fps")
drag1 = -k_m * rel1 * rel1.magnitude()
print(f"drag1.x:                {drag1.x:.15e}  fps²")
print(f"drag1.y:                {drag1.y:.15e}  fps²")
print(f"drag1.z:                {drag1.z:.15e}  fps²")
print(f"a1.x:                   {a1.x:.15e}  fps²")
print(f"a1.y:                   {a1.y:.15e}  fps²")
print(f"a1.z:                   {a1.z:.15e}  fps²")
print(f"a1.x/FT_PER_M:          {a1.x/FT_PER_M:.15e}  m/s²")
print(f"a1.y/FT_PER_M:          {a1.y/FT_PER_M:.15e}  m/s²")
print(f"a1.z/FT_PER_M:          {a1.z/FT_PER_M:.15e}  m/s²")

# k2
v2 = velocity_vector + 0.5 * delta_time * a1
rel2 = v2 - wind_vector
a2 = acceleration(rel2, v2)

print()
print("--- k2 ---")
print(f"v2.x:                   {v2.x:.15e}  fps")
print(f"v2.y:                   {v2.y:.15e}  fps")
print(f"v2.z:                   {v2.z:.15e}  fps")
print(f"rel2.x:                 {rel2.x:.15e}  fps")
print(f"a2.x:                   {a2.x:.15e}  fps²")
print(f"a2.y:                   {a2.y:.15e}  fps²")
print(f"a2.x/FT_PER_M:          {a2.x/FT_PER_M:.15e}  m/s²")
print(f"a2.y/FT_PER_M:          {a2.y/FT_PER_M:.15e}  m/s²")

# k3
v3 = velocity_vector + 0.5 * delta_time * a2
rel3 = v3 - wind_vector
a3 = acceleration(rel3, v3)

print()
print("--- k3 ---")
print(f"v3.x:                   {v3.x:.15e}  fps")
print(f"v3.y:                   {v3.y:.15e}  fps")
print(f"a3.x:                   {a3.x:.15e}  fps²")
print(f"a3.y:                   {a3.y:.15e}  fps²")
print(f"a3.x/FT_PER_M:          {a3.x/FT_PER_M:.15e}  m/s²")
print(f"a3.y/FT_PER_M:          {a3.y/FT_PER_M:.15e}  m/s²")

# k4
v4 = velocity_vector + delta_time * a3
rel4 = v4 - wind_vector
a4 = acceleration(rel4, v4)

print()
print("--- k4 ---")
print(f"v4.x:                   {v4.x:.15e}  fps")
print(f"v4.y:                   {v4.y:.15e}  fps")
print(f"a4.x:                   {a4.x:.15e}  fps²")
print(f"a4.y:                   {a4.y:.15e}  fps²")
print(f"a4.x/FT_PER_M:          {a4.x/FT_PER_M:.15e}  m/s²")
print(f"a4.y/FT_PER_M:          {a4.y/FT_PER_M:.15e}  m/s²")

# Final update
new_velocity = velocity_vector + (a1 + 2*a2 + 2*a3 + a4) * (delta_time / 6.0)
new_range = range_vector + (v1 + 2*v2 + 2*v3 + v4) * (delta_time / 6.0)

print()
print("=== AFTER STEP 1 (imperial) ===")
print(f"new_range.x:            {new_range.x:.15e}  ft")
print(f"new_range.y:            {new_range.y:.15e}  ft")
print(f"new_range.z:            {new_range.z:.15e}  ft")
print(f"new_velocity.x:         {new_velocity.x:.15e}  fps")
print(f"new_velocity.y:         {new_velocity.y:.15e}  fps")
print(f"new_velocity.z:         {new_velocity.z:.15e}  fps")

print()
print("=== AFTER STEP 1 (metric, /3.2808399) ===")
print(f"new_range.x:            {new_range.x/FT_PER_M:.15e}  m")
print(f"new_range.y:            {new_range.y/FT_PER_M:.15e}  m")
print(f"new_range.z:            {new_range.z/FT_PER_M:.15e}  m")
print(f"new_velocity.x:         {new_velocity.x/FT_PER_M:.15e}  m/s")
print(f"new_velocity.y:         {new_velocity.y/FT_PER_M:.15e}  m/s")
print(f"new_velocity.z:         {new_velocity.z/FT_PER_M:.15e}  m/s")
print(f"new_speed:              {new_velocity.magnitude()/FT_PER_M:.15e}  m/s")

mcp.sh

#!/bin/bash
PORT="$1"
SESSION_ID="$2"
shift 2
BODY="$*"

curl -s -X POST "http://localhost:${PORT}/mcp" \
  -H "Content-Type: application/json" \
  -H "Accept: application/json, text/event-stream" \
  -H "Mcp-Session-Id: $SESSION_ID" \
  -d "$BODY" 2>&1 | grep '^data:' | sed 's/^data: //'