diff --git a/julia-port/.gitignore b/julia-port/.gitignore
new file mode 100644
index 0000000..ba39cc5
--- /dev/null
+++ b/julia-port/.gitignore
@@ -0,0 +1 @@
+Manifest.toml
diff --git a/julia-port/Project.toml b/julia-port/Project.toml
new file mode 100644
index 0000000..1f95125
--- /dev/null
+++ b/julia-port/Project.toml
@@ -0,0 +1,14 @@
+authors = ["Dane Sabo <yourstruly@danesabo.com>"]
+
+[deps]
+LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
+MatrixEquations = "99c1a7ee-ab34-5fd5-8076-27c950a045f4"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+ReachabilityAnalysis = "1e97bd63-91d1-579d-8e8d-501d2b57c93f"
+
+[compat]
+OrdinaryDiffEq = "6.111.0"
+julia = "1.10"
diff --git a/julia-port/README.md b/julia-port/README.md
new file mode 100644
index 0000000..55e46a4
--- /dev/null
+++ b/julia-port/README.md
@@ -0,0 +1,40 @@
+# Julia port
+
+Parallel implementation of the plant model (`../plant-model/`) in Julia,
+intended as an eventual landing spot for reachability if we outgrow the
+MATLAB / hand-rolled tooling.
+
+## Status
+
+- `src/pke_params.jl`, `src/pke_th_rhs.jl`, `src/pke_linearize.jl` —
+  functional, match MATLAB term-for-term.
+- `controllers/controllers.jl` — all five modes ported
+  (null, shutdown, heatup, operation, scram).  LQR factory included via
+  MatrixEquations.jl.
+- `scripts/sim_sanity.jl` — closed-loop simulation matches MATLAB to 3
+  decimals on the 100% → 80% `Q_sg` step.  ✅ passing.
+- `scripts/reach_operation.jl` — ❌ stub.  ReachabilityAnalysis.jl
+  algorithms blow up on this stiff, badly-conditioned system.  See the
+  file header for diagnosis and planned fix (state rescaling).
+
+## When to prefer Julia over MATLAB
+
+Today: nowhere.  The sanity path exists so we don't regret the eventual
+port.
+
+Once the reach scaling is resolved:
+- Nonlinear reach for the P controller in operation mode (CORA /
+  JuliaReach territory where MATLAB's linearization doesn't suffice).
+- Heatup reach with its time-varying reference.
+- Parametric studies where MATLAB license fees / CI friction matter.
+
+## Running
+
+```bash
+cd julia-port
+julia --project=. -e 'using Pkg; Pkg.instantiate()'
+julia --project=. scripts/sim_sanity.jl
+```
+
+First run will precompile the dependency stack (OrdinaryDiffEq,
+ReachabilityAnalysis, LazySets — a few minutes).
diff --git a/julia-port/controllers/controllers.jl b/julia-port/controllers/controllers.jl
new file mode 100644
index 0000000..62707bf
--- /dev/null
+++ b/julia-port/controllers/controllers.jl
@@ -0,0 +1,49 @@
+"""
+Mode controllers — signatures match the MATLAB side:
+    u = ctrl_<mode>(t, x, plant, ref)
+
+Pure functions.  `ref` can be any NamedTuple of setpoints; unused fields
+are ignored.  For heatup, `ref` must provide `T_start`, `T_target`, `ramp_rate`.
+"""
+
+ctrl_null(t, x, plant, ref) = 0.0
+
+ctrl_shutdown(t, x, plant, ref) = -5.0 * plant.beta
+ctrl_scram(t, x, plant, ref)    = -8.0 * plant.beta
+
+function ctrl_operation(t, x, plant, ref)
+    Kp = 1e-4
+    T_avg = x[9]
+    return Kp * (ref.T_avg - T_avg)
+end
+
+function ctrl_heatup(t, x, plant, ref)
+    Kp = 1e-4
+    T_f   = x[8]
+    T_avg = x[9]
+    u_ff = -plant.alpha_f * (T_f   - plant.T_f0) -
+            plant.alpha_c * (T_avg - plant.T_c0)
+    T_ref = min(ref.T_start + ref.ramp_rate * t, ref.T_target)
+    u_unsat = u_ff + Kp * (T_ref - T_avg)
+    u_min = get(ref, :u_min, -5 * plant.beta)
+    u_max = get(ref, :u_max,  0.5 * plant.beta)
+    return clamp(u_unsat, u_min, u_max)
+end
+
+"""
+    ctrl_operation_lqr_factory(plant; Q_lqr, R_lqr)
+
+Returns a closure `ctrl(t, x, plant_ignored, ref_ignored)` with the LQR
+gain baked in.  Pattern chosen so the user can specify Q/R from the
+calling script and get a pure function to pass to the ODE solver.
+
+Depends on MatrixEquations.jl for `arec` (algebraic Riccati).
+"""
+function ctrl_operation_lqr_factory(plant, A, B; Q_lqr, R_lqr)
+    x_op = pke_initial_conditions(plant)
+    X_ric, _, _ = MatrixEquations.arec(A, B, R_lqr, Q_lqr)
+    K = (R_lqr \ B') * X_ric   # 1x10
+    return function (t, x, plant_ignored, ref_ignored)
+        return -(K * (x - x_op))[1]
+    end, K, x_op
+end
diff --git a/julia-port/scripts/reach_operation.jl b/julia-port/scripts/reach_operation.jl
new file mode 100644
index 0000000..f4e8770
--- /dev/null
+++ b/julia-port/scripts/reach_operation.jl
@@ -0,0 +1,90 @@
+#!/usr/bin/env julia
+#
+# reach_operation.jl — operation-mode linear reach in Julia.  **STUB**.
+#
+# Attempt to reproduce reachability/reach_operation.m using
+# ReachabilityAnalysis.jl.
+#
+# STATUS: does NOT yet produce a useful result.  The closed-loop system
+# is strongly stiff (eigvals spanning 0.012 to 65 /s) and has states
+# with magnitudes differing by ~10 orders of magnitude (precursors C_i
+# ~ 1e5 due to 1/Lambda, temperatures ~ 300).  LGG09 / GLGM06 / BFFPSV18
+# all blow up numerically over the 600s horizon, returning envelopes
+# ranging from 2757 K to 1e37 K instead of the known-tight 0.03 K.
+#
+# Likely fix: diagonal state scaling S such that A_cl_scaled = S A_cl S^{-1}
+# has O(1) entries, run reach in scaled coordinates, then invert S.
+# Also worth trying: ASB07 (adaptive step) or the Taylor-model schemes.
+# Left as future work — the hand-rolled MATLAB zonotope in
+# reachability/reach_linear.m gives a valid result today, so the Julia
+# port is priority-B.
+
+using Pkg
+Pkg.activate(joinpath(@__DIR__, ".."))
+
+using LinearAlgebra
+using MatrixEquations
+using ReachabilityAnalysis, LazySets
+
+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)
+
+# --- Linearize around the operating point ---
+A, B, B_w, _, _, _ = pke_linearize(plant)
+
+# --- LQR gain (same Q, R as the MATLAB side) ---
+Q_lqr = Diagonal([1.0, 1e-3, 1e-3, 1e-3, 1e-3, 1e-3, 1e-3, 1e-2, 1e2, 1.0])
+R_lqr = 1e6 * ones(1, 1)
+X_ric, _, _ = arec(A, reshape(B, :, 1), R_lqr, Matrix(Q_lqr))
+K   = (R_lqr \ reshape(B, 1, :) * X_ric)    # 1 x 10
+A_cl = A - reshape(B, :, 1) * K
+
+println("\n=== Julia closed-loop spectrum ===")
+ev = eigvals(A_cl)
+println("  max Re(eig) = ", round(maximum(real.(ev)); sigdigits=4))
+println("  min Re(eig) = ", round(minimum(real.(ev)); sigdigits=4))
+
+# --- Reach-set problem: dx/dt = A_cl dx + B_w w, dx(0) ∈ X0, w ∈ W ---
+delta_entry = [0.01 * x_op[1];
+               0.001 .* abs.(x_op[2:7]);
+               0.1; 0.1; 0.1]
+X0 = Hyperrectangle(zeros(10), delta_entry)
+
+Q_nom = plant.P0
+w_lo, w_hi = -0.15 * plant.P0, 0.0
+W  = Interval(w_lo, w_hi)
+B_w_col = reshape(B_w, :, 1)
+
+# Encode the bounded disturbance as a bounded input:
+#   dx/dt = A_cl x + B_w u,   u ∈ W.
+# Direct constructor since the @system macro's dialect is package-version sensitive.
+sys = ConstrainedLinearControlContinuousSystem(A_cl, B_w_col, Universe(10), W)
+prob = InitialValueProblem(sys, X0)
+sol  = solve(prob; T=600.0, alg=GLGM06(; δ=0.5, max_order=20))
+
+# Extract T_c envelope via support function queries.
+flow   = flowpipe(sol)
+e9     = zeros(10); e9[9] = 1.0
+T_c_hi = [ρ(e9,  set(R)) for R in flow]
+T_c_lo = [-ρ(-e9, set(R)) for R in flow]
+
+println("\n=== Operation-mode reach envelope on T_c (deviation from T_c0) ===")
+println("  min dT_c = ", round(minimum(T_c_lo); digits=4), " K")
+println("  max dT_c = ", round(maximum(T_c_hi); digits=4), " K")
+println("  safety band |dT_c| <= 5.0 K")
+if maximum(abs.([T_c_lo; T_c_hi])) <= 5.0
+    println("  OK: Julia reach set inside safety band.")
+else
+    println("  *** VIOLATION ***")
+end
+
+# Save the envelope for later comparison
+using MAT
+matwrite(joinpath(@__DIR__, "..", "..", "reachability", "julia_reach_operation.mat"),
+         Dict("T_c_lo" => T_c_lo, "T_c_hi" => T_c_hi,
+              "A_cl" => A_cl, "K" => Matrix(K), "delta_entry" => delta_entry))
+println("\nSaved envelope to reachability/julia_reach_operation.mat")
diff --git a/julia-port/scripts/sim_sanity.jl b/julia-port/scripts/sim_sanity.jl
new file mode 100644
index 0000000..2ce524f
--- /dev/null
+++ b/julia-port/scripts/sim_sanity.jl
@@ -0,0 +1,40 @@
+#!/usr/bin/env julia
+#
+# sim_sanity.jl — verify the Julia port matches the MATLAB result.
+#
+# Reproduces main.m Run 2 (ctrl_operation under 100% -> 80% Q_sg step)
+# and prints the final state, which should agree with MATLAB to ~1e-3.
+
+using Pkg
+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"))
+
+plant = pke_params()
+x0    = pke_initial_conditions(plant)
+
+Q_sg = t -> plant.P0 * (1.0 - 0.2 * (t >= 30))
+ref  = (; T_avg = plant.T_c0)
+
+function rhs!(dx, x, p, t)
+    u = ctrl_operation(t, x, plant, ref)
+    pke_th_rhs!(dx, x, t, plant, Q_sg, u)
+end
+
+prob = ODEProblem(rhs!, x0, (0.0, 600.0))
+sol  = solve(prob, Rodas5(); reltol=1e-8, abstol=1e-10)
+
+xf = sol.u[end]
+CtoF(T) = T * 9/5 + 32
+println("\n=== Julia port sanity — ctrl_operation under 100% -> 80% Q_sg step ===")
+println("  Final t = ", sol.t[end])
+println("  n       = $(round(xf[1]; digits=4))   (expect ~0.800)")
+println("  T_f     = $(round(CtoF(xf[8]); digits=2)) F   (expect ~616.6)")
+println("  T_avg   = $(round(CtoF(xf[9]); digits=2)) F   (expect ~587.8)")
+println("  T_cold  = $(round(CtoF(xf[10]); digits=2)) F   (expect ~561.4)")
+u_final = ctrl_operation(sol.t[end], xf, plant, ref)
+println("  u       = $(round(u_final/plant.beta; digits=4)) \$   (expect ~-0.0068)")
diff --git a/julia-port/src/pke_linearize.jl b/julia-port/src/pke_linearize.jl
new file mode 100644
index 0000000..01fb716
--- /dev/null
+++ b/julia-port/src/pke_linearize.jl
@@ -0,0 +1,37 @@
+"""
+    pke_linearize(plant; x_star=nothing, u_star=0.0, Q_star=nothing)
+
+Numerical Jacobians via central finite differences.  Returns
+`(A, B, B_w, x_star, u_star, Q_star)` such that for small (dx, du, dw):
+
+    dx/dt ≈ A dx + B du + B_w dw,   where w = Q_sg.
+
+Defaults: `x_star` = operating-point steady state, `u_star = 0`,
+`Q_star = P0`.
+"""
+function pke_linearize(plant; x_star=nothing, u_star=0.0, Q_star=nothing)
+    x_star === nothing && (x_star = pke_initial_conditions(plant))
+    Q_star === nothing && (Q_star = plant.P0)
+
+    n = length(x_star)
+    eps_rel = 1e-6
+    eps_abs = 1e-8
+
+    f = (x, u, w) -> pke_th_rhs(x, 0.0, plant, t -> w, u)
+
+    A = zeros(n, n)
+    for k in 1:n
+        h = max(eps_rel * abs(x_star[k]), eps_abs)
+        xp = copy(x_star);  xp[k] += h
+        xm = copy(x_star);  xm[k] -= h
+        A[:, k] = (f(xp, u_star, Q_star) - f(xm, u_star, Q_star)) ./ (2h)
+    end
+
+    h = max(eps_rel * abs(u_star), eps_abs)
+    B = (f(x_star, u_star + h, Q_star) - f(x_star, u_star - h, Q_star)) ./ (2h)
+
+    h = max(eps_rel * abs(Q_star), 1.0)
+    B_w = (f(x_star, u_star, Q_star + h) - f(x_star, u_star, Q_star - h)) ./ (2h)
+
+    return A, B, B_w, x_star, u_star, Q_star
+end
diff --git a/julia-port/src/pke_params.jl b/julia-port/src/pke_params.jl
new file mode 100644
index 0000000..a1db046
--- /dev/null
+++ b/julia-port/src/pke_params.jl
@@ -0,0 +1,53 @@
+"""
+    pke_params()
+
+Return a NamedTuple with all plant parameters and derived steady-state
+conditions for the PKE + thermal-hydraulics model.  See
+../plant-model/pke_params.m for physical rationale on each value.
+
+Kept as a NamedTuple rather than a mutable struct so that reachability
+tools (`@taylorize` in ReachabilityAnalysis.jl) can inline the constants.
+"""
+function pke_params()
+    # --- Neutronics ---
+    Lambda = 1e-4
+    beta_i = [0.000215, 0.001424, 0.001274, 0.002568, 0.000748, 0.000273]
+    lambda_i = [0.0124, 0.0305, 0.111, 0.301, 1.14, 3.01]
+    beta = sum(beta_i)
+
+    # --- Thermal-hydraulic ---
+    P0   = 1000e6
+    M_f  = 50000.0
+    c_f  = 300.0
+    M_c  = 20000.0
+    c_c  = 5450.0
+    hA   = 5e7
+    W    = 5000.0
+    M_sg = 30000.0
+
+    # --- Reactivity feedback coefficients ---
+    alpha_f = -2.5e-5
+    alpha_c = -1.0e-4
+
+    # --- Derived steady-state (full-power equilibrium) ---
+    T_cold0 = 290.0
+    dT_core = P0 / (W * c_c)
+    T_hot0  = T_cold0 + dT_core
+    T_c0    = (T_hot0 + T_cold0) / 2
+    T_f0    = T_c0 + P0 / hA
+
+    return (; Lambda, beta_i, lambda_i, beta, P0, M_f, c_f, M_c, c_c, hA,
+             W, M_sg, alpha_f, alpha_c, T_cold0, T_hot0, T_c0, T_f0, dT_core)
+end
+
+"""
+    pke_initial_conditions(plant)
+
+Full-power steady state: n=1, precursors at equilibrium, temperatures at
+(T_f0, T_c0, T_cold0).  Returns a 10-element Vector.
+"""
+function pke_initial_conditions(plant)
+    n0 = 1.0
+    C0 = (plant.beta_i ./ (plant.lambda_i .* plant.Lambda)) .* n0
+    return [n0; C0; plant.T_f0; plant.T_c0; plant.T_cold0]
+end
diff --git a/julia-port/src/pke_th_rhs.jl b/julia-port/src/pke_th_rhs.jl
new file mode 100644
index 0000000..034c7ea
--- /dev/null
+++ b/julia-port/src/pke_th_rhs.jl
@@ -0,0 +1,59 @@
+"""
+    pke_th_rhs!(dx, x, t, plant, Q_sg, u)
+
+In-place ODE RHS for the coupled PKE + thermal-hydraulics model.
+Matches ../plant-model/pke_th_rhs.m term-for-term.
+
+State x = [n, C1..C6, T_f, T_c, T_cold].
+
+Arguments:
+- `dx`  : 10-element mutable output
+- `x`   : 10-element state
+- `t`   : time [s]
+- `plant`: parameter NamedTuple from `pke_params`
+- `Q_sg`: callable `Q_sg(t)` returning SG heat removal [W]
+- `u`   : scalar external reactivity [dk/k]
+"""
+function pke_th_rhs!(dx, x, t, plant, Q_sg, u)
+    n      = x[1]
+    T_f    = x[8]
+    T_c    = x[9]
+    T_cold = x[10]
+    T_hot  = 2 * T_c - T_cold
+
+    Qsg_t = Q_sg(t)
+
+    rho = u +
+          plant.alpha_f * (T_f - plant.T_f0) +
+          plant.alpha_c * (T_c - plant.T_c0)
+
+    # Neutronics
+    dx[1] = (rho - plant.beta) / plant.Lambda * n +
+            plant.lambda_i[1]*x[2] + plant.lambda_i[2]*x[3] +
+            plant.lambda_i[3]*x[4] + plant.lambda_i[4]*x[5] +
+            plant.lambda_i[5]*x[6] + plant.lambda_i[6]*x[7]
+
+    # Precursors
+    for i in 1:6
+        dx[1+i] = (plant.beta_i[i] / plant.Lambda) * n - plant.lambda_i[i] * x[1+i]
+    end
+
+    # Thermal-hydraulics
+    dx[8]  = (plant.P0 * n - plant.hA * (T_f - T_c)) / (plant.M_f * plant.c_f)
+    dx[9]  = (plant.hA * (T_f - T_c) - 2 * plant.W * plant.c_c * (T_c - T_cold)) /
+             (plant.M_c * plant.c_c)
+    dx[10] = (plant.W * plant.c_c * (T_hot - T_cold) - Qsg_t) /
+             (plant.M_sg * plant.c_c)
+
+    return nothing
+end
+
+"""
+Convenience: allocating version (returns dx as a new vector).  Useful in
+scripts and tests; for hot-loop reachability use the in-place version.
+"""
+function pke_th_rhs(x, t, plant, Q_sg, u)
+    dx = similar(x)
+    pke_th_rhs!(dx, x, t, plant, Q_sg, u)
+    return dx
+end