From 2bbb1871cc7ef8cfbdc6a6b6f2bd201b26198bd7 Mon Sep 17 00:00:00 2001
From: Dane Sabo <yourstruly@danesabo.com>
Date: Tue, 21 Apr 2026 20:24:48 -0400
Subject: [PATCH] refactor: scripts subdivision + TOML configs + results/ split
 + presentation outline
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Architecture restructure from morning review:

1. code/scripts/ subdivided into sim/, reach/, barrier/, plot/.
   Easier nav; `barrier/` is the natural place for SOS scale-up scripts.
2. Heatup PJ reach variants consolidated behind TOML configs.
   reach_heatup_pj.jl now takes `--config path/to/config.toml`;
   configs/heatup/baseline.toml (wide entry, from predicates.json) and
   configs/heatup/tight.toml (narrow entry, reproduces all-6-halfspaces
   discharged result). Old reach_heatup_pj_tight.jl and
   reach_heatup_pj_tight_full.jl deleted (superseded).
3. Reach output .mat files moved from reachability/ to results/.
   reachability/ now = specs + docs; results/ = ephemeral outputs
   (gitignored *.mat). README added.
4. OVERNIGHT_NOTES.md archived to claude_memory/2026-04-20-21-overnight-
   session-summary.md (date range in the filename makes the history clearer).

All include() / Pkg.activate() paths in scripts updated for the new
depth. Smoke tests pass (reach_operation.jl generates its .mat in
the new results/ location; sim_sanity.jl matches MATLAB).

Presentation outline for the 20-min prelim talk landed in
presentations/prelim-presentation/outline.md. 14-slide assertion-
evidence format targeting OT-informed cybersecurity audience. Each
slide: one declarative assertion + one figure. Outline includes
which figures already exist and which need to be created, timing
checkpoints, cybersecurity angle to emphasize, and Q&A prep.

New config configs/heatup/with_steam_dump.toml + its companion
scripts/reach/reach_heatup_pj_sd.jl (12-state RHS with Q_sg as an
augmented bounded parameter x[10] and time as x[11]). Kicks off
point 3 from morning review.

Next up: scram X_entry expansion (morning point 2) — LOCA scenario
+ union of mode reach envelopes.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 app/predicate_explorer.jl                     |   4 +-
 ...2026-04-20-21-overnight-session-summary.md |   0
 code/configs/heatup/baseline.toml             |  22 +
 code/configs/heatup/tight.toml                |  27 ++
 code/configs/heatup/with_steam_dump.toml      |  40 ++
 .../{ => barrier}/barrier_compare_OL_CL.jl    |  10 +-
 .../scripts/{ => barrier}/barrier_lyapunov.jl |  10 +-
 .../{ => barrier}/barrier_polytopic.jl        |  12 +-
 code/scripts/{ => barrier}/barrier_sos_2d.jl  |   8 +-
 code/scripts/{ => plot}/plot_reach_tubes.jl   |  12 +-
 .../{ => reach}/reach_heatup_nonlinear.jl     |   4 +-
 code/scripts/reach/reach_heatup_pj.jl         | 269 ++++++++++++
 code/scripts/reach/reach_heatup_pj_sd.jl      | 178 ++++++++
 code/scripts/{ => reach}/reach_operation.jl   |  16 +-
 code/scripts/{ => reach}/reach_scram_pj.jl    |   4 +-
 code/scripts/reach_heatup_pj.jl               | 243 -----------
 code/scripts/reach_heatup_pj_tight.jl         | 109 -----
 code/scripts/reach_heatup_pj_tight_full.jl    | 145 -------
 code/scripts/{ => sim}/main_mode_sweep.jl     |  18 +-
 code/scripts/{ => sim}/sim_heatup.jl          |   2 +-
 code/scripts/{ => sim}/sim_sanity.jl          |   8 +-
 code/scripts/{ => sim}/validate_pj.jl         |  12 +-
 presentations/prelim-presentation/outline.md  | 395 ++++++++++++++++++
 results/.gitignore                            |   1 +
 results/README.md                             |  32 ++
 25 files changed, 1024 insertions(+), 557 deletions(-)
 rename OVERNIGHT_NOTES.md => claude_memory/2026-04-20-21-overnight-session-summary.md (100%)
 create mode 100644 code/configs/heatup/baseline.toml
 create mode 100644 code/configs/heatup/tight.toml
 create mode 100644 code/configs/heatup/with_steam_dump.toml
 rename code/scripts/{ => barrier}/barrier_compare_OL_CL.jl (90%)
 rename code/scripts/{ => barrier}/barrier_lyapunov.jl (92%)
 rename code/scripts/{ => barrier}/barrier_polytopic.jl (94%)
 rename code/scripts/{ => barrier}/barrier_sos_2d.jl (96%)
 rename code/scripts/{ => plot}/plot_reach_tubes.jl (94%)
 rename code/scripts/{ => reach}/reach_heatup_nonlinear.jl (98%)
 create mode 100644 code/scripts/reach/reach_heatup_pj.jl
 create mode 100644 code/scripts/reach/reach_heatup_pj_sd.jl
 rename code/scripts/{ => reach}/reach_operation.jl (91%)
 rename code/scripts/{ => reach}/reach_scram_pj.jl (98%)
 delete mode 100644 code/scripts/reach_heatup_pj.jl
 delete mode 100644 code/scripts/reach_heatup_pj_tight.jl
 delete mode 100644 code/scripts/reach_heatup_pj_tight_full.jl
 rename code/scripts/{ => sim}/main_mode_sweep.jl (89%)
 rename code/scripts/{ => sim}/sim_heatup.jl (98%)
 rename code/scripts/{ => sim}/sim_sanity.jl (83%)
 rename code/scripts/{ => sim}/validate_pj.jl (91%)
 create mode 100644 presentations/prelim-presentation/outline.md
 create mode 100644 results/.gitignore
 create mode 100644 results/README.md

diff --git a/app/predicate_explorer.jl b/app/predicate_explorer.jl
index d1dacc8..913ca50 100644
--- a/app/predicate_explorer.jl
+++ b/app/predicate_explorer.jl
@@ -324,7 +324,7 @@ no more hand-maintained traceability table.
 
 # ╔═╡ 37d6b212-9f00-4684-9f91-50c7e17cbd62
 begin
-    reach_mat_path = joinpath(@__DIR__, "..", "reachability", "reach_operation_result.mat")
+    reach_mat_path = joinpath(@__DIR__, "..", "results", "reach_operation_result.mat")
     have_reach_mat = isfile(reach_mat_path)
     if have_reach_mat
         reach_mat = matread(reach_mat_path)
@@ -463,7 +463,7 @@ extend past the ~10 s prompt-neutron stiffness wall.
 
 # ╔═╡ 45e8962a-873e-48d3-aac4-70ea0cf36c97
 begin
-    pj_mat_path = joinpath(@__DIR__, "..", "reachability", "reach_heatup_pj_result.mat")
+    pj_mat_path = joinpath(@__DIR__, "..", "results", "reach_heatup_pj_result.mat")
     if isfile(pj_mat_path)
         pj_mat = matread(pj_mat_path)
         md"""**Loaded** `reach_heatup_pj_result.mat` — pending plot rendering."""
diff --git a/OVERNIGHT_NOTES.md b/claude_memory/2026-04-20-21-overnight-session-summary.md
similarity index 100%
rename from OVERNIGHT_NOTES.md
rename to claude_memory/2026-04-20-21-overnight-session-summary.md
diff --git a/code/configs/heatup/baseline.toml b/code/configs/heatup/baseline.toml
new file mode 100644
index 0000000..dd0539c
--- /dev/null
+++ b/code/configs/heatup/baseline.toml
@@ -0,0 +1,22 @@
+# Baseline heatup-mode reach configuration.
+# Wide X_entry as specified in predicates.json mode_boundaries.
+# Matches the original reach_heatup_pj.jl defaults.
+
+name = "baseline"
+description = "Wide X_entry from mode_boundaries.q_heatup.X_entry_polytope"
+
+# Load X_entry from predicates.json (true) or from this file (false).
+use_predicates_entry = true
+
+[tmjets]
+orderT = 4
+orderQ = 2
+abstol = 1e-9
+maxsteps = 100000
+
+[probes]
+horizons_seconds = [60.0, 300.0, 1800.0, 5400.0]
+
+[output]
+save_per_step = true
+result_file = "reach_heatup_pj_baseline.mat"
diff --git a/code/configs/heatup/tight.toml b/code/configs/heatup/tight.toml
new file mode 100644
index 0000000..df3dce9
--- /dev/null
+++ b/code/configs/heatup/tight.toml
@@ -0,0 +1,27 @@
+# Tight X_entry heatup reach: T_c in [285, 291] (6 K vs 14 K baseline).
+# Produces the clean "all 6 inv1_holds halfspaces discharged at T=300s"
+# result from 2026-04-20 overnight session.
+
+name = "tight"
+description = "Narrow X_entry on T_c/T_f/T_cold; n in [1e-3, 2e-3]"
+
+use_predicates_entry = false
+
+[entry]
+n_range = [1.0e-3, 2.0e-3]
+T_f_range_C = [285.0, 291.0]
+T_c_range_C = [285.0, 291.0]
+T_cold_range_C = [278.0, 285.0]
+
+[tmjets]
+orderT = 4
+orderQ = 2
+abstol = 1e-9
+maxsteps = 100000
+
+[probes]
+horizons_seconds = [60.0, 300.0]
+
+[output]
+save_per_step = true
+result_file = "reach_heatup_pj_tight.mat"
diff --git a/code/configs/heatup/with_steam_dump.toml b/code/configs/heatup/with_steam_dump.toml
new file mode 100644
index 0000000..a00400e
--- /dev/null
+++ b/code/configs/heatup/with_steam_dump.toml
@@ -0,0 +1,40 @@
+# Heatup reach with a bounded secondary-side steam-dump Q_sg.
+#
+# Instead of Q_sg ≡ 0 (original assumption), treat Q_sg as a bounded
+# disturbance in [0, 0.05·P_0].  Physical interpretation: operator
+# opens/closes the secondary-side steam dump to manage primary
+# temperature during the ramp; exact value not known, but bounded
+# by atmospheric-dump capacity (~5% of P_0 rated).
+#
+# The reach script picks up Q_sg as an augmented state x[11] with
+# dx[11] = 0, entry box covering [0, 0.05*P_0].
+
+name = "with_steam_dump"
+description = "Tight X_entry + bounded Q_sg ∈ [0, 0.05·P_0] as disturbance"
+
+use_predicates_entry = false
+steam_dump_enabled = true
+
+[entry]
+n_range = [1.0e-3, 2.0e-3]
+T_f_range_C = [285.0, 291.0]
+T_c_range_C = [285.0, 291.0]
+T_cold_range_C = [278.0, 285.0]
+
+[entry.steam_dump]
+# Q_sg bounds in fractions of P_0.
+Q_lo_fraction_P0 = 0.0
+Q_hi_fraction_P0 = 0.05
+
+[tmjets]
+orderT = 4
+orderQ = 2
+abstol = 1e-9
+maxsteps = 100000
+
+[probes]
+horizons_seconds = [60.0, 300.0]
+
+[output]
+save_per_step = true
+result_file = "reach_heatup_pj_with_steam_dump.mat"
diff --git a/code/scripts/barrier_compare_OL_CL.jl b/code/scripts/barrier/barrier_compare_OL_CL.jl
similarity index 90%
rename from code/scripts/barrier_compare_OL_CL.jl
rename to code/scripts/barrier/barrier_compare_OL_CL.jl
index 5a6102d..5c33280 100644
--- a/code/scripts/barrier_compare_OL_CL.jl
+++ b/code/scripts/barrier/barrier_compare_OL_CL.jl
@@ -8,17 +8,17 @@
 # anisotropy, not controller tuning.
 
 using Pkg
-Pkg.activate(joinpath(@__DIR__, ".."))
+Pkg.activate(joinpath(@__DIR__, "..", ".."))
 
 using Printf
 using LinearAlgebra
 using MatrixEquations
 using JSON
 
