Add Astra example: ascent+descent validation, hazard zone & abort

README.md

@@ -55,7 +55,11 @@ trochia against a real measured flight: its predicted apogee lands within ~1 % o
 the altimeter reading (39.1 m vs 39 m). The
 [psas-launch12](examples/psas-launch12/) example scales that up to a real ~5 km
 high-power flight (validated apogee) and computes the contingency landing zones —
-nominal recovery vs parachute failure vs CATO.
+nominal recovery vs parachute failure vs CATO. The
+[astra](examples/astra/) example validates a real ~3.25 km flight on **both** the
+ascent and the descent (its recovery failed, so the measured fall is near
+ballistic — exactly what trochia can check) and lays out the hazard zone
+(警戒区域) and abort (途中破談) footprints.
 
 ## Python tooling (uv)
 

examples/README.md

@@ -11,8 +11,8 @@ $ ./fetch-engine.sh        # LARKSPUR-XP.300 (20191020_01.eng), used by
                            # single-trajectory, landing-dispersion and parachute
 ```
 
-`validation-estes-viking` fetches the Estes A8 from ThrustCurve.org with its own
-`fetch-engine.sh`.
+`validation-estes-viking` (Estes A8) and `astra` (Cesaroni 4263L1350-P) fetch
+their own motors with their own `fetch-engine.sh`.
 
 You also need the `trochia` binary — build it from the repo root (see the top
 [README](../README.md)); it ends up at `build/bin/trochia`. trochia always reads
@@ -28,5 +28,6 @@ run each example from its own directory.
 | [parachute](parachute/)                   | parachute vs ballistic descent — how recovery changes the landing point |
 | [validation-estes-viking](validation-estes-viking/) | accuracy check: predicted vs **measured** apogee of a real flight |
 | [psas-launch12](psas-launch12/)           | km-scale real flight: apogee validation + contingency landing zones |
+| [astra](astra/)                           | km-scale real flight: ascent **+ descent** validation, hazard zone (警戒区域) & abort (途中破談) |
 
 See each directory's `README.md` for how to run it and the resulting plot.

examples/astra/README.md

@@ -0,0 +1,114 @@
+# Astra — km-scale real-flight validation (ascent **and** descent) + hazard zone & abort
+
+**Use case:** a real, instrumented flight (~3.25 km apogee) used to (a) validate
+trochia against measured data on **both** the ascent and the descent, and
+(b) compute the hazard area / landing dispersion (警戒区域 / 落下分散) and the
+abort (途中破談) footprints — a full launch-safety analysis on a single rocket.
+
+## The flight
+
+**Astra**, built by **Faraday Rocketry UPV**, flown at **EuRoC 2022** (European
+Rocketry Challenge, Ponte de Sor, Portugal) on a commercial **Cesaroni
+4263L1350-P** solid motor. Measured apogee **3249 m at t ≈ 24 s**; the recovery
+**did not deploy cleanly**, so the rocket came down close to ballistic
+(~100 m/s) — which, conveniently, lets us validate trochia's *descent* leg
+against real data, not just the apogee.
+
+The rocket parameters, the measured trajectory (`measured.dat`,
+`measured-track.dat`) and the drag curve all come from the RocketPy "Astra"
+example:
+- rocket model + measured flight data: <https://github.com/RocketPy-Team/RocketPy/tree/master/data/rockets/astra>
+- worked example notebook: `docs/examples/astra_flight_sim.ipynb` in that repo
+- motor: <https://www.thrustcurve.org/motors/Cesaroni/4263L1350-P/>
+
+## Run
+
+```sh
+./run.sh         # fetch motor, run validation + hazard/abort sims, write the plots
+```
+
+`run.sh` needs the `trochia` binary at `build/bin/trochia` (build from the repo
+root) and, for the map, network access + `uv`.
+
+## Validation vs the measured flight
+
+`config-validation.toml` reconstructs the actual flight in still air and is
+plotted against the rocket's own measured altitude history:
+
+![validation](validation.png)
+
+| quantity | measured | trochia | 
+|---|---|---|
+| apogee | 3249 m | 3260 m (**+0.3 %**) |
+| time to apogee | ~24 s | 25.0 s |
+| descent speed (near ground) | ~100 m/s | ~108 m/s |
+
+- **Ascent (left)** is the validation: with a single representative drag
+  coefficient the whole climb and the apogee match within a fraction of a
+  percent.
+- **Descent (right)** is the part that's usually *not* checkable. Because Astra's
+  recovery only partially deployed, the measured descent is close to ballistic
+  and trochia reproduces it to within ~70 m down to ~145 m AGL, where the
+  measured record ends (see *recovery modelling* below).
+
+### Why a constant Cd works here
+
+trochia uses a single constant drag coefficient, but Astra's real drag varies
+with Mach (subsonic ~0.58, a transonic peak of 0.857 near Mach 0.89, which the
+rocket reaches around burnout). Sweeping the constant Cd in still air shows the
+measured apogee selects an **effective constant Cd ≈ 0.55** — the constant that
+reproduces the measured apogee, i.e. the transonic drag rise folded into one
+number (not the arithmetic mean of the Cd–Mach curve):
+
+![cd-sensitivity](cd-sensitivity.png)
+
+## Hazard zone (落下分散) + abort (途中破談)
+
+`config.toml` flies the validated ascent model under a wind-direction sweep
+(5 m/s, 0–360°) for four contingency branches, giving one landing footprint
+each:
+
+| scenario | meaning | recovery |
+|---|---|---|
+| `nominal` | intended flight | dual deploy (drogue at apogee + main at 450 m AGL) |
+| `recovery-failure` | total recovery failure | none → ballistic (worst-case impact) |
+| `motor-cutoff` | thrust quits mid-burn (t = 1.5 s) | airframe intact, recovers |
+| `cato` | structural failure (t = 1.0 s) | thrust stops + no recovery |
+
+![footprints](footprints.png)
+
+On the real launch site (the measured GPS track, up to where its record ends at
+~145 m AGL, is drawn for reference):
+
+![footprints map](footprints-map.png)
+
+The launch is tilted (84° elevation, 133° heading), so every footprint is biased
+**downrange toward 133°**. Reading the safety story off the plot:
+
+- **dual deploy** drifts the *most* — it spends minutes under canopy, so wind
+  exposure dominates and the zone is the widest, though the landing is gentle;
+- **recovery failure** (no recovery at all) is a tight, downrange *lawn-dart* —
+  it falls too fast to drift much, but hits hard. The real flight's partial
+  recovery (~100 m/s) was a milder version of this; a clean total failure would
+  be faster and tighter still;
+- **motor cutoff** stays closer in (lower apogee, still recovers);
+- **CATO** debris stays right by the pad.
+
+## Modelling notes / caveats
+
+- **Representative Cd.** 0.55 is the effective constant that reproduces the
+  measured apogee; trochia cannot reproduce the *timing* of the transonic drag
+  rise, only its net effect on apogee (which is what `cd-sensitivity.png`
+  calibrates).
+- **CATO is single-body.** trochia models the `cato` branch as one intact body
+  with thrust cut and no recovery, not a fragmenting debris cloud; read its
+  footprint as "where the (single) vehicle comes down", not a debris dispersion.
+- **Recovery modelling.** Astra's recovery hardware drag is not published. The
+  validation run models the as-flown partial recovery as an *effective* drag
+  `cd·area ≈ 0.015 m²` calibrated to the measured descent; `recovery-failure`
+  in the hazard analysis is a stricter, fully ballistic fall (faster still).
+- **Mass / motor.** trochia's `mass` is the airframe **without** the loaded
+  motor (9.10 kg); the motor (loaded 3.57 kg, 2.02 kg propellant) comes from the
+  `.eng`, giving a 12.67 kg liftoff mass that matches the RocketPy model.
+- **CG / CP / inertia / Cna** are estimated from the RocketPy geometry and are
+  second-order for these (apogee- and descent-dominated) results.

examples/astra/cd-apogee.dat

@@ -0,0 +1,8 @@
+# Cd apogee[m] t_apogee[s]
+0.45 3529.8 26.2
+0.50 3388.3 25.5
+0.55 3260.1 25.0
+0.58 3188.8 24.5
+0.65 3036.5 23.8
+0.75 2847.6 23.0
+0.85 2685.5 22.0

examples/astra/cd-sensitivity.gp

@@ -0,0 +1,21 @@
+# Apogee sensitivity to the constant drag coefficient.
+# cd-apogee.dat ("Cd apogee[m] t_apogee[s]") is produced by run.sh, which sweeps
+# Cd in still air. The measured apogee picks out the representative Cd.
+set datafile separator whitespace
+set terminal pngcairo size 760,560 enhanced font "sans,11"
+set output "cd-sensitivity.png"
+
+set title "Astra apogee vs constant Cd\nthe measured apogee picks the effective constant Cd = 0.55"
+set xlabel "constant drag coefficient Cd [-]"
+set ylabel "apogee [m]"
+set grid
+set key top right
+
+MEAS = 3249.0
+set arrow from graph 0, first MEAS to graph 1, first MEAS nohead lw 2 dt 2 lc rgb "#e08214"
+set label "measured 3249 m" at graph 0.03, first MEAS+90 tc rgb "#b35806"
+set arrow from 0.55, graph 0 to 0.55, graph 1 nohead dt 3 lc rgb "#888888"
+set label "Cd=0.55" at 0.56, graph 0.10 tc rgb "#555555"
+
+plot "cd-apogee.dat" using 1:2 with linespoints pt 7 ps 1.2 lw 2 lc rgb "#31688e" \
+       title "trochia (still air)"

examples/astra/config-validation.toml

@@ -0,0 +1,48 @@
+# Astra validation: reconstruct the actual EuRoC'22 flight and compare to the
+# rocket's own measured trajectory (measured.dat, from RocketPy's flight_data).
+#
+# Still air, single flight. The descent uses an *effective* recovery drag
+# calibrated to the measured fall: Astra's recovery did not deploy cleanly, so
+# the rocket came down at ~100 m/s (tumbling / partial deployment) -- far slower
+# than a clean nose-first ballistic fall (~185 m/s) but far faster than a working
+# main. cd*area ~ 0.015 m^2 reproduces the measured altitude history.
+#
+#   ./fetch-engine.sh
+#   ../../build/bin/trochia config-validation.toml
+#   gnuplot validation.gp        # -> validation.png
+
+[simulation]
+dt			= 0.005
+output.dt	= 0.25
+output.dir	= "output-validation"
+timeout		= 120
+
+[launcher]
+length		= 12             # 12 m rail (RocketPy Flight rail_length)
+elevation	= 84             # inclination 84 deg from horizontal
+azimuth		= 133            # heading 133 deg
+
+[wind]
+model		= "power"
+ground.dir	= [ 0.0 ]
+ground.speed= 0.001          # still air (validate the vertical profile)
+
+[rocket]
+name	= "Astra (Faraday Rocketry UPV)"
+type	= "single"
+
+[[rocket.stage]]
+engine	= "Cesaroni_4263L1350-P.eng"   # Cesaroni 4263L1350-P, fetched above
+length	= 2.52               # nose tip (2.5214) to tail (~0), RocketPy geometry
+diameter= 0.094              # radius 0.047 m
+mass	= 9.10               # dry 10.6462 kg (incl. casing) minus casing 1.546 kg
+I0		= 7.15               # loaded pitch inertia (dry 6.70 + propellant, est.)
+If		= 6.70               # empty pitch inertia (RocketPy dry inertia 6.6986)
+lcg0	= 1.21               # loaded CG from nose (est. from RocketPy masses)
+lcgf	= 1.126              # empty CG from nose (CoM-without-motor 1.3957 from tail)
+lcgp	= 1.66               # propellant CG from nose (motor at aft)
+lcp		= 1.60               # CP from nose (near fins, est.)
+Cd		= 0.55               # representative constant Cd (see README; sweep in run.sh)
+Cna		= 14.0               # estimated (4 fins + nose)
+# as-flown effective recovery: ~100 m/s descent (cd*area ~ 0.015 m^2)
+parachute= { condition = "top", delay = 0.0, cd = 1.0, diameter = 0.14 }

examples/astra/config.toml

@@ -0,0 +1,62 @@
+# Astra hazard zone (警戒区域 / 落下分散) + abort (途中破談) analysis.
+#
+# The validated ascent model (see config-validation.toml) is flown under a wind
+# direction sweep for several contingency branches, giving one landing footprint
+# per scenario -- the core of a launch-safety / hazard-area analysis.
+#
+#   ./fetch-engine.sh
+#   ../../build/bin/trochia config.toml          # all scenarios x wind dirs
+#   gnuplot footprints.gp                        # -> footprints.png
+#   uv run footprints-map.py                     # -> footprints-map.png
+
+[simulation]
+dt			= 0.005
+output.dt	= 0.5
+output.dir	= "output"
+timeout		= 400          # a dual-deploy descent from ~3.2 km takes minutes
+
+[launcher]
+length		= 12
+elevation	= 84
+azimuth		= 133          # downrange toward heading 133 deg
+
+[wind]
+model		= "power"
+ground.dir	= [ 0.0, 30.0, 60.0, 90.0, 120.0, 150.0, 180.0, 210.0, 240.0, 270.0, 300.0, 330.0, 360.0 ]
+ground.speed= 5.0
+
+[rocket]
+name	= "Astra (Faraday Rocketry UPV)"
+type	= "single"
+
+[[rocket.stage]]
+engine	= "Cesaroni_4263L1350-P.eng"
+length	= 2.52
+diameter= 0.094
+mass	= 9.10
+I0		= 7.15
+If		= 6.70
+lcg0	= 1.21
+lcgf	= 1.126
+lcgp	= 1.66
+lcp		= 1.60
+Cd		= 0.55
+Cna		= 14.0
+# intended recovery: dual deploy -- drogue at apogee, main at 450 m AGL.
+parachute= [
+	{ condition = "top",  delay = 0.0, cd = 1.5, diameter = 0.5 },   # drogue (~24 m/s)
+	{ altitude = 450.0,                cd = 1.5, diameter = 1.5 },   # main  (~8 m/s)
+]
+
+# contingency branches -> one landing footprint each
+[[scenario]]
+name = "nominal"                  # dual deploy works
+[[scenario]]
+name = "recovery-failure"         # no recovery -> ballistic (worst case impact)
+recovery = "none"
+[[scenario]]
+name = "motor-cutoff"             # thrust quits mid-burn, airframe intact (recovers)
+motor_cutoff = 1.5
+[[scenario]]
+name = "cato"                     # structural failure: thrust stops + no recovery
+cato = 1.0

examples/astra/fetch-engine.sh

@@ -0,0 +1,10 @@
+#!/usr/bin/env bash
+# Download the Cesaroni 4263L1350-P motor thrust curve (RASP .eng).
+# This is the exact file used by RocketPy's "Astra" model, so the motor matches
+# the rocket parameters in config.toml.
+set -euo pipefail
+cd "$(dirname "$0")"
+URL="https://raw.githubusercontent.com/RocketPy-Team/RocketPy/master/data/motors/cesaroni/Cesaroni_4263L1350-P.eng"
+echo "downloading Cesaroni_4263L1350-P.eng ..."
+curl -fsSL "$URL" -o Cesaroni_4263L1350-P.eng
+echo "saved to $(pwd)/Cesaroni_4263L1350-P.eng"

examples/astra/footprints-map.py

@@ -0,0 +1,97 @@
+#!/usr/bin/env python
+"""Overlay the contingency landing footprints + the real GPS track on the map.
+
+Reads each scenario's ground-hit points (output/<scenario>/84/ghp.csv), converts
+the local ENU coordinates to lat/lon (pymap3d), and draws one footprint per
+scenario on OpenStreetMap tiles (staticmap). The rocket's measured GPS track
+(measured-track.dat, up to ~145 m AGL where the record ends) is drawn too.
+
+Run via uv after trochia:   uv run footprints-map.py   # -> footprints-map.png
+"""
+import csv
+
+import pymap3d as pm
+from PIL import ImageDraw, ImageFont
+from staticmap import StaticMap, Line, CircleMarker
+
+# Launch site: Ponte de Sor, Portugal (EuRoC). First GPS fix in flight_data.
+LAUNCH_LAT, LAUNCH_LON, LAUNCH_ALT = 39.3898, -8.290, 160.0
+
+SCENARIOS = [
+    ("nominal",          "nominal (dual deploy)",      "#31688e"),
+    ("recovery-failure", "recovery failure (ballistic)", "#cc0000"),
+    ("motor-cutoff",     "motor cutoff @ 1.5 s",       "#dd8800"),
+    ("cato",             "CATO @ 1.0 s",               "#4ac16d"),
+]
+
+
+def enu_to_lonlat(e, n):
+    x, y, z = pm.enu2ecef(e, n, 0.0, LAUNCH_LAT, LAUNCH_LON, LAUNCH_ALT)
+    lat, lon, _ = pm.ecef2geodetic(x, y, z)
+    return lon, lat
+
+
+def loop(scenario):
+    pts = []
+    with open(f"output/{scenario}/84/ghp.csv") as f:
+        for row in csv.DictReader(f, skipinitialspace=True):
+            pts.append(enu_to_lonlat(float(row["ghp_e"]), float(row["ghp_n"])))
+    return pts
+
+
+def measured_track():
+    pts = []
+    with open("measured-track.dat") as f:
+        for line in f:
+            if line.startswith("#"):
+                continue
+            _, lat, lon, _ = line.split()
+            lat, lon = float(lat), float(lon)
+            if abs(lat) < 1.0 and abs(lon) < 1.0:
+                continue  # GPS dropout fix (lat=lon~0); skip
+            pts.append((lon, lat))
+    return pts
+
+
+def find_font(size):
+    for p in ("/usr/share/fonts/TTF/Roboto-Regular.ttf",
+              "/usr/share/fonts/truetype/dejavu/DejaVuSans.ttf"):
+        try:
+            return ImageFont.truetype(p, size)
+        except OSError:
+            continue
+    return ImageFont.load_default()
+
+
+def main():
+    m = StaticMap(950, 950, url_template="https://tile.openstreetmap.org/{z}/{x}/{y}.png")
+    for name, _, color in SCENARIOS:
+        m.add_line(Line(loop(name), color, 3))
+    # measured GPS track (white casing + dark line), ends at ~145 m AGL
+    track = measured_track()
+    m.add_line(Line(track, "white", 5))
+    m.add_line(Line(track, "#222222", 2))
+    m.add_marker(CircleMarker(enu_to_lonlat(0.0, 0.0), "black", 13))
+    m.add_marker(CircleMarker(enu_to_lonlat(0.0, 0.0), "white", 7))
+
+    img = m.render().convert("RGB")
+    d = ImageDraw.Draw(img, "RGBA")
+    font, title = find_font(15), find_font(16)
+    rows = SCENARIOS + [(None, "measured GPS track", "#222222")]
+    w, h = 250, 24 * (len(rows) + 1) + 34
+    d.rectangle([10, 10, 10 + w, 10 + h], fill=(255, 255, 255, 225), outline=(0, 0, 0, 255))
+    d.text((20, 18), "contingency footprints", font=title, fill=(0, 0, 0))
+    for i, (_, label, color) in enumerate(rows):
+        y = 44 + i * 24
+        d.line([(20, y + 8), (50, y + 8)], fill=color, width=4)
+        d.text((58, y), label, font=font, fill=(0, 0, 0))
+    y = 44 + len(rows) * 24
+    d.ellipse([30, y + 2, 40, y + 12], fill="black")
+    d.text((58, y), "launch (Ponte de Sor)", font=font, fill=(0, 0, 0))
+
+    img.save("footprints-map.png")
+    print("saved footprints-map.png")
+
+
+if __name__ == "__main__":
+    main()