-include(joinpath(@__DIR__, "..", "src", "pke_params.jl"))
-include(joinpath(@__DIR__, "..", "src", "pke_th_rhs.jl"))
-include(joinpath(@__DIR__, "..", "src", "pke_linearize.jl"))
-include(joinpath(@__DIR__, "..", "src", "load_predicates.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_params.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_th_rhs.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_linearize.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "load_predicates.jl"))
 
 plant = pke_params()
 x_op  = pke_initial_conditions(plant)
diff --git a/code/scripts/barrier_lyapunov.jl b/code/scripts/barrier/barrier_lyapunov.jl
similarity index 92%
rename from code/scripts/barrier_lyapunov.jl
rename to code/scripts/barrier/barrier_lyapunov.jl
index 7de835c..dcd3926 100644
--- a/code/scripts/barrier_lyapunov.jl
+++ b/code/scripts/barrier/barrier_lyapunov.jl
@@ -10,17 +10,17 @@
 # This is the structural anisotropy limitation, not a code bug.
 
 using Pkg
-Pkg.activate(joinpath(@__DIR__, ".."))
+Pkg.activate(joinpath(@__DIR__, "..", ".."))
 
 using Printf
 using LinearAlgebra
 using MatrixEquations
 using JSON
 
-include(joinpath(@__DIR__, "..", "src", "pke_params.jl"))
-include(joinpath(@__DIR__, "..", "src", "pke_th_rhs.jl"))
-include(joinpath(@__DIR__, "..", "src", "pke_linearize.jl"))
-include(joinpath(@__DIR__, "..", "src", "load_predicates.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_params.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_th_rhs.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_linearize.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "load_predicates.jl"))
 
 plant = pke_params()
 x_op  = pke_initial_conditions(plant)
diff --git a/code/scripts/barrier_polytopic.jl b/code/scripts/barrier/barrier_polytopic.jl
similarity index 94%
rename from code/scripts/barrier_polytopic.jl
rename to code/scripts/barrier/barrier_polytopic.jl
index f1e3706..605f2d7 100644
--- a/code/scripts/barrier_polytopic.jl
+++ b/code/scripts/barrier/barrier_polytopic.jl
@@ -21,7 +21,7 @@
 # Uses JuMP + HiGHS (open-source, ships with no license trouble).
 
 using Pkg
-Pkg.activate(joinpath(@__DIR__, ".."))
+Pkg.activate(joinpath(@__DIR__, "..", ".."))
 
 using Printf
 using LinearAlgebra
@@ -31,10 +31,10 @@ using JSON
 using JuMP
 using HiGHS
 
-include(joinpath(@__DIR__, "..", "src", "pke_params.jl"))
-include(joinpath(@__DIR__, "..", "src", "pke_th_rhs.jl"))
-include(joinpath(@__DIR__, "..", "src", "pke_linearize.jl"))
-include(joinpath(@__DIR__, "..", "src", "load_predicates.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_params.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_th_rhs.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_linearize.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "load_predicates.jl"))
 
 plant = pke_params()
 x_op  = pke_initial_conditions(plant)
@@ -55,7 +55,7 @@ b_inv2 = inv2.b_poly        # 6
 
 # --- Precursor bounds from reach-tube envelope ---
 # Read reach_operation_result.mat; take min/max of X_lo, X_hi on precursors.
-reach_mat_path = joinpath(@__DIR__, "..", "..", "reachability",
+reach_mat_path = joinpath(@__DIR__, "..", "..", "..", "reachability",
                           "reach_operation_result.mat")
 reach = matread(reach_mat_path)
 X_lo = reach["X_lo"]; X_hi = reach["X_hi"]
diff --git a/code/scripts/barrier_sos_2d.jl b/code/scripts/barrier/barrier_sos_2d.jl
similarity index 96%
rename from code/scripts/barrier_sos_2d.jl
rename to code/scripts/barrier/barrier_sos_2d.jl
index 2cabae5..7ad302c 100644
--- a/code/scripts/barrier_sos_2d.jl
+++ b/code/scripts/barrier/barrier_sos_2d.jl
@@ -24,7 +24,7 @@
 # Using SOS multipliers σ_i(x), w-dependence via lossless-disturbance bound.
 
 using Pkg
-Pkg.activate(joinpath(@__DIR__, ".."))
+Pkg.activate(joinpath(@__DIR__, "..", ".."))
 
 using LinearAlgebra
 using MatrixEquations
@@ -32,9 +32,9 @@ using DynamicPolynomials
 using SumOfSquares
 using CSDP
 
-include(joinpath(@__DIR__, "..", "src", "pke_params.jl"))
-include(joinpath(@__DIR__, "..", "src", "pke_th_rhs.jl"))
-include(joinpath(@__DIR__, "..", "src", "pke_linearize.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_params.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_th_rhs.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_linearize.jl"))
 
 plant = pke_params()
 x_op  = pke_initial_conditions(plant)
diff --git a/code/scripts/plot_reach_tubes.jl b/code/scripts/plot/plot_reach_tubes.jl
similarity index 94%
rename from code/scripts/plot_reach_tubes.jl
rename to code/scripts/plot/plot_reach_tubes.jl
index 76de00c..998107b 100644
--- a/code/scripts/plot_reach_tubes.jl
+++ b/code/scripts/plot/plot_reach_tubes.jl
@@ -22,17 +22,17 @@
 #   julia --project=. scripts/plot_reach_tubes.jl heatup_pj
 
 using Pkg
-Pkg.activate(joinpath(@__DIR__, ".."))
+Pkg.activate(joinpath(@__DIR__, "..", ".."))
 
 using MAT
 using Plots
 gr()
 
-include(joinpath(@__DIR__, "..", "src", "pke_params.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_params.jl"))
 const PLANT = pke_params()
 
 function plot_tubes_operation()
-    mat_path = joinpath(@__DIR__, "..", "..", "reachability", "reach_operation_result.mat")
+    mat_path = joinpath(@__DIR__, "..", "..", "..", "results", "reach_operation_result.mat")
     d = matread(mat_path)
 
     T_vec = vec(d["T"])                        # time grid
@@ -75,7 +75,7 @@ function plot_tubes_operation()
 end
 
 function plot_tubes_heatup_pj()
-    mat_path = joinpath(@__DIR__, "..", "..", "reachability",
+    mat_path = joinpath(@__DIR__, "..", "..", "..", "reachability",
                        "reach_heatup_pj_tight_full.mat")
     d = matread(mat_path)
 
@@ -144,7 +144,7 @@ function _plot_common(t, Tc_lo, Tc_hi, Th_lo, Th_hi, Tco_lo, Tco_hi,
 
     fig = plot(p1, p2, p3, p4, layout=(2, 2), size=(1300, 800),
                plot_title=title_stem)
-    figdir = joinpath(@__DIR__, "..", "..", "docs", "figures")
+    figdir = joinpath(@__DIR__, "..", "..", "..", "docs", "figures")
     isdir(figdir) || mkpath(figdir)
     outpath = joinpath(figdir, outname)
     savefig(fig, outpath)
@@ -157,7 +157,7 @@ if which_plot in ("operation", "both")
     plot_tubes_operation()
 end
 if which_plot in ("heatup_pj", "both")
-    mat_path = joinpath(@__DIR__, "..", "..", "reachability",
+    mat_path = joinpath(@__DIR__, "..", "..", "..", "reachability",
                        "reach_heatup_pj_tight_full.mat")
     if isfile(mat_path)
         plot_tubes_heatup_pj()
diff --git a/code/scripts/reach_heatup_nonlinear.jl b/code/scripts/reach/reach_heatup_nonlinear.jl
similarity index 98%
rename from code/scripts/reach_heatup_nonlinear.jl
rename to code/scripts/reach/reach_heatup_nonlinear.jl
index 330f3fc..0ef2c0d 100644
--- a/code/scripts/reach_heatup_nonlinear.jl
+++ b/code/scripts/reach/reach_heatup_nonlinear.jl
@@ -26,7 +26,7 @@
 # Time carried as augmented state x[11].
 
 using Pkg
-Pkg.activate(joinpath(@__DIR__, ".."))
+Pkg.activate(joinpath(@__DIR__, "..", ".."))
 
 using LinearAlgebra
 using ReachabilityAnalysis, LazySets
@@ -101,7 +101,7 @@ const KP_HEATUP    = 1e-4
 end
 
 # --- Build X_entry from predicates.json ---
-pred_path = joinpath(@__DIR__, "..", "..", "reachability", "predicates.json")
+pred_path = joinpath(@__DIR__, "..", "..", "..", "reachability", "predicates.json")
 pred_raw  = JSON.parsefile(pred_path)
 entry     = pred_raw["mode_boundaries"]["q_heatup"]["X_entry_polytope"]
 
diff --git a/code/scripts/reach/reach_heatup_pj.jl b/code/scripts/reach/reach_heatup_pj.jl
new file mode 100644
index 0000000..f5645f6
--- /dev/null
+++ b/code/scripts/reach/reach_heatup_pj.jl
@@ -0,0 +1,269 @@
+#!/usr/bin/env julia
+#
+# reach_heatup_pj.jl — nonlinear reach on heatup, prompt-jump model.
+#
+# Reduced from 10-state to 9-state (n is algebraic). Removes the Λ⁻¹
+# stiffness that capped the full-state reach at ~10 s.
+#
+# State (10D with augmented time):
+#   x[1..6] = C_1..C_6   (delayed-neutron precursors)
+#   x[7]    = T_f
+#   x[8]    = T_c
+#   x[9]    = T_cold
+#   x[10]   = t (augmented time, dt/dt = 1)
+#
+# n is algebraic:  n = Λ·Σ λ_i C_i / (β - ρ),  ρ = K_p·(T_ref - T_c).
+#
+# Configuration-driven: pass a TOML config path as the first CLI arg,
+# or omit for the baseline config.
+#
+#   julia --project=. scripts/reach/reach_heatup_pj.jl                  # baseline
+#   julia --project=. scripts/reach/reach_heatup_pj.jl configs/heatup/tight.toml
+#
+# Configs live in code/configs/heatup/*.toml. See baseline.toml for
+# the full schema.
+
+using Pkg
+Pkg.activate(joinpath(@__DIR__, "..", ".."))
+
+using LinearAlgebra
+using ReachabilityAnalysis, LazySets
+using JSON
+using MAT
+using TOML
+
+# --- Plant constants (must match pke_params) ---
+const LAMBDA  = 1e-4
+const BETA_1, BETA_2, BETA_3, BETA_4, BETA_5, BETA_6 =
+    0.000215, 0.001424, 0.001274, 0.002568, 0.000748, 0.000273
+const BETA    = BETA_1 + BETA_2 + BETA_3 + BETA_4 + BETA_5 + BETA_6
+const LAM_1, LAM_2, LAM_3, LAM_4, LAM_5, LAM_6 =
+    0.0124, 0.0305, 0.111, 0.301, 1.14, 3.01
+
+const P0   = 1e9
+const M_F, C_F, M_C, C_C, HA, W_M, M_SG =
+    50000.0, 300.0, 20000.0, 5450.0, 5e7, 5000.0, 30000.0
+
+const T_COLD0 = 290.0
+const DT_CORE = P0 / (W_M * C_C)
+const T_HOT0  = T_COLD0 + DT_CORE
+const T_C0    = (T_HOT0 + T_COLD0) / 2
+const T_F0    = T_C0 + P0 / HA
+
+const T_STANDBY    = T_C0 - 33.333333
+const RAMP_RATE_CS = 28.0 / 3600
+const KP_HEATUP    = 1e-4
+
+# --- Taylorized heatup PJ RHS ---
+@taylorize function rhs_heatup_pj_taylor!(dx, x, p, t)
+    rho = KP_HEATUP * (T_STANDBY + RAMP_RATE_CS * x[10] - x[8])
+    sum_lam_C = LAM_1*x[1] + LAM_2*x[2] + LAM_3*x[3] +
+                LAM_4*x[4] + LAM_5*x[5] + LAM_6*x[6]
+    denom = BETA - rho
+    n = LAMBDA * sum_lam_C / denom
+    inv_factor = sum_lam_C / denom
+
+    dx[1] = BETA_1 * inv_factor - LAM_1 * x[1]
+    dx[2] = BETA_2 * inv_factor - LAM_2 * x[2]
+    dx[3] = BETA_3 * inv_factor - LAM_3 * x[3]
+    dx[4] = BETA_4 * inv_factor - LAM_4 * x[4]
+    dx[5] = BETA_5 * inv_factor - LAM_5 * x[5]
+    dx[6] = BETA_6 * inv_factor - LAM_6 * x[6]
+
+    dx[7]  = (P0 * n - HA * (x[7] - x[8])) / (M_F * C_F)
+    dx[8]  = (HA * (x[7] - x[8]) - 2 * W_M * C_C * (x[8] - x[9])) / (M_C * C_C)
+    dx[9]  = (2 * W_M * C_C * (x[8] - x[9])) / (M_SG * C_C)
+
+    dx[10] = one(x[1])
+    return nothing
+end
+
+# --- Config loader ---
+function load_config(config_path)
+    if isfile(config_path)
+        return TOML.parsefile(config_path)
+    else
+        error("Config file not found: $config_path")
+    end
+end
+
+function build_entry_box(config)
+    if get(config, "use_predicates_entry", false)
+        pred_path = joinpath(@__DIR__, "..", "..", "..", "reachability", "predicates.json")
+        pred_raw  = JSON.parsefile(pred_path)
+        entry     = pred_raw["mode_boundaries"]["q_heatup"]["X_entry_polytope"]
+        n_lo, n_hi           = entry["n_range"]
+        T_f_lo, T_f_hi       = entry["T_f_range_C"]
+        T_c_lo, T_c_hi       = entry["T_c_range_C"]
+        T_cold_lo, T_cold_hi = entry["T_cold_range_C"]
+    else
+        e = config["entry"]
+        n_lo, n_hi           = e["n_range"]
+        T_f_lo, T_f_hi       = e["T_f_range_C"]
+        T_c_lo, T_c_hi       = e["T_c_range_C"]
+        T_cold_lo, T_cold_hi = e["T_cold_range_C"]
+    end
+
+    n_mid = 0.5 * (n_lo + n_hi)
+    C_mid = [BETA_1/(LAM_1*LAMBDA), BETA_2/(LAM_2*LAMBDA),
+             BETA_3/(LAM_3*LAMBDA), BETA_4/(LAM_4*LAMBDA),
+             BETA_5/(LAM_5*LAMBDA), BETA_6/(LAM_6*LAMBDA)] .* n_mid
+
+    x_lo = [C_mid[1]*(n_lo/n_mid), C_mid[2]*(n_lo/n_mid),
+            C_mid[3]*(n_lo/n_mid), C_mid[4]*(n_lo/n_mid),
+            C_mid[5]*(n_lo/n_mid), C_mid[6]*(n_lo/n_mid),
+            T_f_lo, T_c_lo, T_cold_lo, 0.0]
+    x_hi = [C_mid[1]*(n_hi/n_mid), C_mid[2]*(n_hi/n_mid),
+            C_mid[3]*(n_hi/n_mid), C_mid[4]*(n_hi/n_mid),
+            C_mid[5]*(n_hi/n_mid), C_mid[6]*(n_hi/n_mid),
+            T_f_hi, T_c_hi, T_cold_hi, 0.0]
+
+    return Hyperrectangle(low=x_lo, high=x_hi),
+           (n_lo=n_lo, n_hi=n_hi,
+            T_f_lo=T_f_lo, T_f_hi=T_f_hi,
+            T_c_lo=T_c_lo, T_c_hi=T_c_hi,
+            T_cold_lo=T_cold_lo, T_cold_hi=T_cold_hi)
+end
+
+# --- Per-step envelope extraction ---
+function extract_envelopes(flow_hr)
+    n_steps = length(flow_hr)
+    t_arr     = zeros(n_steps)
+    Tc_lo_ts  = zeros(n_steps); Tc_hi_ts  = zeros(n_steps)
+    Tf_lo_ts  = zeros(n_steps); Tf_hi_ts  = zeros(n_steps)
+    Tco_lo_ts = zeros(n_steps); Tco_hi_ts = zeros(n_steps)
+    n_lo_ts   = zeros(n_steps); n_hi_ts   = zeros(n_steps)
+    rho_lo_ts = zeros(n_steps); rho_hi_ts = zeros(n_steps)
+    for (k, R) in enumerate(flow_hr)
+        s = set(R)
+        t_arr[k]    = high(s, 10)
+        Tc_lo_ts[k] = low(s, 8);  Tc_hi_ts[k] = high(s, 8)
+        Tf_lo_ts[k] = low(s, 7);  Tf_hi_ts[k] = high(s, 7)
+        Tco_lo_ts[k] = low(s, 9); Tco_hi_ts[k] = high(s, 9)
+        sumLC_lo = LAM_1*low(s,1) + LAM_2*low(s,2) + LAM_3*low(s,3) +
+                   LAM_4*low(s,4) + LAM_5*low(s,5) + LAM_6*low(s,6)
+        sumLC_hi = LAM_1*high(s,1) + LAM_2*high(s,2) + LAM_3*high(s,3) +
+                   LAM_4*high(s,4) + LAM_5*high(s,5) + LAM_6*high(s,6)
+        t_hi_here = high(s, 10)
+        t_lo_here = low(s, 10)
+        Tref_lo = T_STANDBY + RAMP_RATE_CS * t_lo_here
+        Tref_hi = T_STANDBY + RAMP_RATE_CS * t_hi_here
+        rho_lo_here = KP_HEATUP * (Tref_lo - high(s, 8))
+        rho_hi_here = KP_HEATUP * (Tref_hi - low(s, 8))
+        rho_lo_ts[k] = rho_lo_here
+        rho_hi_ts[k] = rho_hi_here
+        denom_lo = BETA - rho_hi_here
+        denom_hi = BETA - rho_lo_here
+        if denom_lo > 0
+            n_lo_ts[k] = LAMBDA * sumLC_lo / denom_hi
+            n_hi_ts[k] = LAMBDA * sumLC_hi / denom_lo
+        end
+    end
+    return (; t_arr,
+              Tc_lo_ts, Tc_hi_ts, Tf_lo_ts, Tf_hi_ts,
+              Tco_lo_ts, Tco_hi_ts, n_lo_ts, n_hi_ts,
+              rho_lo_ts, rho_hi_ts)
+end
+
+# --- Main ---
+function main()
+    default_config = joinpath(@__DIR__, "..", "..", "configs", "heatup", "baseline.toml")
+    config_path = length(ARGS) > 0 ? ARGS[1] : default_config
+    # Allow a config path relative to repo root or code/.
+    if !isfile(config_path)
+        alt = joinpath(@__DIR__, "..", "..", config_path)
+        isfile(alt) && (config_path = alt)
+    end
+    config = load_config(config_path)
+
+    println("\n=== Heatup PJ reach — config: $(config["name"]) ===")
+    println("  $(get(config, "description", ""))")
+
+    X0, entry_info = build_entry_box(config)
+    println("  X_entry: n ∈ [$(entry_info.n_lo), $(entry_info.n_hi)], " *
+            "T_c ∈ [$(entry_info.T_c_lo), $(entry_info.T_c_hi)] °C")
+
+    tmjets_cfg = config["tmjets"]
+    probes = config["probes"]["horizons_seconds"]
+
+    results = Dict{Float64, Any}()
+
+    for T_probe in probes
+        println("\n--- Probe T = $T_probe s ($(round(T_probe/60; digits=1)) min) ---")
+        sys  = BlackBoxContinuousSystem(rhs_heatup_pj_taylor!, 10)
+        prob = InitialValueProblem(sys, X0)
+        try
+            alg = TMJets(
+                orderT   = tmjets_cfg["orderT"],
+                orderQ   = tmjets_cfg["orderQ"],
+                abstol   = tmjets_cfg["abstol"],
+                maxsteps = tmjets_cfg["maxsteps"],
+            )
+            t_start = time()
+            sol = solve(prob; T=Float64(T_probe), alg=alg)
+            elapsed = time() - t_start
+            flow = flowpipe(sol)
+            n_sets = length(flow)
+            println("  TMJets: $n_sets reach-sets, wall $(round(elapsed; digits=1)) s")
+
+            flow_hr = overapproximate(flow, Hyperrectangle)
+            env = extract_envelopes(flow_hr)
+
+            println("  n envelope:   [$(round(minimum(env.n_lo_ts); sigdigits=4)), $(round(maximum(env.n_hi_ts); sigdigits=4))]")
+            println("  T_c envelope: [$(round(minimum(env.Tc_lo_ts); digits=2)), $(round(maximum(env.Tc_hi_ts); digits=2))] °C")
+            println("  T_f envelope: [$(round(minimum(env.Tf_lo_ts); digits=2)), $(round(maximum(env.Tf_hi_ts); digits=2))] °C")
+            println("  T_cold env:   [$(round(minimum(env.Tco_lo_ts); digits=2)), $(round(maximum(env.Tco_hi_ts); digits=2))] °C")
+            println("  rho env:      [$(round(minimum(env.rho_lo_ts); sigdigits=4)), $(round(maximum(env.rho_hi_ts); sigdigits=4))]")
+
+            results[T_probe] = (status="OK", n_sets=n_sets, elapsed=elapsed, env=env)
+        catch err
+            msg = sprint(showerror, err)
+            println("  FAILED: ", first(msg, 300))
+            results[T_probe] = (status="FAILED", err=first(msg, 300))
+            break
+        end
+    end
+
+    println("\n=== Summary ===")
+    for T_probe in probes
+        haskey(results, T_probe) || continue
+        r = results[T_probe]
+        if r.status == "OK"
+            println("  T = $(T_probe) s: OK, $(r.n_sets) sets, $(round(r.elapsed; digits=1))s wall")
+        else
+            println("  T = $(T_probe) s: FAILED")
+        end
+    end
+
+    # --- Save ---
+    if get(config, "output", Dict()) |> (o -> get(o, "save_per_step", false))
+        result_file = config["output"]["result_file"]
+        mat_out = joinpath(@__DIR__, "..", "..", "..", "results", result_file)
+        saved = Dict{String, Any}(
+            "config_name"   => config["name"],
+            "probe_horizons" => collect(probes),
+            "beta"          => BETA,
+            "Kp"            => KP_HEATUP,
+            "T_c0"          => T_C0,
+            "T_cold0"       => T_COLD0,
+            "T_standby"     => T_STANDBY,
+        )
+        for T_probe in probes
+            haskey(results, T_probe) || continue
+            r = results[T_probe]
+            r.status == "OK" || continue
+            pre = "T_$(Int(T_probe))_"
+            env = r.env
+            saved[pre * "t_arr"]     = env.t_arr
+            saved[pre * "Tc_lo_ts"]  = env.Tc_lo_ts;  saved[pre * "Tc_hi_ts"]  = env.Tc_hi_ts
+            saved[pre * "Tf_lo_ts"]  = env.Tf_lo_ts;  saved[pre * "Tf_hi_ts"]  = env.Tf_hi_ts
+            saved[pre * "Tco_lo_ts"] = env.Tco_lo_ts; saved[pre * "Tco_hi_ts"] = env.Tco_hi_ts
+            saved[pre * "n_lo_ts"]   = env.n_lo_ts;   saved[pre * "n_hi_ts"]   = env.n_hi_ts
+            saved[pre * "rho_lo_ts"] = env.rho_lo_ts; saved[pre * "rho_hi_ts"] = env.rho_hi_ts
+        end
+        matwrite(mat_out, saved)
+        println("\nSaved per-step envelopes to $mat_out")
+    end
+end
+
+main()
diff --git a/code/scripts/reach/reach_heatup_pj_sd.jl b/code/scripts/reach/reach_heatup_pj_sd.jl
new file mode 100644
index 0000000..01cb0ee
--- /dev/null
+++ b/code/scripts/reach/reach_heatup_pj_sd.jl
@@ -0,0 +1,178 @@
+#!/usr/bin/env julia
+#
+# reach_heatup_pj_sd.jl — heatup PJ reach WITH bounded Q_sg steam dump.
+#
+# Adds an augmented state x[11] = Q_sg as a bounded parameter (dx[11] = 0).
+# The reach propagates the entry-box range of Q_sg forward; at each reach
+# set, the Q_sg extent is the disturbance envelope the controller has to
+# reject.  Physical story: operator controls secondary steam dump; actual
+# value unknown but bounded in [0, 0.05·P_0].
+#
+# Time is carried as x[12] = t (dt/dt = 1).
+#
+# Diagnostic script — uses the tight-entry + steam-dump config.
+
+using Pkg
+Pkg.activate(joinpath(@__DIR__, "..", ".."))
+
+using LinearAlgebra
+using ReachabilityAnalysis, LazySets
+using JSON
+using MAT
+using TOML
+
+# Plant constants.
+const LAMBDA  = 1e-4
+const BETA_1, BETA_2, BETA_3, BETA_4, BETA_5, BETA_6 =
+    0.000215, 0.001424, 0.001274, 0.002568, 0.000748, 0.000273
+const BETA    = BETA_1 + BETA_2 + BETA_3 + BETA_4 + BETA_5 + BETA_6
+const LAM_1, LAM_2, LAM_3, LAM_4, LAM_5, LAM_6 =
+    0.0124, 0.0305, 0.111, 0.301, 1.14, 3.01
+
+const P0   = 1e9
+const M_F, C_F, M_C, C_C, HA, W_M, M_SG =
+    50000.0, 300.0, 20000.0, 5450.0, 5e7, 5000.0, 30000.0
+
+const T_COLD0 = 290.0
+const DT_CORE = P0 / (W_M * C_C)
+const T_HOT0  = T_COLD0 + DT_CORE
+const T_C0    = (T_HOT0 + T_COLD0) / 2
+const T_F0    = T_C0 + P0 / HA
+
+const T_STANDBY    = T_C0 - 33.333333
+const RAMP_RATE_CS = 28.0 / 3600
+const KP_HEATUP    = 1e-4
+
+# 12-state RHS: [C_1..C_6, T_f, T_c, T_cold, Q_sg, t]
+# dx[10] = 0   (Q_sg is a bounded parameter)
+# dx[11] = 1   (time)
+# Hmm — indexing conflict with original version.  Reorganize:
+#   x[1..6]  = C_1..C_6
+#   x[7]     = T_f
+#   x[8]     = T_c
+#   x[9]     = T_cold
+#   x[10]    = Q_sg (constant-parameter, dx[10]=0)
+#   x[11]    = t (dx[11]=1)
+@taylorize function rhs_heatup_sd_taylor!(dx, x, p, t)
+    rho = KP_HEATUP * (T_STANDBY + RAMP_RATE_CS * x[11] - x[8])
+    sum_lam_C = LAM_1*x[1] + LAM_2*x[2] + LAM_3*x[3] +
+                LAM_4*x[4] + LAM_5*x[5] + LAM_6*x[6]
+    denom = BETA - rho
+    n = LAMBDA * sum_lam_C / denom
+    inv_factor = sum_lam_C / denom
+
+    dx[1] = BETA_1 * inv_factor - LAM_1 * x[1]
+    dx[2] = BETA_2 * inv_factor - LAM_2 * x[2]
+    dx[3] = BETA_3 * inv_factor - LAM_3 * x[3]
+    dx[4] = BETA_4 * inv_factor - LAM_4 * x[4]
+    dx[5] = BETA_5 * inv_factor - LAM_5 * x[5]
+    dx[6] = BETA_6 * inv_factor - LAM_6 * x[6]
+
+    dx[7]  = (P0 * n - HA * (x[7] - x[8])) / (M_F * C_F)
+    dx[8]  = (HA * (x[7] - x[8]) - 2 * W_M * C_C * (x[8] - x[9])) / (M_C * C_C)
+    dx[9]  = (2 * W_M * C_C * (x[8] - x[9]) - x[10]) / (M_SG * C_C)  # Q_sg now enters!
+
+    dx[10] = zero(x[1])   # Q_sg constant
+    dx[11] = one(x[1])    # time
+
+    return nothing
+end
+
+function main()
+    config_path = length(ARGS) > 0 ? ARGS[1] :
+        joinpath(@__DIR__, "..", "..", "configs", "heatup", "with_steam_dump.toml")
+    if !isfile(config_path)
+        alt = joinpath(@__DIR__, "..", "..", config_path)
+        isfile(alt) && (config_path = alt)
+    end
+    config = TOML.parsefile(config_path)
+
+    println("\n=== Heatup PJ reach with steam dump — config: $(config["name"]) ===")
+    println("  $(get(config, "description", ""))")
+
+    e = config["entry"]
+    n_lo, n_hi           = e["n_range"]
+    T_f_lo, T_f_hi       = e["T_f_range_C"]
+    T_c_lo, T_c_hi       = e["T_c_range_C"]
+    T_cold_lo, T_cold_hi = e["T_cold_range_C"]
+    sd = e["steam_dump"]
+    Q_lo = sd["Q_lo_fraction_P0"] * P0
+    Q_hi = sd["Q_hi_fraction_P0"] * P0
+
+    n_mid = 0.5 * (n_lo + n_hi)
+    C_mid = [BETA_1/(LAM_1*LAMBDA), BETA_2/(LAM_2*LAMBDA),
+             BETA_3/(LAM_3*LAMBDA), BETA_4/(LAM_4*LAMBDA),
+             BETA_5/(LAM_5*LAMBDA), BETA_6/(LAM_6*LAMBDA)] .* n_mid
+
+    x_lo = [C_mid[1]*(n_lo/n_mid), C_mid[2]*(n_lo/n_mid),
+            C_mid[3]*(n_lo/n_mid), C_mid[4]*(n_lo/n_mid),
+            C_mid[5]*(n_lo/n_mid), C_mid[6]*(n_lo/n_mid),
+            T_f_lo, T_c_lo, T_cold_lo, Q_lo, 0.0]
+    x_hi = [C_mid[1]*(n_hi/n_mid), C_mid[2]*(n_hi/n_mid),
+            C_mid[3]*(n_hi/n_mid), C_mid[4]*(n_hi/n_mid),
+            C_mid[5]*(n_hi/n_mid), C_mid[6]*(n_hi/n_mid),
+            T_f_hi, T_c_hi, T_cold_hi, Q_hi, 0.0]
+
+    X0 = Hyperrectangle(low=x_lo, high=x_hi)
+    println("  n ∈ [$n_lo, $n_hi], T_c ∈ [$T_c_lo, $T_c_hi] °C")
+    println("  Q_sg (steam dump) ∈ [$(Q_lo/1e6), $(Q_hi/1e6)] MW")
+
+    t_cfg = config["tmjets"]
+    probes = config["probes"]["horizons_seconds"]
+
+    results = Dict{Float64, Any}()
+    for T_probe in probes
+        println("\n--- Probe T = $T_probe s ---")
+        sys  = BlackBoxContinuousSystem(rhs_heatup_sd_taylor!, 11)
+        prob = InitialValueProblem(sys, X0)
+        try
+            alg = TMJets(orderT=t_cfg["orderT"], orderQ=t_cfg["orderQ"],
+                         abstol=t_cfg["abstol"], maxsteps=t_cfg["maxsteps"])
+            t0 = time()
+            sol = solve(prob; T=Float64(T_probe), alg=alg)
+            elapsed = time() - t0
+            flow_hr = overapproximate(flowpipe(sol), Hyperrectangle)
+            n_sets = length(flow_hr)
+            println("  TMJets: $n_sets reach-sets in $(round(elapsed; digits=1)) s")
+
+            Tc_lo_env = minimum(low(set(R), 8) for R in flow_hr)
+            Tc_hi_env = maximum(high(set(R), 8) for R in flow_hr)
+            Tco_lo_env = minimum(low(set(R), 9) for R in flow_hr)
+            Tco_hi_env = maximum(high(set(R), 9) for R in flow_hr)
+            Tf_lo_env = minimum(low(set(R), 7) for R in flow_hr)
+            Tf_hi_env = maximum(high(set(R), 7) for R in flow_hr)
+            println("  T_c envelope:   [$(round(Tc_lo_env; digits=2)), $(round(Tc_hi_env; digits=2))] °C")
+            println("  T_f envelope:   [$(round(Tf_lo_env; digits=2)), $(round(Tf_hi_env; digits=2))] °C")
+            println("  T_cold env:     [$(round(Tco_lo_env; digits=2)), $(round(Tco_hi_env; digits=2))] °C")
+
+            # Compare low-T_avg trip (280 °C).
+            low_trip_ok = Tc_lo_env >= 280.0
+            println("  Low-T_avg trip (T_c ≥ 280): $(low_trip_ok ? "✅" : "× loose")")
+
+            results[T_probe] = (status="OK", Tc=(Tc_lo_env, Tc_hi_env),
+                                Tf=(Tf_lo_env, Tf_hi_env),
+                                Tcold=(Tco_lo_env, Tco_hi_env),
+                                low_trip_ok=low_trip_ok)
+        catch err
+            println("  FAILED: ", first(sprint(showerror, err), 300))
+            results[T_probe] = (status="FAILED",)
+            break
+        end
+    end
+
+    mat_out = joinpath(@__DIR__, "..", "..", "..", "results", config["output"]["result_file"])
+    saved = Dict{String, Any}("config_name" => config["name"],
+                              "Q_lo" => Q_lo, "Q_hi" => Q_hi)
+    for (T_probe, r) in results
+        if r.status == "OK"
+            pre = "T_$(Int(T_probe))_"
+            saved[pre * "Tc_lo"] = r.Tc[1]; saved[pre * "Tc_hi"] = r.Tc[2]
+            saved[pre * "Tf_lo"] = r.Tf[1]; saved[pre * "Tf_hi"] = r.Tf[2]
+            saved[pre * "Tcold_lo"] = r.Tcold[1]; saved[pre * "Tcold_hi"] = r.Tcold[2]
+        end
+    end
+    matwrite(mat_out, saved)
+    println("\nSaved to $mat_out")
+end
+
+main()
diff --git a/code/scripts/reach_operation.jl b/code/scripts/reach/reach_operation.jl
similarity index 91%
rename from code/scripts/reach_operation.jl
rename to code/scripts/reach/reach_operation.jl
index 16640e4..73581d8 100644
--- a/code/scripts/reach_operation.jl
+++ b/code/scripts/reach/reach_operation.jl
@@ -11,7 +11,7 @@
 # Linear-reach results here are an approximate-but-useful gut check.
 
 using Pkg
-Pkg.activate(joinpath(@__DIR__, ".."))
+Pkg.activate(joinpath(@__DIR__, "..", ".."))
 
 using Printf
 using LinearAlgebra
@@ -20,11 +20,11 @@ using Plots
 using JSON
 using MAT
 
-include(joinpath(@__DIR__, "..", "src", "pke_params.jl"))
-include(joinpath(@__DIR__, "..", "src", "pke_th_rhs.jl"))
-include(joinpath(@__DIR__, "..", "src", "pke_linearize.jl"))
-include(joinpath(@__DIR__, "..", "src", "load_predicates.jl"))
-include(joinpath(@__DIR__, "..", "src", "reach_linear.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_params.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_th_rhs.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_linearize.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "load_predicates.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "reach_linear.jl"))
 
 plant = pke_params()
 x_op  = pke_initial_conditions(plant)
@@ -101,7 +101,7 @@ end
 # --- Plots: T_c reach tube, two views ---
 CtoF(T) = T * 9/5 + 32
 delta_safe_Tc = pred.constants.t_avg_in_range_halfwidth_C
-figdir = joinpath(@__DIR__, "..", "..", "docs", "figures")
+figdir = joinpath(@__DIR__, "..", "..", "..", "docs", "figures")
 isdir(figdir) || mkpath(figdir)
 
 p_safety = plot(T, CtoF.(X_nom[9, :]), lw=1.2, color=:red,
@@ -140,7 +140,7 @@ plot!(fig_n, T, X_hi[1, :], fillrange=X_lo[1, :], fillalpha=0.3,
 savefig(fig_n, joinpath(figdir, "reach_operation_n.png"))
 
 # --- Save result ---
-matfile = joinpath(@__DIR__, "..", "..", "reachability", "reach_operation_result.mat")
+matfile = joinpath(@__DIR__, "..", "..", "..", "results", "reach_operation_result.mat")
 matwrite(matfile, Dict("T" => collect(T), "R_lo" => R_lo, "R_hi" => R_hi,
                        "center" => Cnom, "X_lo" => X_lo, "X_hi" => X_hi,
                        "X_nom" => X_nom, "K" => Matrix(K), "A_cl" => A_cl,
diff --git a/code/scripts/reach_scram_pj.jl b/code/scripts/reach/reach_scram_pj.jl
similarity index 98%
rename from code/scripts/reach_scram_pj.jl
rename to code/scripts/reach/reach_scram_pj.jl
index d51aa65..9b3d4f4 100644
--- a/code/scripts/reach_scram_pj.jl
+++ b/code/scripts/reach/reach_scram_pj.jl
@@ -9,7 +9,7 @@
 # 9-state PJ model (10D with augmented time).
 
 using Pkg
-Pkg.activate(joinpath(@__DIR__, ".."))
+Pkg.activate(joinpath(@__DIR__, "..", ".."))
 
 using LinearAlgebra
 using ReachabilityAnalysis, LazySets
@@ -140,7 +140,7 @@ for T_probe in (10.0, 30.0, 60.0)
     end
 end
 
-mat_out = joinpath(@__DIR__, "..", "..", "reachability", "reach_scram_pj_result.mat")
+mat_out = joinpath(@__DIR__, "..", "..", "..", "results", "reach_scram_pj_result.mat")
 saved = Dict{String, Any}("probe_horizons" => collect((10.0, 30.0, 60.0)))
 for T_probe in (10.0, 30.0, 60.0)
     haskey(results, T_probe) || continue
diff --git a/code/scripts/reach_heatup_pj.jl b/code/scripts/reach_heatup_pj.jl
deleted file mode 100644
index 5ddd4bb..0000000
--- a/code/scripts/reach_heatup_pj.jl
+++ /dev/null
@@ -1,243 +0,0 @@
-#!/usr/bin/env julia
-#
-# reach_heatup_pj.jl — nonlinear reach on heatup, prompt-jump model.
-#
-# Reduced from 10-state to 9-state (n is algebraic).  Removes the
-# Λ⁻¹ stiffness that capped the full-state reach at ~10 s.  We push
-# horizons up: 60 s, 300 s, 1800 s, 5400 s, full T_max = 18000 s (5 hr).
-#
-# State (10D with augmented time):
-#   x[1..6] = C_1..C_6   (delayed-neutron precursors)
-#   x[7]    = T_f
-#   x[8]    = T_c
-#   x[9]    = T_cold
-#   x[10]   = t (augmented time, dt/dt = 1)
-#
-# n is algebraic:  n = Λ · Σ λ_i C_i / (β - ρ),  ρ = K_p · (T_ref - T_c).
-
-using Pkg
-Pkg.activate(joinpath(@__DIR__, ".."))
-
-using LinearAlgebra
-using ReachabilityAnalysis, LazySets
-using JSON
-using MAT
-
-# --- Inlined plant constants (must match pke_params) ---
-const LAMBDA  = 1e-4
-const BETA_1, BETA_2, BETA_3, BETA_4, BETA_5, BETA_6 =
-    0.000215, 0.001424, 0.001274, 0.002568, 0.000748, 0.000273
-const BETA    = BETA_1 + BETA_2 + BETA_3 + BETA_4 + BETA_5 + BETA_6
-const LAM_1, LAM_2, LAM_3, LAM_4, LAM_5, LAM_6 =
-    0.0124, 0.0305, 0.111, 0.301, 1.14, 3.01
-
-const P0   = 1e9
-const M_F, C_F, M_C, C_C, HA, W_M, M_SG =
-    50000.0, 300.0, 20000.0, 5450.0, 5e7, 5000.0, 30000.0
-
-# Note: feedback-linearization in ctrl_heatup_unsat cancels the alpha
-# terms exactly, so under closed-loop the effective rho is just Kp·e.
-# We don't need ALPHA_F, ALPHA_C in the reach RHS as a result.
-
-const T_COLD0 = 290.0
-const DT_CORE = P0 / (W_M * C_C)
-const T_HOT0  = T_COLD0 + DT_CORE
-const T_C0    = (T_HOT0 + T_COLD0) / 2
-const T_F0    = T_C0 + P0 / HA
-
-const T_STANDBY    = T_C0 - 33.333333
-const RAMP_RATE_CS = 28.0 / 3600
-const KP_HEATUP    = 1e-4
-
-# --- Taylorized PJ heatup RHS, 10D with augmented time ---
-@taylorize function rhs_heatup_pj_taylor!(dx, x, p, t)
-    # rho_total under closed-loop feedback linearization = Kp · (T_ref - T_c).
-    rho = KP_HEATUP * (T_STANDBY + RAMP_RATE_CS * x[10] - x[8])
-
-    # Algebraic prompt-jump n.
-    sum_lam_C = LAM_1*x[1] + LAM_2*x[2] + LAM_3*x[3] +
-                LAM_4*x[4] + LAM_5*x[5] + LAM_6*x[6]
-    denom = BETA - rho
-    n = LAMBDA * sum_lam_C / denom
-    inv_factor = sum_lam_C / denom
-
-    # Precursor balance under PJ.
-    dx[1] = BETA_1 * inv_factor - LAM_1 * x[1]
-    dx[2] = BETA_2 * inv_factor - LAM_2 * x[2]
-    dx[3] = BETA_3 * inv_factor - LAM_3 * x[3]
-    dx[4] = BETA_4 * inv_factor - LAM_4 * x[4]
-    dx[5] = BETA_5 * inv_factor - LAM_5 * x[5]
-    dx[6] = BETA_6 * inv_factor - LAM_6 * x[6]
-
-    # Thermal — n appears algebraically in fuel eq.
-    dx[7] = (P0 * n - HA * (x[7] - x[8])) / (M_F * C_F)
-    dx[8] = (HA * (x[7] - x[8]) - 2 * W_M * C_C * (x[8] - x[9])) / (M_C * C_C)
-    dx[9] = (2 * W_M * C_C * (x[8] - x[9])) / (M_SG * C_C)
-
-    # Augmented time.
-    dx[10] = one(x[1])
-
-    return nothing
-end
-
-# --- Build X_entry (PJ form: no n) from predicates.json ---
-pred_path = joinpath(@__DIR__, "..", "..", "reachability", "predicates.json")
-pred_raw  = JSON.parsefile(pred_path)
-entry     = pred_raw["mode_boundaries"]["q_heatup"]["X_entry_polytope"]
-
-n_lo, n_hi           = entry["n_range"]
-T_f_lo, T_f_hi       = entry["T_f_range_C"]
-T_c_lo, T_c_hi       = entry["T_c_range_C"]
-T_cold_lo, T_cold_hi = entry["T_cold_range_C"]
-
-n_mid = 0.5 * (n_lo + n_hi)
-C_mid_1 = (BETA_1 / (LAM_1 * LAMBDA)) * n_mid
-C_mid_2 = (BETA_2 / (LAM_2 * LAMBDA)) * n_mid
-C_mid_3 = (BETA_3 / (LAM_3 * LAMBDA)) * n_mid
-C_mid_4 = (BETA_4 / (LAM_4 * LAMBDA)) * n_mid
-C_mid_5 = (BETA_5 / (LAM_5 * LAMBDA)) * n_mid
-C_mid_6 = (BETA_6 / (LAM_6 * LAMBDA)) * n_mid
-
-# C_i scale linearly with n; sweep across the n_lo..n_hi band.
-x_lo = [C_mid_1 * (n_lo / n_mid); C_mid_2 * (n_lo / n_mid);
-        C_mid_3 * (n_lo / n_mid); C_mid_4 * (n_lo / n_mid);
-        C_mid_5 * (n_lo / n_mid); C_mid_6 * (n_lo / n_mid);
-        T_f_lo; T_c_lo; T_cold_lo;
-        0.0]
-x_hi = [C_mid_1 * (n_hi / n_mid); C_mid_2 * (n_hi / n_mid);
-        C_mid_3 * (n_hi / n_mid); C_mid_4 * (n_hi / n_mid);
-        C_mid_5 * (n_hi / n_mid); C_mid_6 * (n_hi / n_mid);
-        T_f_hi; T_c_hi; T_cold_hi;
-        0.0]
-
-X0 = Hyperrectangle(low=x_lo, high=x_hi)
-
-println("\n=== Nonlinear heatup reach, prompt-jump model ===")
-println("  X_entry (n-implied range): n ∈ [$(n_lo), $(n_hi)]")
-println("  T_c ∈ [$(T_c_lo), $(T_c_hi)] °C")
-
-T_RAMP_END = (T_C0 - T_STANDBY) / RAMP_RATE_CS
-println("  T_ramp_end = $(round(T_RAMP_END; digits=0)) s ($(round(T_RAMP_END/60; digits=1)) min)")
-println("  Probing horizons up to T_max(heatup) = 18000 s (5 hr).")
-
-# Probe at increasing horizons.  Stop early if any probe fails.
-probe_horizons = (60.0, 300.0, 1800.0, 5400.0)
-
-results = Dict{Float64, Any}()
-for T_probe in probe_horizons
-    println("\n--- Probe T = $T_probe s ($(round(T_probe/60; digits=1)) min) ---")
-    sys  = BlackBoxContinuousSystem(rhs_heatup_pj_taylor!, 10)
-    prob = InitialValueProblem(sys, X0)
-    try
-        alg = TMJets(orderT=4, orderQ=2, abstol=1e-9, maxsteps=100000)
-        t_start = time()
-        sol = solve(prob; T=T_probe, alg=alg)
-        elapsed = time() - t_start
-        flow = flowpipe(sol)
-        n_sets = length(flow)
-        println("  TMJets: $n_sets reach-sets, wall-time $(round(elapsed; digits=1)) s")
-
-        flow_hr = overapproximate(flow, Hyperrectangle)
-        n_steps = length(flow_hr)
-        # Per-timestep envelopes for plotting (saved below).
-        t_arr     = zeros(n_steps)
-        Tc_lo_ts  = zeros(n_steps); Tc_hi_ts  = zeros(n_steps)
-        Tf_lo_ts  = zeros(n_steps); Tf_hi_ts  = zeros(n_steps)
-        Tco_lo_ts = zeros(n_steps); Tco_hi_ts = zeros(n_steps)
-        n_lo_ts   = zeros(n_steps); n_hi_ts   = zeros(n_steps)
-        rho_lo_ts = zeros(n_steps); rho_hi_ts = zeros(n_steps)
-        for (k, R) in enumerate(flow_hr)
-            s = set(R)
-            t_arr[k]    = high(s, 10)  # time-state upper bound ≈ current time
-            Tc_lo_ts[k] = low(s, 8);  Tc_hi_ts[k] = high(s, 8)
-            Tf_lo_ts[k] = low(s, 7);  Tf_hi_ts[k] = high(s, 7)
-            Tco_lo_ts[k] = low(s, 9); Tco_hi_ts[k] = high(s, 9)
-            # Algebraic n reconstruction at this reach-set.
-            sumLC_lo = LAM_1*low(s,1) + LAM_2*low(s,2) + LAM_3*low(s,3) +
-                       LAM_4*low(s,4) + LAM_5*low(s,5) + LAM_6*low(s,6)
-            sumLC_hi = LAM_1*high(s,1) + LAM_2*high(s,2) + LAM_3*high(s,3) +
-                       LAM_4*high(s,4) + LAM_5*high(s,5) + LAM_6*high(s,6)
-            rho_lo_here = KP_HEATUP * (T_STANDBY - high(s, 8))
-            rho_hi_here = KP_HEATUP * (T_C0 - low(s, 8))
-            rho_lo_ts[k] = rho_lo_here
-            rho_hi_ts[k] = rho_hi_here
-            denom_lo = BETA - rho_hi_here
-            denom_hi = BETA - rho_lo_here
-            if denom_lo > 0
-                n_lo_ts[k] = LAMBDA * sumLC_lo / denom_hi
-                n_hi_ts[k] = LAMBDA * sumLC_hi / denom_lo
-            end
-        end
-        # Also global-envelope values (for backward compat).
-        Tc_lo_env = minimum(Tc_lo_ts); Tc_hi_env = maximum(Tc_hi_ts)
-        Tf_lo_env = minimum(Tf_lo_ts); Tf_hi_env = maximum(Tf_hi_ts)
-        Tcold_lo  = minimum(Tco_lo_ts); Tcold_hi = maximum(Tco_hi_ts)
-        n_env_lo = minimum(n_lo_ts); n_env_hi = maximum(n_hi_ts)
-
-        println("  n envelope (reconstructed): [$(round(n_env_lo; sigdigits=4)), $(round(n_env_hi; sigdigits=4))]")
-        println("  T_f envelope:               [$(round(Tf_lo_env; digits=2)), $(round(Tf_hi_env; digits=2))] °C")
-        println("  T_c envelope:               [$(round(Tc_lo_env; digits=2)), $(round(Tc_hi_env; digits=2))] °C")
-        println("  T_cold envelope:            [$(round(Tcold_lo; digits=2)), $(round(Tcold_hi; digits=2))] °C")
-        results[T_probe] = (status="OK", n_sets=n_sets, elapsed=elapsed,
-                            Tc=(Tc_lo_env, Tc_hi_env),
-                            Tf=(Tf_lo_env, Tf_hi_env),
-                            Tcold=(Tcold_lo, Tcold_hi),
-                            n=(n_env_lo, n_env_hi),
-                            t_arr=t_arr,
-                            Tc_lo_ts=Tc_lo_ts, Tc_hi_ts=Tc_hi_ts,
-                            Tf_lo_ts=Tf_lo_ts, Tf_hi_ts=Tf_hi_ts,
-                            Tco_lo_ts=Tco_lo_ts, Tco_hi_ts=Tco_hi_ts,
-                            n_lo_ts=n_lo_ts, n_hi_ts=n_hi_ts,
-                            rho_lo_ts=rho_lo_ts, rho_hi_ts=rho_hi_ts)
-    catch err
-        msg = sprint(showerror, err)
-        println("  FAILED: ", first(msg, 300))
-        results[T_probe] = (status="FAILED", err=first(msg, 300))
-        break
-    end
-end
-
-println("\n=== Summary ===")
-for T_probe in probe_horizons
-    haskey(results, T_probe) || continue
-    r = results[T_probe]
-    if r.status == "OK"
-        println("  T = $(T_probe) s: OK, $(r.n_sets) reach-sets, $(round(r.elapsed; digits=1))s wall")
-    else
-        println("  T = $(T_probe) s: FAILED")
-    end
-end
-
-# Save the longest-successful probe's envelope arrays for the app.
-mat_out = joinpath(@__DIR__, "..", "..", "reachability", "reach_heatup_pj_result.mat")
-saved = Dict{String, Any}()
-saved["probe_horizons"] = collect(probe_horizons)
-for T_probe in probe_horizons
-    haskey(results, T_probe) || continue
-    r = results[T_probe]
-    if r.status == "OK"
-        pre = "T_$(Int(T_probe))_"
-        saved[pre * "Tc_lo"]    = r.Tc[1]
-        saved[pre * "Tc_hi"]    = r.Tc[2]
-        saved[pre * "Tf_lo"]    = r.Tf[1]
-        saved[pre * "Tf_hi"]    = r.Tf[2]
-        saved[pre * "Tcold_lo"] = r.Tcold[1]
-        saved[pre * "Tcold_hi"] = r.Tcold[2]
-        saved[pre * "n_lo"]     = r.n[1]
-        saved[pre * "n_hi"]     = r.n[2]
-        # Per-timestep envelopes
-        saved[pre * "t_arr"]     = r.t_arr
-        saved[pre * "Tc_lo_ts"]  = r.Tc_lo_ts
-        saved[pre * "Tc_hi_ts"]  = r.Tc_hi_ts
-        saved[pre * "Tf_lo_ts"]  = r.Tf_lo_ts
-        saved[pre * "Tf_hi_ts"]  = r.Tf_hi_ts
-        saved[pre * "Tco_lo_ts"] = r.Tco_lo_ts
-        saved[pre * "Tco_hi_ts"] = r.Tco_hi_ts
-        saved[pre * "n_lo_ts"]   = r.n_lo_ts
-        saved[pre * "n_hi_ts"]   = r.n_hi_ts
-        saved[pre * "rho_lo_ts"] = r.rho_lo_ts
-        saved[pre * "rho_hi_ts"] = r.rho_hi_ts
-    end
-end
-matwrite(mat_out, saved)
-println("\nSaved envelope summaries to $mat_out")
diff --git a/code/scripts/reach_heatup_pj_tight.jl b/code/scripts/reach_heatup_pj_tight.jl
deleted file mode 100644
index a677f31..0000000
--- a/code/scripts/reach_heatup_pj_tight.jl
+++ /dev/null
@@ -1,109 +0,0 @@
-#!/usr/bin/env julia
-#
-# reach_heatup_pj_tight.jl — heatup PJ reach with a tighter X_entry.
-#
-# The default X_entry (from predicates.json) has T_c ∈ [281, 295] — 14 K
-# wide.  The baseline PJ reach at T=300s produced T_c envelope
-# [272.4, 295.0], violating the low-T_avg trip at 280.
-#
-# Hypothesis: entry-box width is contributing to tube growth.  Try
-# T_c ∈ [285, 291] (6 K) and T_f matched, see if the lower envelope
-# rises above 280.
-
-using Pkg
-Pkg.activate(joinpath(@__DIR__, ".."))
-
-using LinearAlgebra
-using ReachabilityAnalysis, LazySets
-using MAT
-
-# Same constants as reach_heatup_pj.jl.
-const LAMBDA  = 1e-4
-const BETA_1, BETA_2, BETA_3, BETA_4, BETA_5, BETA_6 =
-    0.000215, 0.001424, 0.001274, 0.002568, 0.000748, 0.000273
-const BETA    = BETA_1 + BETA_2 + BETA_3 + BETA_4 + BETA_5 + BETA_6
-const LAM_1, LAM_2, LAM_3, LAM_4, LAM_5, LAM_6 =
-    0.0124, 0.0305, 0.111, 0.301, 1.14, 3.01
-const P0   = 1e9
-const M_F, C_F, M_C, C_C, HA, W_M, M_SG =
-    50000.0, 300.0, 20000.0, 5450.0, 5e7, 5000.0, 30000.0
-const T_COLD0 = 290.0
-const DT_CORE = P0 / (W_M * C_C)
-const T_HOT0  = T_COLD0 + DT_CORE
-const T_C0    = (T_HOT0 + T_COLD0) / 2
-const T_F0    = T_C0 + P0 / HA
-
-const T_STANDBY    = T_C0 - 33.333333
-const RAMP_RATE_CS = 28.0 / 3600
-const KP_HEATUP    = 1e-4
-
-@taylorize function rhs_heatup_pj_tight!(dx, x, p, t)
-    rho = KP_HEATUP * (T_STANDBY + RAMP_RATE_CS * x[10] - x[8])
-    sum_lam_C = LAM_1*x[1] + LAM_2*x[2] + LAM_3*x[3] +
-                LAM_4*x[4] + LAM_5*x[5] + LAM_6*x[6]
-    denom = BETA - rho
-    n = LAMBDA * sum_lam_C / denom
-    inv_factor = sum_lam_C / denom
-    dx[1] = BETA_1 * inv_factor - LAM_1 * x[1]
-    dx[2] = BETA_2 * inv_factor - LAM_2 * x[2]
-    dx[3] = BETA_3 * inv_factor - LAM_3 * x[3]
-    dx[4] = BETA_4 * inv_factor - LAM_4 * x[4]
-    dx[5] = BETA_5 * inv_factor - LAM_5 * x[5]
-    dx[6] = BETA_6 * inv_factor - LAM_6 * x[6]
-    dx[7] = (P0 * n - HA * (x[7] - x[8])) / (M_F * C_F)
-    dx[8] = (HA * (x[7] - x[8]) - 2 * W_M * C_C * (x[8] - x[9])) / (M_C * C_C)
-    dx[9] = (2 * W_M * C_C * (x[8] - x[9])) / (M_SG * C_C)
-    dx[10] = one(x[1])
-    return nothing
-end
-
-# Tighter X_entry on T_c and T_f specifically.
-n_lo, n_hi         = 1.0e-3, 2.0e-3    # narrower n
-T_f_lo, T_f_hi     = 285.0, 291.0      # was 275–295
-T_c_lo, T_c_hi     = 285.0, 291.0      # was 281–295
-T_cold_lo, T_cold_hi = 278.0, 285.0    # was 270–281 (shifted up)
-
-n_mid = 0.5 * (n_lo + n_hi)
-C_mid = [BETA_1/(LAM_1*LAMBDA), BETA_2/(LAM_2*LAMBDA),
-         BETA_3/(LAM_3*LAMBDA), BETA_4/(LAM_4*LAMBDA),
-         BETA_5/(LAM_5*LAMBDA), BETA_6/(LAM_6*LAMBDA)] .* n_mid
-
-x_lo = [C_mid[1]*(n_lo/n_mid), C_mid[2]*(n_lo/n_mid),
-        C_mid[3]*(n_lo/n_mid), C_mid[4]*(n_lo/n_mid),
-        C_mid[5]*(n_lo/n_mid), C_mid[6]*(n_lo/n_mid),
-        T_f_lo, T_c_lo, T_cold_lo, 0.0]
-x_hi = [C_mid[1]*(n_hi/n_mid), C_mid[2]*(n_hi/n_mid),
-        C_mid[3]*(n_hi/n_mid), C_mid[4]*(n_hi/n_mid),
-        C_mid[5]*(n_hi/n_mid), C_mid[6]*(n_hi/n_mid),
-        T_f_hi, T_c_hi, T_cold_hi, 0.0]
-
-X0 = Hyperrectangle(low=x_lo, high=x_hi)
-
-println("\n=== Heatup PJ reach with TIGHTENED X_entry ===")
-println("  T_c ∈ [$(T_c_lo), $(T_c_hi)] (width 6 K, was 14 K)")
-println("  T_f ∈ [$(T_f_lo), $(T_f_hi)]")
-println("  n-implied ∈ [$(n_lo), $(n_hi)]")
-
-for T_probe in (60.0, 300.0)
-    println("\n--- Probe T = $T_probe s ---")
-    sys  = BlackBoxContinuousSystem(rhs_heatup_pj_tight!, 10)
-    prob = InitialValueProblem(sys, X0)
-    try
-        alg = TMJets(orderT=4, orderQ=2, abstol=1e-9, maxsteps=100000)
-        t_start = time()
-        sol = solve(prob; T=T_probe, alg=alg)
-        elapsed = time() - t_start
-        flow = flowpipe(sol)
-        flow_hr = overapproximate(flow, Hyperrectangle)
-        Tc_lo_env = minimum(low(set(R), 8)  for R in flow_hr)
-        Tc_hi_env = maximum(high(set(R), 8) for R in flow_hr)
-        Tf_lo_env = minimum(low(set(R), 7)  for R in flow_hr)
-        Tf_hi_env = maximum(high(set(R), 7) for R in flow_hr)
-        println("  $(length(flow_hr)) sets in $(round(elapsed; digits=1))s")
-        println("  T_c envelope: [$(round(Tc_lo_env; digits=2)), $(round(Tc_hi_env; digits=2))] °C")
-        println("  T_f envelope: [$(round(Tf_lo_env; digits=2)), $(round(Tf_hi_env; digits=2))] °C")
-        println("  Low-T_avg trip (T_c ≥ 280): $(Tc_lo_env >= 280 ? "✅ DISCHARGED" : "× still loose")")
-    catch err
-        println("  FAILED: ", first(sprint(showerror, err), 200))
-    end
-end
diff --git a/code/scripts/reach_heatup_pj_tight_full.jl b/code/scripts/reach_heatup_pj_tight_full.jl
deleted file mode 100644
index 5808efe..0000000
--- a/code/scripts/reach_heatup_pj_tight_full.jl
+++ /dev/null
@@ -1,145 +0,0 @@
-#!/usr/bin/env julia
-#
-# reach_heatup_pj_tight_full.jl — tight-entry heatup PJ reach,
-# per-timestep envelopes saved for plotting (T_c, T_h, T_cold, rho, n).
-#
-# Identical dynamics to reach_heatup_pj_tight.jl but keeps the full
-# tube data so we can overlay T_c and T_h (and infer power) on one plot.
-
-using Pkg
-Pkg.activate(joinpath(@__DIR__, ".."))
-
-using LinearAlgebra
-using ReachabilityAnalysis, LazySets
-using MAT
-
-const LAMBDA  = 1e-4
-const BETA_1, BETA_2, BETA_3, BETA_4, BETA_5, BETA_6 =
-    0.000215, 0.001424, 0.001274, 0.002568, 0.000748, 0.000273
-const BETA    = BETA_1 + BETA_2 + BETA_3 + BETA_4 + BETA_5 + BETA_6
-const LAM_1, LAM_2, LAM_3, LAM_4, LAM_5, LAM_6 =
-    0.0124, 0.0305, 0.111, 0.301, 1.14, 3.01
-const P0   = 1e9
-const M_F, C_F, M_C, C_C, HA, W_M, M_SG =
-    50000.0, 300.0, 20000.0, 5450.0, 5e7, 5000.0, 30000.0
-const T_COLD0 = 290.0
-const DT_CORE = P0 / (W_M * C_C)
-const T_HOT0  = T_COLD0 + DT_CORE
-const T_C0    = (T_HOT0 + T_COLD0) / 2
-const T_F0    = T_C0 + P0 / HA
-
-const T_STANDBY    = T_C0 - 33.333333
-const RAMP_RATE_CS = 28.0 / 3600
-const KP_HEATUP    = 1e-4
-
-@taylorize function rhs_heatup_pj_tf!(dx, x, p, t)
-    rho = KP_HEATUP * (T_STANDBY + RAMP_RATE_CS * x[10] - x[8])
-    sum_lam_C = LAM_1*x[1] + LAM_2*x[2] + LAM_3*x[3] +
-                LAM_4*x[4] + LAM_5*x[5] + LAM_6*x[6]
-    denom = BETA - rho
-    n = LAMBDA * sum_lam_C / denom
-    inv_factor = sum_lam_C / denom
-    dx[1] = BETA_1 * inv_factor - LAM_1 * x[1]
-    dx[2] = BETA_2 * inv_factor - LAM_2 * x[2]
-    dx[3] = BETA_3 * inv_factor - LAM_3 * x[3]
-    dx[4] = BETA_4 * inv_factor - LAM_4 * x[4]
-    dx[5] = BETA_5 * inv_factor - LAM_5 * x[5]
-    dx[6] = BETA_6 * inv_factor - LAM_6 * x[6]
-    dx[7] = (P0 * n - HA * (x[7] - x[8])) / (M_F * C_F)
-    dx[8] = (HA * (x[7] - x[8]) - 2 * W_M * C_C * (x[8] - x[9])) / (M_C * C_C)
-    dx[9] = (2 * W_M * C_C * (x[8] - x[9])) / (M_SG * C_C)
-    dx[10] = one(x[1])
-    return nothing
-end
-
-# Tight X_entry.
-n_lo, n_hi         = 1.0e-3, 2.0e-3
-T_f_lo, T_f_hi     = 285.0, 291.0
-T_c_lo, T_c_hi     = 285.0, 291.0
-T_cold_lo, T_cold_hi = 278.0, 285.0
-
-n_mid = 0.5 * (n_lo + n_hi)
-C_mid = [BETA_1/(LAM_1*LAMBDA), BETA_2/(LAM_2*LAMBDA),
-         BETA_3/(LAM_3*LAMBDA), BETA_4/(LAM_4*LAMBDA),
-         BETA_5/(LAM_5*LAMBDA), BETA_6/(LAM_6*LAMBDA)] .* n_mid
-
-x_lo = [C_mid[1]*(n_lo/n_mid), C_mid[2]*(n_lo/n_mid),
-        C_mid[3]*(n_lo/n_mid), C_mid[4]*(n_lo/n_mid),
-        C_mid[5]*(n_lo/n_mid), C_mid[6]*(n_lo/n_mid),
-        T_f_lo, T_c_lo, T_cold_lo, 0.0]
-x_hi = [C_mid[1]*(n_hi/n_mid), C_mid[2]*(n_hi/n_mid),
-        C_mid[3]*(n_hi/n_mid), C_mid[4]*(n_hi/n_mid),
-        C_mid[5]*(n_hi/n_mid), C_mid[6]*(n_hi/n_mid),
-        T_f_hi, T_c_hi, T_cold_hi, 0.0]
-X0 = Hyperrectangle(low=x_lo, high=x_hi)
-
-println("\n=== Tight-entry heatup PJ reach — saving full per-step envelopes ===")
-
-T_probe = 300.0
-sys  = BlackBoxContinuousSystem(rhs_heatup_pj_tf!, 10)
-prob = InitialValueProblem(sys, X0)
-alg  = TMJets(orderT=4, orderQ=2, abstol=1e-9, maxsteps=100000)
-t_start = time()
-sol = solve(prob; T=T_probe, alg=alg)
-println("  Wall time: $(round(time() - t_start; digits=1)) s")
-flow_hr = overapproximate(flowpipe(sol), Hyperrectangle)
-n_steps = length(flow_hr)
-println("  $n_steps reach-sets")
-
-t_arr     = zeros(n_steps)
-Tc_lo_ts  = zeros(n_steps); Tc_hi_ts  = zeros(n_steps)
-Tf_lo_ts  = zeros(n_steps); Tf_hi_ts  = zeros(n_steps)
-Tco_lo_ts = zeros(n_steps); Tco_hi_ts = zeros(n_steps)
-n_lo_ts   = zeros(n_steps); n_hi_ts   = zeros(n_steps)
-rho_lo_ts = zeros(n_steps); rho_hi_ts = zeros(n_steps)
-
-for (k, R) in enumerate(flow_hr)
-    s = set(R)
-    t_arr[k]    = high(s, 10)
-    Tc_lo_ts[k] = low(s, 8);  Tc_hi_ts[k] = high(s, 8)
-    Tf_lo_ts[k] = low(s, 7);  Tf_hi_ts[k] = high(s, 7)
-    Tco_lo_ts[k] = low(s, 9); Tco_hi_ts[k] = high(s, 9)
-    sumLC_lo = LAM_1*low(s,1) + LAM_2*low(s,2) + LAM_3*low(s,3) +
-               LAM_4*low(s,4) + LAM_5*low(s,5) + LAM_6*low(s,6)
-    sumLC_hi = LAM_1*high(s,1) + LAM_2*high(s,2) + LAM_3*high(s,3) +
-               LAM_4*high(s,4) + LAM_5*high(s,5) + LAM_6*high(s,6)
-    # Ramp-ref bounds at this reach-set's t.
-    t_hi_here = high(s, 10)
-    t_lo_here = low(s, 10)
-    # T_ref at these times (monotone increasing):
-    Tref_lo = T_STANDBY + RAMP_RATE_CS * t_lo_here
-    Tref_hi = T_STANDBY + RAMP_RATE_CS * t_hi_here
-    rho_lo_here = KP_HEATUP * (Tref_lo - high(s, 8))
-    rho_hi_here = KP_HEATUP * (Tref_hi - low(s, 8))
-    rho_lo_ts[k] = rho_lo_here
-    rho_hi_ts[k] = rho_hi_here
-    denom_lo = BETA - rho_hi_here
-    denom_hi = BETA - rho_lo_here
-    if denom_lo > 0
-        n_lo_ts[k] = LAMBDA * sumLC_lo / denom_hi
-        n_hi_ts[k] = LAMBDA * sumLC_hi / denom_lo
-    end
-end
-
-println("  T_c envelope:   [$(round(minimum(Tc_lo_ts); digits=2)), $(round(maximum(Tc_hi_ts); digits=2))] °C")
-println("  T_hot envelope: [$(round(minimum(2 .* Tc_lo_ts .- Tco_hi_ts); digits=2)), $(round(maximum(2 .* Tc_hi_ts .- Tco_lo_ts); digits=2))] °C")
-println("  T_cold env:     [$(round(minimum(Tco_lo_ts); digits=2)), $(round(maximum(Tco_hi_ts); digits=2))] °C")
-println("  rho envelope:   [$(round(minimum(rho_lo_ts); sigdigits=4)), $(round(maximum(rho_hi_ts); sigdigits=4))]")
-println("  rho / beta:     [$(round(minimum(rho_lo_ts)/BETA; digits=3)), $(round(maximum(rho_hi_ts)/BETA; digits=3))]")
-
-mat_out = joinpath(@__DIR__, "..", "..", "reachability", "reach_heatup_pj_tight_full.mat")
-matwrite(mat_out, Dict(
-    "t_arr"     => t_arr,
-    "Tc_lo_ts"  => Tc_lo_ts,  "Tc_hi_ts"  => Tc_hi_ts,
-    "Tf_lo_ts"  => Tf_lo_ts,  "Tf_hi_ts"  => Tf_hi_ts,
-    "Tco_lo_ts" => Tco_lo_ts, "Tco_hi_ts" => Tco_hi_ts,
-    "n_lo_ts"   => n_lo_ts,   "n_hi_ts"   => n_hi_ts,
-    "rho_lo_ts" => rho_lo_ts, "rho_hi_ts" => rho_hi_ts,
-    "T_probe"   => T_probe,
-    "beta"      => BETA,
-    "Kp"        => KP_HEATUP,
-    "T_c0"      => T_C0,
-    "T_cold0"   => T_COLD0,
-    "T_standby" => T_STANDBY,
-))
-println("\nSaved full per-step envelopes to $mat_out")
diff --git a/code/scripts/main_mode_sweep.jl b/code/scripts/sim/main_mode_sweep.jl
similarity index 89%
rename from code/scripts/main_mode_sweep.jl
rename to code/scripts/sim/main_mode_sweep.jl
index b2de003..9a7c228 100644
--- a/code/scripts/main_mode_sweep.jl
+++ b/code/scripts/sim/main_mode_sweep.jl
@@ -7,7 +7,7 @@
 # names (overwrites MATLAB outputs — the Julia versions take over).
 
 using Pkg
-Pkg.activate(joinpath(@__DIR__, ".."))
+Pkg.activate(joinpath(@__DIR__, "..", ".."))
 
 using Printf
 using LinearAlgebra
@@ -16,18 +16,18 @@ using Plots
 using MatrixEquations
 using JSON
 
-include(joinpath(@__DIR__, "..", "src", "pke_params.jl"))
-include(joinpath(@__DIR__, "..", "src", "pke_th_rhs.jl"))
-include(joinpath(@__DIR__, "..", "src", "pke_linearize.jl"))
-include(joinpath(@__DIR__, "..", "src", "pke_solver.jl"))
-include(joinpath(@__DIR__, "..", "src", "plot_pke_results.jl"))
-include(joinpath(@__DIR__, "..", "controllers", "controllers.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_params.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_th_rhs.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_linearize.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_solver.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "plot_pke_results.jl"))
+include(joinpath(@__DIR__, "..", "..", "controllers", "controllers.jl"))
 
 # --- Plant + predicates ---
 plant = pke_params()
 
 # Load T_standby from predicates.json for the hot-standby IC + heatup ref.
-pred_path = joinpath(@__DIR__, "..", "..", "reachability", "predicates.json")
+pred_path = joinpath(@__DIR__, "..", "..", "..", "reachability", "predicates.json")
 pred_raw  = JSON.parsefile(pred_path)
 T_standby = plant.T_c0 + pred_raw["derived"]["T_standby_offset_C"]
 
@@ -52,7 +52,7 @@ ctrl_op_lqr, K_lqr, _ = ctrl_operation_lqr_factory(
     plant, A, B; Q_lqr=Matrix(Q_lqr), R_lqr=R_lqr)
 
 # --- Figure output directory ---
-figdir = joinpath(@__DIR__, "..", "..", "docs", "figures")
+figdir = joinpath(@__DIR__, "..", "..", "..", "docs", "figures")
 isdir(figdir) || mkpath(figdir)
 
 # --- Run each mode ---
diff --git a/code/scripts/sim_heatup.jl b/code/scripts/sim/sim_heatup.jl
similarity index 98%
rename from code/scripts/sim_heatup.jl
rename to code/scripts/sim/sim_heatup.jl
index 4381b33..ac12d46 100644
--- a/code/scripts/sim_heatup.jl
+++ b/code/scripts/sim/sim_heatup.jl
@@ -9,7 +9,7 @@
 # trajectory that reach tubes can be checked against.
 
 using Pkg
-Pkg.activate(joinpath(@__DIR__, ".."))
+Pkg.activate(joinpath(@__DIR__, "..", ".."))
 
 using OrdinaryDiffEq
 
diff --git a/code/scripts/sim_sanity.jl b/code/scripts/sim/sim_sanity.jl
similarity index 83%
rename from code/scripts/sim_sanity.jl
rename to code/scripts/sim/sim_sanity.jl
index 2ce524f..244add1 100644
--- a/code/scripts/sim_sanity.jl
+++ b/code/scripts/sim/sim_sanity.jl
@@ -6,13 +6,13 @@
 # and prints the final state, which should agree with MATLAB to ~1e-3.
 
 using Pkg
-Pkg.activate(joinpath(@__DIR__, ".."))
+Pkg.activate(joinpath(@__DIR__, "..", ".."))
 
 using OrdinaryDiffEq
 
-include(joinpath(@__DIR__, "..", "src", "pke_params.jl"))
-include(joinpath(@__DIR__, "..", "src", "pke_th_rhs.jl"))
-include(joinpath(@__DIR__, "..", "controllers", "controllers.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_params.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_th_rhs.jl"))
+include(joinpath(@__DIR__, "..", "..", "controllers", "controllers.jl"))
 
 plant = pke_params()
 x0    = pke_initial_conditions(plant)
diff --git a/code/scripts/validate_pj.jl b/code/scripts/sim/validate_pj.jl
similarity index 91%
rename from code/scripts/validate_pj.jl
rename to code/scripts/sim/validate_pj.jl
index 894af04..719bbf9 100644
--- a/code/scripts/validate_pj.jl
+++ b/code/scripts/sim/validate_pj.jl
@@ -11,17 +11,17 @@
 # reconstructed n vs the dynamic n.  Quantify the error explicitly.
 
 using Pkg
-Pkg.activate(joinpath(@__DIR__, ".."))
+Pkg.activate(joinpath(@__DIR__, "..", ".."))
 
 using Printf
 using LinearAlgebra
 using OrdinaryDiffEq
 using Plots
 
-include(joinpath(@__DIR__, "..", "src", "pke_params.jl"))
-include(joinpath(@__DIR__, "..", "src", "pke_th_rhs.jl"))
-include(joinpath(@__DIR__, "..", "src", "pke_th_rhs_pj.jl"))
-include(joinpath(@__DIR__, "..", "controllers", "controllers.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_params.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_th_rhs.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_th_rhs_pj.jl"))
+include(joinpath(@__DIR__, "..", "..", "controllers", "controllers.jl"))
 
 plant = pke_params()
 T_standby = plant.T_c0 - 33.333333
@@ -88,7 +88,7 @@ for t_check in (1.0, 5.0, 10.0, 60.0, 300.0, 1200.0, 3000.0)
 end
 
 # --- Plot trajectory overlay ---
-figdir = joinpath(@__DIR__, "..", "..", "docs", "figures")
+figdir = joinpath(@__DIR__, "..", "..", "..", "docs", "figures")
 isdir(figdir) || mkpath(figdir)
 
 CtoF(T) = T*9/5 + 32
diff --git a/presentations/prelim-presentation/outline.md b/presentations/prelim-presentation/outline.md
new file mode 100644
index 0000000..001f140
--- /dev/null
+++ b/presentations/prelim-presentation/outline.md
@@ -0,0 +1,395 @@
+# Preliminary Example Presentation — Assertion-Evidence Outline
+
+**Format:** Assertion-evidence (Alley). Each slide: one declarative sentence
+at the top, one piece of visual evidence below. No bullet soup.
+
+**Duration:** 20 minutes. ~12 slides at ~1.5 min each (60 s typical, up to
+2 min for the heavier ones).
+
+**Audience:** OT-informed cybersecurity experts. Mostly CS, some control-theory
+familiarity, very little reactor-physics background. Can assume fluency with:
+LTL, automata, model-checking, reachability (as a concept), SMT/SAT. Must
+explain: PKE, reactivity, precursors, hybrid systems.
+
+**Running thread:** FRET natural-language requirement → LTL → synthesized
+discrete controller → continuous plant → per-mode reach-avoid proof
+→ hybrid correctness by composition.
+
+**Design principles:**
+- **Plots over bullets.** Every result slide anchors on one figure.
+- **Physical intuition before math.** Reactor basics before PKE.
+- **Honest limitations boxed on each result slide.** Audience are cyber
+  folks — they respect limits more than triumphs.
+- **CS vocabulary by default, engineering terms defined inline.**
+- **End with the seam**, not with a victory lap. The thesis question is
+  "how do discrete proofs and continuous proofs compose?" not "we verified
+  everything."
+
+---
+
+## Slide 1 — Title + hook
+
+**Assertion:** Verified hybrid control for nuclear startup is a real problem
+with no deployed solution.
+
+**Evidence:** Title block + a single image showing a control-room board of
+old-school analog gauges adjacent to a digital SCADA screen. Caption: "The
+gap we're filling."
+
+**Speak:** Introduce self + advisor + NRC fellowship context. Name the problem:
+nuclear reactor operations are 70-80% human-error-driven at severe-accident
+level (Chernobyl, TMI, Fukushima all primarily human factors). Plants run on
+paper procedures that have no formal verification. The most advanced work to
+date (HARDENS, NRC-funded) verified the discrete trip system in isolation and
+stopped there. This talk: a preliminary example of bridging discrete
+requirements all the way to continuous-state safety proofs, with the seam
+between them as the research contribution.
+
+**Reference:** thesis §2, ¶1-4. Cites Kemeny1979, Hogberg_2013, Kiniry2024.
+
+**Figures to make:** find or compose the control-room image.
+
+---
+
+## Slide 2 — The hybrid-systems gap
+
+**Assertion:** Existing tools verify either discrete *or* continuous systems,
+never the seam between them.
+
+**Evidence:** Two-column diagram. Left column: discrete (FRET, Spot,
+SAL, HARDENS, TLA+). Right column: continuous (MATLAB, Modelica, CORA,
+Flow*, HyTech). A dashed line between them labeled "the bridge problem."
+
+**Speak:** HARDENS got to TRL 2-3 on the discrete side alone. Continuous-side
+reach tools exist but assume the discrete controller is given as input. What
+nobody has done: produce an *end-to-end* proof where the discrete controller's
+transitions actually fire in physical time on the modeled plant.
+
+**Reference:** thesis §2.2 "Formal Methods", §2.2.1 HARDENS, §2.2.2 Temporal
+Logic. Thesis §2.3 (if exists — continuous methods).
+
+---
+
+## Slide 3 — Running example: the PWR_HYBRID_3 controller
+
+**Assertion:** Our target is a 4-mode hybrid controller that takes a PWR
+from hot standby to full power and back.
+
+**Evidence:** The DRC state diagram (from `fret-pipeline/diagrams/PWR_HYBRID_3_DRC_states.png`).
+Four nodes — shutdown, heatup, operation, scram — with labeled transitions.
+
+**Speak:** One sentence of nuclear physics: a PWR has neutron population
+(power), coolant temperatures at three locations, and precursor concentrations
+that determine delayed-neutron generation. Startup is shutdown → heatup →
+operation; any fault triggers scram. Mode transitions gated by boolean
+conditions like `t_avg_in_range ∧ p_above_crit ∧ inv1_holds`.
+
+Four modes means four *physical* behaviors the plant has to exhibit for
+the discrete logic to be sound. This example stresses every layer of the
+bridge we're building.
+
+**Evidence path:** `fret-pipeline/diagrams/PWR_HYBRID_3_DRC_states.png`,
+`fret-pipeline/pwr_hybrid_3.json` for one example requirement.
+
+---
+
+## Slide 4 — Layer 1: FRET → LTL → synthesized DRC
+
+**Assertion:** Natural-language requirements become a verified discrete
+controller automatically.
+
+**Evidence:** One-panel figure: left half = snippet of FRET natural language
+("Upon `control_mode = q_heatup ∧ t_avg_in_range ∧ p_above_crit ∧ inv1_holds`
+DRC shall at the next timepoint satisfy `control_mode = q_operation`"), right
+half = the corresponding LTL, bottom = arrow labeled `ltlsynt` pointing to
+the synthesized state-diagram node.
+
+**Speak:** FRET compiles natural-language requirements into LTL (Spot
+ecosystem). `ltlsynt` turns LTL into an AIGER circuit representing the
+minimal correct discrete controller. This is well-trodden ground in CS-land;
+our contribution is *using it* on reactor operating procedures — so far
+a formal-methods-free domain.
+
+**Evidence path:** `fret-pipeline/README.md`, `pwr_hybrid_3.json` sample
+requirement, `fret-pipeline/circuits/PWR_HYBRID_3_DRC.aag` AIGER circuit.
+
+---
+
+## Slide 5 — Layer 2: the plant model
+
+**Assertion:** A 10-state PKE + thermal-hydraulics model is faithful enough
+for reach analysis and tractable enough to actually compute with.
+
+**Evidence:** The state-vector diagram: `x = [n, C_1, ..., C_6, T_f, T_c, T_cold]`
+with arrows showing the coupling — neutron balance → fuel heating → coolant
+transport → feedback to reactivity.
+
+**Speak:** Point-kinetic equations + three thermal nodes. 10 states, 1
+control input (rod worth u), 1 disturbance (steam-generator heat removal Q_sg).
+Stiff system: prompt-neutron timescale is 10⁻⁴ s, thermal timescales are
+seconds to minutes. That stiffness becomes critical later — it's what forces
+the singular-perturbation move.
+
+Five controllers: one per DRC mode (shutdown, heatup, operation under P, under
+LQR, scram). All Julia, all in `code/src/` + `code/controllers/`.
+
+**Evidence path:** `code/src/pke_th_rhs.jl`, `journal/journal.pdf` entry
+2026-04-17 for derivation.
+
+**Figures to make:** state-vector coupling diagram (can be ASCII → inkscape).
+
+---
+
+## Slide 6 — Layer 3: the reach-avoid obligation per mode
+
+**Assertion:** Discrete correctness transfers to continuous correctness
+via one reach-avoid obligation per mode.
+
+**Evidence:** Taxonomy table from `reachability/WALKTHROUGH.md`:
+
+| Mode | Kind | Obligation |
+|---|---|---|
+| shutdown | equilibrium | stay in X_safe forever |
+| heatup | transition | from X_entry, reach X_exit in [T_min, T_max] while maintaining inv1 |
+| operation | equilibrium | stay in X_safe forever under bounded Q_sg |
+| scram | transition | from any mode, drive n safely subcritical in T_max |
+
+**Speak:** This is the structural insight. Equilibrium modes → forever-invariance
+proofs. Transition modes → reach-avoid proofs with time budgets. The DRC
+transitions fire iff the reach set enters the right exit predicate. If
+every mode's obligation is discharged, the hybrid system is correct by
+composition. The methodological contribution of the chapter.
+
+**Reference:** `reachability/WALKTHROUGH.md` §1, `reachability/predicates.json`
+→ `mode_boundaries`.
+
+---
+
+## Slide 7 — Operation mode: discharge all 6 safety halfspaces (the win)
+
+**Assertion:** Under LQR feedback and ±15% Q_sg variation, the operation-mode
+reach tube discharges every safety halfspace with orders of magnitude of margin.
+
+**Evidence:** The two-panel figure `docs/figures/reach_operation_tubes.png`.
+Left: temperature tubes (T_c, T_hot, T_cold) overlaid; the near-zero width
+of T_c visually demonstrates LQR contraction. Right: ΔT_core tube with a
+right-axis in MW showing the power is tight around nominal.
+
+Per-halfspace margin table (inset or second slide if space allows):
+
+| Halfspace | Limit | Reach max | Margin |
+|---|---:|---:|---:|
+| fuel_centerline | 1200 °C | 328.9 | **+871 K** |
+| t_avg_high_trip | 320 °C | 308.4 | **+11.6 K** |
+| n_high_trip | 1.15 | 1.012 | **+0.14** |
+
+**Speak:** LQR contracts the regulated direction (T_c) from 0.1 K to 0.03 K
+halfwidth. Uncontrolled directions (precursors) expand but don't appear in
+any safety predicate. Tightest margin is the high-flux trip at 12% headroom.
+
+**Limitations box:** linear-reach of a nonlinear plant — *approximate, not
+sound* yet for the physical plant. That's what the next slides address.
+
+**Evidence path:** `code/scripts/reach/reach_operation.jl`,
+`docs/figures/reach_operation_tubes.png`.
+
+---
+
+## Slide 8 — Quadratic Lyapunov barrier fails structurally
+
+**Assertion:** The canonical quadratic Lyapunov barrier cannot certify this
+plant — *even with LQR inside the barrier* — because of plant anisotropy.
+
+**Evidence:** Bar chart comparing open-loop vs LQR-closed-loop Lyapunov bound
+on `n_high_trip` direction: OL is 27 million × nominal, CL is 1,242 × nominal.
+Title: "LQR helps 20,000× but both bounds are physically meaningless."
+
+**Speak:** Plant has prompt timescale 10⁻⁴ s vs thermal timescale ~10 s.
+Lyapunov matrix P inherits that 10⁴ spread. Resulting ellipsoid stretches
+obscenely along the fast directions when projected to physical coordinates.
+No amount of controller tuning fixes this — it's set by the plant's own
+physics. The ellipsoid geometry cannot match the slab-shaped safe region.
+**This is the motivation for polynomial (SOS) or polytopic barriers.**
+
+**Reference:** `code/scripts/barrier/barrier_compare_OL_CL.jl`, journal entry
+2026-04-20 (evening).
+
+---
+
+## Slide 9 — Prompt-jump reduction makes nonlinear reach tractable
+
+**Assertion:** Singular-perturbation reduction of the fast neutronics gives
+us 30× more reach horizon and a rigorous O(Λ) error bound via Tikhonov.
+
+**Evidence:** Side-by-side: left panel = validate_pj_heatup.png (empirical:
+PJ and full-state trajectories overlay perfectly over 50 minutes). Right panel
+= the Tikhonov bound written out mathematically:
+$$|x(t) - x_{\mathrm{PJ}}(t)| \leq C \cdot \Lambda = \mathcal{O}(10^{-4})$$
+
+**Speak:** Setting dn/dt = 0 and solving algebraically for n gives us
+`n = Λ·Σλᵢ·Cᵢ / (β - ρ)`. State drops from 10 to 9, removes the Λ⁻¹
+stiffness. Reach tool (TMJets) can take steps on thermal timescale instead
+of prompt timescale.
+
+The magic: we *prove* the PJ approximation's validity condition (`β - ρ > 0`
+with margin) **as part of the same reach obligation** via the
+`prompt_critical_margin_heatup` halfspace in `inv1_holds`. No hand-waving —
+the reach machinery proves both the safety claim and the approximation's
+soundness together.
+
+**Evidence path:** `code/src/pke_th_rhs_pj.jl`, `code/scripts/sim/validate_pj.jl`,
+`journal/` entry 2026-04-21 "Tikhonov bound".
+
+---
+
+## Slide 10 — First sound nonlinear heatup reach (the second win)
+
+**Assertion:** With PJ + a tightened entry box, nonlinear Taylor-model reach
+discharges all 6 `inv1_holds` safety halfspaces over the first 300 s of heatup.
+
+**Evidence:** The four-panel `reach_heatup_pj_tubes.png`: temperature tubes
+overlaid, ΔT_core, ρ in dollars (stays between -0.25 $ and -0.05 $ —
+sub-prompt-critical), n decaying.
+
+**Speak:** TMJets Taylor-model integration. 12,932 reach-sets, 200 s wall time.
+Tube **stable** — T_c envelope `[281.05, 291.0]` identical at 60 s and 300 s,
+meaning the controller holds the state inside the tube indefinitely. This
+is the first sound (modulo O(Λ) PJ error) nonlinear reach-avoid artifact
+for this plant.
+
+**Limitations box:** 300 s of a 5-hour `T_max`. Step budget caps horizon; entry
+refinement (Blanchini-style) is the path to hours.
+
+**Evidence path:** `code/scripts/reach/reach_heatup_pj.jl configs/heatup/tight.toml`,
+`docs/figures/reach_heatup_pj_tubes.png`.
+
+---
+
+## Slide 11 — Degree-4 SOS polynomial barrier, a working proof of concept
+
+**Assertion:** A degree-4 polynomial barrier certificate — the structure that
+quadratic Lyapunov cannot provide — works on the operation-mode projection.
+
+**Evidence:** Equation block showing the symbolic polynomial
+(abbreviated — the actual 13-term polynomial from `code/scripts/barrier/barrier_sos_2d.jl`).
+Side caption: "CSDP returned OPTIMAL. Solve: ~3 seconds."
+
+**Speak:** Polynomials of degree 4 can localize to slab-shaped safety regions
+in a way degree-2 (quadratic) cannot. The Prajna-Jadbabaie formulation: find
+`B(x)` such that `B ≤ 0` on entry, `B ≥ 0` on unsafe, and `∇B·f ≤ 0` on
+`{B=0}`. Sum-of-squares programming reduces this to an SDP. 2D projection
+proof of concept; scaling to full 10-state requires bigger solver (Mosek or SCS).
+
+**Reference:** `code/scripts/barrier/barrier_sos_2d.jl`, journal entry
+2026-04-21.
+
+---
+
+## Slide 12 — The hybrid correctness story, assembled
+
+**Assertion:** All four pieces — FRET/synthesis, plant model, reach tubes,
+SOS/polytopic barriers — compose into a hybrid system correctness argument.
+
+**Evidence:** The composition picture. DRC state diagram at center; each
+DRC node has an arrow to a "per-mode obligation box" labeled with its proof
+artifact (tube or barrier or certificate). Arrows between nodes = transitions,
+labeled with the predicate polytope. Outer dashed box: "hybrid system
+closed-loop safety."
+
+**Speak:** What's been built, what's been proven, what's next. Transition
+correctness between modes is the next thesis-blocking piece — each mode's
+`X_exit` has to be inside the next mode's `X_entry` (with margin). That's an
+inclusion check, not a reach — but it's the final structural piece.
+
+---
+
+## Slide 13 — Limitations and next steps
+
+**Assertion:** This is a preliminary example, not a deployable system.
+
+**Evidence:** A two-column list:
+
+| What's proven | What's next |
+|---|---|
+| Operation mode safe under ±15% load (approximate, via linear reach) | Nonlinear operation reach (tight SOS barrier, LOCA disturbance) |
+| Heatup PJ reach discharges all 6 safety halfspaces for 300 s | Extend to full 5-hr `T_max`; model steam-dump disturbance |
+| Scram PJ n-decay monotone over 60 s | Expand `X_entry(scram)` to union of all mode tubes + LOCA |
+| Degree-4 SOS barrier on 2D projection | Full 10-state SOS barrier |
+| PJ error rigorously bounded by Tikhonov | Parametric α uncertainty; DNBR correlation |
+
+**Speak:** The gap from prelim example to deployable system is well-defined:
+more states, tighter bounds, real tech-spec numbers, hardware-in-the-loop.
+None of these is a research novelty unto itself — they're engineering.
+The research contribution is the composition framework, demonstrated end-to-end
+on a nontrivial running example. Phase 2 of the thesis is filling in the
+gaps and expanding to multiple plants.
+
+---
+
+## Slide 14 — Questions / acknowledgements
+
+**Assertion:** (none — backup slide)
+
+**Evidence:** Thanks to advisor, committee, NRC fellowship. Links to:
+repo (gitea), journal PDF, thesis in progress.
+
+**Speak:** Open for questions. Expect questions on:
+- "Why not Stateflow/Simulink?" → (have answer prepared; HARDENS used
+  Cryptol for a reason — formal tool integration matters)
+- "What's the relationship to LTLMoP / DragonFly?" → (survey answer)
+- "How would this interact with ML components?" → (out of scope for now,
+  the whole pitch is *no ML in safety-critical loop*)
+- "What's your threat model for cybersecurity?" → (tie back to OT audience —
+  formal methods guarantee the controller's logic; they do not protect the
+  comms layer or implementation-level vulns. Mention HARDENS' focus on
+  "correctness of implementation" vs our focus on "correctness of
+  specification")
+
+---
+
+## Presentation construction notes
+
+### What to build in Beamer later
+- Assertion-evidence template (one sentence at top, centered figure below).
+- Consistent color coding: green for "discharged" / "proved", red for
+  "limitation" / "gap", blue for "discrete-layer", amber for "continuous-layer".
+- The composition diagram on slide 12 — this will take the most Tikz work.
+
+### Figures that need to be created or dressed up for the talk
+- Slide 1: control-room visual.
+- Slide 2: discrete/continuous tools comparison.
+- Slide 3: use `fret-pipeline/diagrams/PWR_HYBRID_3_DRC_states.png`.
+- Slide 4: FRET → LTL → AIGER panel.
+- Slide 5: 10-state coupling diagram.
+- Slide 7, 10: reuse `docs/figures/reach_*_tubes.png`.
+- Slide 8: bar chart (new, easy matplotlib).
+- Slide 9: reuse `validate_pj_heatup.png`.
+- Slide 11: latex-typeset polynomial + CSDP log snippet.
+- Slide 12: composition diagram (Tikz, will take time).
+
+### Cybersecurity angle to emphasize
+- The point the OT audience will care about: this kind of verification
+  **constrains what the controller can physically do**. Even if an attacker
+  gets past authentication and can inject arbitrary commands, the DRC-plus-
+  reach-certified envelope limits how bad things can get *within the
+  physical plant*. That's a different assurance axis than the usual
+  comms-security one and complements it.
+- Formal methods as defense-in-depth: they catch bugs *before* deployment,
+  which reduces attack surface more than any runtime defense.
+- The PJ reduction + Tikhonov approach might be of interest for other
+  safety-critical stiff systems (power grid, aerospace).
+
+### Things to NOT do
+- Don't get lost in reactor-physics detail. One sentence per physical
+  concept; get them to the CS content fast.
+- Don't show code unless it's a slide about *why* the code structure
+  matters. Code screenshots are terrible evidence.
+- Don't oversell. Honest limitations at every stage builds trust with
+  a skeptical audience.
+- Don't use more than 2 bullet points on any slide. Alley rules.
+
+### Timing checkpoints
+- Slide 6 (mode-obligation taxonomy) by minute 8.
+- Slide 10 (first nonlinear reach result) by minute 13.
+- Slide 12 (composition story) by minute 17.
+- Slide 13 (limitations) by minute 19.
diff --git a/results/.gitignore b/results/.gitignore
new file mode 100644
index 0000000..e63f75c
--- /dev/null
+++ b/results/.gitignore
@@ -0,0 +1 @@
+*.mat
diff --git a/results/README.md b/results/README.md
new file mode 100644
index 0000000..c31e786
--- /dev/null
+++ b/results/README.md
@@ -0,0 +1,32 @@
+# results/
+
+Output artifacts from reach / barrier / validation runs.  All files
+in this directory are **gitignored** — they're regenerated by running
+the corresponding script.
+
+## What lives here
+
+Named `<analysis>_<variant>.mat` (MATLAB v7 format, read/written via
+Julia's MAT.jl for historical convenience — switching to JLD2 is a
+deferred item, see `journal/` for context).
+
+Generated by (as of 2026-04-21):
+
+| File | Generated by |
+|---|---|
+| `reach_operation_result.mat`      | `code/scripts/reach/reach_operation.jl` |
+| `reach_heatup_pj_tight_full.mat`  | `code/scripts/reach/reach_heatup_pj_tight_full.jl` |
+| `reach_scram_pj_result.mat`       | `code/scripts/reach/reach_scram_pj.jl` |
+
+Consumed by:
+
+- `code/scripts/plot/plot_reach_tubes.jl` — tube overlay figures.
+- `app/predicate_explorer.jl` — live per-halfspace margin rendering.
+
+## Why this directory exists
+
+`reachability/` used to hold both source-of-truth specs
+(`predicates.json`, README, WALKTHROUGH.md) and ephemeral reach results.
+Mixing stable vs.\ regenerated content made it hard to tell which was
+which at a glance.  Split: `reachability/` = specs + docs,
+`results/` = outputs.