suslu-soken-2026.jl

# Süslü & Söken (2026), J. Spacecraft & Rockets — Article in Advance
# "Convex Optimization of Lunar Ascent Trajectory with Generalized
#  Terminal Conditions" — DOI: 10.2514/1.A36519
#
# Free-time sequential convex programming (SCvx) replication of Sec. IV
# with both terminal-condition formulations:
#   • CTC (Circular Terminal Conditions, eqs. 9–12)  — 4 quadratic eqs.
#   • GTC (Generalized Terminal Conditions, eqs. 16–23) — 8 vector eqs.
#
# Scenarios (Table 1 vehicle: T_max = 10 kN, I_sp = 300 s,
#                              m_dry = 400 kg, m_prop = 1600 kg):
#   1. Inclined circular (a = R_M + 200 km, e = 0, i = 60°, Ω = 280°)
#                                                — GTC + CTC
#   2. Equatorial circular (i = 0°, launch lat = 60° N)
#                                                — GTC + CTC (paper: CTC fails)
#   3. Elliptic (a = 3087.4 km, e = 0.372, i = 40°, Ω = 290°, ω = 90°)
#                                                — GTC + CTC (paper: not run)
#
# Paper-typo notes (checked against the printed equations and the
# appendix Jacobian H, eqs. A5–A42):
#   • Eqs. 16–20 are correct as printed:
#       g₁ = r_y v_z − r_z v_y = h_tx,   g₂ = r_z v_x − r_x v_z = h_ty,
#       g₃ = r_x v_y − r_y v_x = h_tz,
#       g₄ = r_x v_z − r_z v_x = N_tx,   g₅ = r_y v_z − r_z v_y = N_ty.
#     g₅ is algebraically equal to g₁ (since N = k̂ × h ⇒ N_y = h_x);
#     the paper keeps both, so we do too.
#   • Eq. 13 has a typesetting slip: the third component of h prints as
#     "v_y − r_y v_x", its leading r_x orphaned at the end of the second
#     row.  Harmless — eq. 18 gives the intended r_x v_y − r_y v_x.
#   • The eccentricity expansions (eqs. 15, 21–23) contain subscript and
#     sign slips (e.g. eq. 21's middle coefficient v_fy should be v_fz;
#     eqs. 22–23 print +μ r_f-terms where e = (1/μ)(v × h − μ r/‖r‖)
#     requires −).  The appendix H rows (A8 ff.) are consistent with the
#     standard e-vector, which terminal_residual() implements directly
#     via cross products, so none of these slips affect the code.
#   • Appendix A3 shows B[7,4] = +α; dynamics ż = −α u_N requires −α.
#
# Documented deviations from the paper:
#   • Trust-region update: the paper defers to the SCvx update rules of
#     ref. [43] (accept/reject on the ratio ρ of actual-to-predicted
#     reduction of an exact penalty function).  Both that ratio test
#     (tr_update=:ratio) and an accept-every-step stagnation heuristic
#     (tr_update=:stagnation) are implemented; the DEFAULT IS THE
#     HEURISTIC because the faithful ratio test converges too slowly
#     for the paper's 50-iteration budget on this problem — see the
#     scvx() docstring for the measured behavior behind that choice.
#   • Trust region on p: eq. 51 bounds only the control.  A relative
#     trust region on p (scaled with η) is kept because the free-time
#     subproblem is otherwise prone to artificial unboundedness in t_f.
#   • Discretization: eq. 54 is a first-order-hold on the control but
#     explicit Euler in the state, exactly as printed, and is the
#     default; its defect is quantified post-hoc by propagating the
#     converged control through the full nonlinear dynamics.
#     SUSLU_DISC=foh switches to a matrix-exponential FOH (Van Loan
#     block exponential per interval, A frozen at the left node,
#     residual defined against nonlinear propagation of the reference —
#     a multiple-shooting defect).  This is NOT what the paper
#     discretizes, but it is far more accurate at the same N: the
#     nonlinear-prop terminal error drops from km-scale to ≤ 0.03 km
#     across the converging cases.  Two practical notes: (a) the radius
#     in the frozen Jacobian is floored at 0.3 DU because intermediate
#     iterates can place reference nodes at nonphysical near-zero radii
#     where exp(A Δτ) overflows; (b) the Fig.-7 monotonic-altitude
#     variant does not converge under :foh from the default
#     initialization — use the default :euler for that comparison.
#
# Extras beyond the basic replication:
#   • Convergence additionally requires the virtual controls to vanish
#     ("virtual terms must be zero for a converged solution", Sec. III.D).
#   • Equatorial CTC failure is retried at N = 50 (Sec. IV.C).
#   • SUSLU_MC=<n> runs an n-sample Monte Carlo of Table 4/5
#     (paper: 500 samples, GTC 48.2 % / CTC 0 % convergence).
#   • The Fig. 7 lofted-vs-monotonic-altitude comparison is reproduced
#     via a linearized altitude-nondecrease constraint.
#
# Internal scaling: canonical units with DU = R_M and μ = 1 for
# numerical conditioning.  SI is used only for input/output and plots.
#
# Run with:
#   julia --project=. suslu-soken-2026.jl
# Env switches: SUSLU_SOLVER=ECOS|Clarabel, SUSLU_DISC=euler|foh,
#               SUSLU_MC=<n samples>

using JuMP
using ECOS
using Clarabel
using LinearAlgebra
using OrdinaryDiffEq
using Printf
using Random

# Solver selection: set to ECOS.Optimizer or Clarabel.Optimizer.
const SOLVER = get(ENV, "SUSLU_SOLVER", "Clarabel") == "ECOS" ?
               ECOS.Optimizer : Clarabel.Optimizer
const SOLVER_NAME = SOLVER === ECOS.Optimizer ? "ECOS" : "Clarabel"

# Discretization selection: "euler" = paper's eq. 54 (default),
# "foh" = matrix-exponential FOH (documented deviation, see header).
const DISCRETIZATION = let d = Symbol(get(ENV, "SUSLU_DISC", "euler"))
    d in (:euler, :foh) || error("SUSLU_DISC must be 'euler' or 'foh'")
    d
end

# Trust-region update selection (see scvx docstring for why the
# accept-always heuristic is the default in both discretizations).
const TR_UPDATE = let t = Symbol(get(ENV, "SUSLU_TR", "stagnation"))
    t in (:ratio, :stagnation) || error("SUSLU_TR must be 'ratio' or 'stagnation'")
    t
end
ENV["GKSwstype"] = "100"
using Plots

# ─── Physical & vehicle constants (SI) ────────────────────────────────
const G0       = 9.80665                  # m/s²
const MU_SI    = 4.902800066e12           # m³/s² (selenocentric)
const R_MOON   = 1.7374e6                 # m
const OMEGA_M  = 2.6617e-6                # rad/s (sidereal rotation about z)

const T_MAX_SI = 10_000.0                 # N
const I_SP     = 300.0                    # s
const M_DRY    = 400.0                    # kg
const M_PROP   = 1600.0                   # kg
const M_WET    = M_DRY + M_PROP           # 2000 kg
const ALPHA_SI = 1 / (I_SP * G0)          # s/m, mass-flow coefficient

# ─── Canonical units (μ = 1) ──────────────────────────────────────────
const DU = R_MOON
const TU = sqrt(DU^3 / MU_SI)             # ≈ 1059 s
const VU = DU / TU                        # ≈ 1640 m/s
const AU = DU / TU^2                      # ≈ 1.55 m/s²

const ALPHA  = ALPHA_SI * VU              # canonical α
const T_MAX  = T_MAX_SI / AU              # canonical thrust-magnitude (kg)
                                          # u_N ≤ T_MAX · exp(−z) with z = ln m_kg
const OMEGA  = OMEGA_M * TU               # canonical moon angular velocity

const Z_WET = log(M_WET)
const Z_DRY = log(M_DRY)

# ─── Orbital mechanics helpers (canonical units) ─────────────────────
function R_z(θ); c, s = cos(θ), sin(θ); return [c -s 0.0; s c 0.0; 0.0 0.0 1.0]; end
function R_x(θ); c, s = cos(θ), sin(θ); return [1.0 0.0 0.0; 0.0 c -s; 0.0 s c]; end

"""Classical elements → (r, v) in MCI. a in canonical DU, μ = 1, angles in rad."""
function elements_to_rv(a, e, i, Ω, ω, ν)
    p = a * (1 - e^2)
    r_mag = p / (1 + e*cos(ν))
    r_pf  = r_mag * [cos(ν), sin(ν), 0.0]
    v_pf  = sqrt(1/p) * [-sin(ν), e + cos(ν), 0.0]
    Q     = R_z(Ω) * R_x(i) * R_z(ω)
    return Q * r_pf, Q * v_pf
end

"""Target vectors (h, N, e) from a representative (r, v) on the target orbit."""
function target_vectors(r, v)
    h = cross(r, v)
    N = cross([0.0, 0.0, 1.0], h)
    e = cross(v, h) - r / norm(r)
    return h, N, e
end

"""(r, v) → classical orbital elements (a, e, i, Ω, ω) in canonical, μ = 1."""
function rv_to_elements(r, v)
    hv = cross(r, v);  h = norm(hv)
    Nv = cross([0.0, 0.0, 1.0], hv);  Nm = norm(Nv)
    ev = cross(v, hv) - r/norm(r);  e = norm(ev)
    a  = 1 / (2/norm(r) - dot(v, v))
    i  = acos(clamp(hv[3]/h, -1, 1))
    if Nm > 1e-10
        Ω = acos(clamp(Nv[1]/Nm, -1, 1))
        if Nv[2] < 0; Ω = 2π - Ω; end
    else
        Ω = 0.0
    end
    if Nm > 1e-10 && e > 1e-8
        ω = acos(clamp(dot(Nv, ev) / (Nm*e), -1, 1))
        if ev[3] < 0; ω = 2π - ω; end
    else
        ω = 0.0
    end
    return (a=a, e=e, i=i, Ω=Ω, ω=ω)
end

"""Launch point in MCI (= MCMF at t=0): from latitude, longitude (deg), altitude (m).
Returns r (canonical) and v = ω_M × r (canonical, including moon-rotation)."""
function launch_state(lat_deg, lon_deg, alt_m)
    λ = deg2rad(lat_deg); ϕ = deg2rad(lon_deg)
    r_mag = (R_MOON + alt_m) / DU
    r = r_mag * [cos(λ)*cos(ϕ), cos(λ)*sin(ϕ), sin(λ)]
    v = cross([0.0, 0.0, OMEGA], r)
    return r, v
end

"""MCI position → MCMF latitude/longitude (deg) after canonical time t."""
function mci_to_latlon(r, t_canonical)
    θ = OMEGA * t_canonical
    c, s = cos(-θ), sin(-θ)
    r_mcmf = [c*r[1] - s*r[2], s*r[1] + c*r[2], r[3]]
    lat = asind(clamp(r_mcmf[3] / norm(r_mcmf), -1, 1))
    lon = atand(r_mcmf[2], r_mcmf[1])
    return lat, lon
end

# ─── Reference trajectory: gravity-turn integration ──────────────────
"""
Vern6 integration of (r, v, z) under full-thrust gravity-turn.
For t < t_kick the thrust direction is fixed at the initial pitched-east
unit vector T̂₀; afterward it tracks v̂ (the classical gravity-turn
maneuver).  When propellant is exhausted the vehicle coasts; the burnout
switch is localized by a ContinuousCallback on z − z_dry.

The kick duration is an initialization knob, not a paper parameter
(Table 4 samples only pitch, azimuth, and t_f).  The default is 30 s:
with the integrator resolving the kick exactly, a 10 s kick leaves the
nominal equatorial case feasibility-converged but creeping along the
free-t_f direction too slowly to meet eq. 52 within 50 iterations.
(The former hand-rolled RK4 at ~12 s steps never resolved a 10 s kick
anyway, so its effective initialization was different.)

Returns:  x̄ (7×(N+1) reference state),  ū (3×(N+1) reference thrust
acceleration u = T/m in canonical units),  p̄ (canonical t_f).
"""
function gravity_turn_reference(lat_deg, lon_deg, alt_m, t_f_s, N_intervals;
                                pitch_deg=70.0, azimuth_deg=90.0,
                                t_kick_s=30.0)
    r0, v0 = launch_state(lat_deg, lon_deg, alt_m)
    up    = r0 / norm(r0)
    ϕ     = deg2rad(lon_deg); λ = deg2rad(lat_deg)
    east  = [-sin(ϕ), cos(ϕ), 0.0]
    north = cross(up, east)
    az    = deg2rad(azimuth_deg)
    horiz = sin(az)*east + cos(az)*north
    pitch = deg2rad(pitch_deg)
    T̂0    = cos(pitch)*horiz + sin(pitch)*up

    p̄ = t_f_s / TU

    function thrust_dir(v, t_can)
        if t_can < t_kick_s / TU
            return T̂0
        elseif norm(v) > 1e-6
            return v / norm(v)
        else
            return T̂0
        end
    end

    function thrust_accel(x, t_can)
        v = x[4:6]; z = x[7]
        if z > Z_DRY + 1e-6
            u_mag = T_MAX / exp(z)
            return u_mag * thrust_dir(v, t_can), u_mag
        end
        return zeros(3), 0.0
    end

    function rhs!(dx, x, _, t_can)
        u, u_N = thrust_accel(x, t_can)
        dx[1:3] = x[4:6]
        dx[4:6] = u .- x[1:3] ./ norm(x[1:3])^3
        dx[7]   = -ALPHA * u_N
    end

    ts = range(0.0, p̄; length=N_intervals+1)
    # No-op event at burnout so the step lands exactly on the thrust
    # cutoff instead of integrating across the rhs discontinuity.
    burnout = ContinuousCallback((x, t, integ) -> x[7] - (Z_DRY + 1e-6),
                                 nothing; affect_neg! = integ -> nothing,
                                 save_positions = (false, false))
    prob = ODEProblem(rhs!, [r0; v0; Z_WET], (0.0, p̄))
    sol = solve(prob, Vern6(); reltol=1e-10, abstol=1e-10, saveat=ts,
                tstops=[t_kick_s / TU], callback=burnout)
    x̄ = Array(sol)
    x̄[7, :] .= max.(x̄[7, :], Z_DRY)   # floor for stability

    ū = zeros(3, N_intervals+1)
    for k in 1:N_intervals+1
        ū[:, k], _ = thrust_accel(x̄[:, k], ts[k])
    end
    return x̄, ū, p̄
end

# ─── Terminal conditions: nonlinear residual and linearization ────────
"""
Nonlinear terminal residual g(x_f) (zero ⟺ target orbit reached).

mode = :CTC → eqs. 9–12 (4 quadratic equalities)
mode = :GTC → eqs. 16–23 (8 equalities; the e-vector rows use the
standard e = v × h − r/‖r‖ rather than the garbled printed
expansions — see the file header).
"""
function terminal_residual(x_f, h_t, N_t, e_t, mode)
    r = x_f[1:3]; v = x_f[4:6]
    if mode == :CTC
        # For circular target: r_t = |h_t|² (μ = 1, since |h| = √(μ r) ⇒ r = h²),
        # cos i_t = h_tz / |h_t|.
        ht_mag = norm(h_t)
        r_t = ht_mag^2
        cos_it = h_t[3] / ht_mag
        return [dot(r, r) - r_t^2,                              # eq. 9
                dot(v, v) - 1/r_t,                              # eq. 10, μ = 1
                dot(r, v),                                      # eq. 11
                r[1]*v[2] - r[2]*v[1] - sqrt(r_t)*cos_it]       # eq. 12
    else  # :GTC — note N = k̂ × h = (−h_y, h_x, 0)
        h  = cross(r, v)
        ev = cross(v, h) - r / norm(r)
        return [h[1] - h_t[1],                                  # g1 = h_x
                h[2] - h_t[2],                                  # g2 = h_y
                h[3] - h_t[3],                                  # g3 = h_z
                -h[2] - N_t[1],                                 # g4 = N_x
                h[1] - N_t[2],                                  # g5 = N_y (= g1)
                ev[1] - e_t[1],
                ev[2] - e_t[2],
                ev[3] - e_t[3]]
    end
end

"""
Build (H, l) such that  H · x_f + l = 0  is the first-order linearization
of the terminal conditions at the reference final state x̄_f (eq. 46).
"""
function terminal_linearization(x̄_f, h_t, N_t, e_t, mode)
    r̄ = x̄_f[1:3]; v̄ = x̄_f[4:6]
    g = terminal_residual(x̄_f, h_t, N_t, e_t, mode)
    if mode == :CTC
        H = zeros(4, 7)
        # g1 = r·r − r_t²   (eq. 9)
        H[1, 1:3] = 2 .* r̄
        # g2 = v·v − 1/r_t   (eq. 10, μ = 1)
        H[2, 4:6] = 2 .* v̄
        # g3 = r·v            (eq. 11)
        H[3, 1:3] = v̄; H[3, 4:6] = r̄
        # g4 = r_x v_y − r_y v_x − √r_t · cos i_t   (eq. 12)
        H[4, 1] =  v̄[2]; H[4, 2] = -v̄[1]
        H[4, 4] = -r̄[2]; H[4, 5] =  r̄[1]
    else  # :GTC — eight constraints, N = k̂ × h = (−h_y, h_x, 0)
        H = zeros(8, 7)
        set!(row, dr, dv) = (H[row, 1:3] = dr; H[row, 4:6] = dv)
        # g1: r_y v_z − r_z v_y − h_tx
        set!(1, [0.0, v̄[3], -v̄[2]], [0.0, -r̄[3], r̄[2]])
        # g2: r_z v_x − r_x v_z − h_ty
        set!(2, [-v̄[3], 0.0, v̄[1]], [r̄[3], 0.0, -r̄[1]])
        # g3: r_x v_y − r_y v_x − h_tz
        set!(3, [v̄[2], -v̄[1], 0.0], [-r̄[2], r̄[1], 0.0])
        # g4: r_x v_z − r_z v_x − N_tx
        set!(4, [v̄[3], 0.0, -v̄[1]], [-r̄[3], 0.0, r̄[1]])
        # g5: r_y v_z − r_z v_y − N_ty   (same expression as g1)
        set!(5, [0.0, v̄[3], -v̄[2]], [0.0, -r̄[3], r̄[2]])
        # g6–g8: eccentricity-vector components.
        # Using bac-cab,  e = r(v·v) − v(r·v) − r/‖r‖   (μ = 1), so
        #   ∂e_i/∂r_j = δ_ij (v² − 1/‖r‖) − v_i v_j + r_i r_j /‖r‖³
        #   ∂e_i/∂v_j = 2 r_i v_j − δ_ij (r·v) − v_i r_j
        # (matches appendix eqs. A8–A23, A27–A42 elementwise).
        rmag = norm(r̄); rv = dot(r̄, v̄); v2 = dot(v̄, v̄)
        for i in 1:3
            drow = [((i==j) ? (v2 - 1/rmag) : 0.0) - v̄[i]*v̄[j] + r̄[i]*r̄[j]/rmag^3 for j in 1:3]
            vrow = [2*r̄[i]*v̄[j] - ((i==j) ? rv : 0.0) - v̄[i]*r̄[j] for j in 1:3]
            set!(5+i, drow, vrow)
        end
    end
    l = g - H * x̄_f
    return H, l
end

# ─── State Jacobian J_x = ∂f/∂x for the linearized dynamics ──────────
function jac_x(r̄)
    r = norm(r̄)
    J = zeros(7, 7)
    J[1:3, 4:6] = I(3)                          # ṙ = v
    J[4:6, 1:3] = -I(3)/r^3 + 3 .* (r̄*r̄') / r^5   # ∂v̇/∂r = ∂(−r/r³)/∂r
    # ż depends on the slack u_N only — no state derivative
    return J
end

# ─── Nonlinear propagation over one normalized-time interval ─────────
"""
Vern6 integration of the nonlinear dynamics f̂ = p·f over one Δτ interval,
with FOH interpolation of (u, u_N) between the endpoint controls.
"""
function propagate_interval(x0, u0, uN0, u1, uN1, p, Δτ)
    T = p * Δτ
    function f!(dx, x, _, t)
        s  = t / T
        u  = (1 - s) .* u0 .+ s .* u1
        uN = (1 - s)  * uN0 +  s * uN1
        dx[1:3] = x[4:6]
        dx[4:6] = u .- x[1:3] ./ norm(x[1:3])^3
        dx[7]   = -ALPHA * uN
    end
    prob = ODEProblem(f!, copy(x0), (0.0, T))
    sol = solve(prob, Vern6(); reltol=1e-10, abstol=1e-10,
                save_everystep=false, save_start=false)
    return sol.u[end]
end

# ─── Matrix-exponential FOH discretization (deviation, see header) ───
"""
Exact first-order-hold discretization of the linearized dynamics over
one interval, with A frozen at the interval's left node:

  x_{k+1} = A_d x_k + B⁻ (u_k; u_Nk) + B⁺ (u_{k+1}; u_Nk+1) + F_d p + r_d

The Van Loan block exponential of  G = [A I 0; 0 0 I; 0 0 0]·Δτ  yields
A_d = e^{AΔτ} together with the integrals
  I₀ = ∫₀^Δτ e^{A(Δτ−s)} ds   and   I₁ = ∫₀^Δτ e^{A(Δτ−s)} s ds,
from which B⁻ = (I₀ − I₁/Δτ)B and B⁺ = (I₁/Δτ)B.  r_d is chosen so the
model reproduces the *nonlinear* propagation of the reference across
the interval, i.e. the virtual control measures a multiple-shooting
defect instead of an Euler defect.
"""
function foh_discretize(x̄k, ūk, ūNk, ūk1, ūNk1, p̄, Δτ)
    n = 7
    # Floor the radius used in the frozen Jacobian: intermediate SCP
    # iterates can place reference nodes at nonphysical near-zero radii,
    # where exp(A Δτ) ~ e^{p̄Δτ/r³} overflows (the Euler model degrades
    # only polynomially there).  0.3 DU is far below any physical state.
    r̄k = x̄k[1:3]
    if norm(r̄k) < 0.3
        r̄k = r̄k .* (0.3 / norm(r̄k))
    end
    A_c = p̄ .* jac_x(r̄k)
    B_c = p̄ .* J_U
    F_c = [x̄k[4:6]; ūk .- x̄k[1:3] ./ norm(x̄k[1:3])^3; -ALPHA * ūNk]
    G = zeros(3n, 3n)
    G[1:n, 1:n]        = A_c
    G[1:n, n+1:2n]     = I(n)
    G[n+1:2n, 2n+1:3n] = I(n)
    E  = exp(G .* Δτ)
    A_d = E[1:n, 1:n]
    I0  = E[1:n, n+1:2n]
    I1  = E[1:n, 2n+1:3n]
    Bm  = (I0 .- I1 ./ Δτ) * B_c
    Bp  = (I1 ./ Δτ) * B_c
    F_d = I0 * F_c
    xprop = propagate_interval(x̄k, ūk, ūNk, ūk1, ūNk1, p̄, Δτ)
    r_d = xprop - A_d * x̄k - Bm * [ūk; ūNk] - Bp * [ūk1; ūNk1] - F_d * p̄
    return A_d, Bm, Bp, F_d, r_d
end

# Control Jacobian (constant): ∂f/∂[u; u_N] = [0 0; I 0; 0 −α]
const J_U = begin
    M = zeros(7, 4)
    M[4:6, 1:3] = I(3)
    M[7, 4] = -ALPHA
    M
end

# ─── SCvx subproblem ────────────────────────────────────────────────
"""
Solve one SCvx subproblem at reference (x̄, ū, ūN, p̄).

LCvx (eqs. 39, 48) replaces ‖T‖ = u_N with ‖u‖ ≤ u_N and linearizes the
upper bound u_N ≤ T_max e^(−z̄)(1 − (z − z̄))  +  ν_s.
Dynamics use eq. 54 (FOH average-of-endpoints, explicit Euler in the
state) by default; discretization=:foh switches to the
matrix-exponential FOH of foh_discretize() — a documented deviation.
Trust region (eq. 51) caps ‖u − ū‖∞ ≤ η.  Terminal constraint
H x_N + l + ν_f = 0 (eq. 49).  Objective per eq. 55 (mass-max + L1
virtual-control penalty).

monotonic_alt=true adds the Sec. IV.C / Fig. 7 variant: altitude is
forced nondecreasing via the linearized radial constraint
r̄ₖ·(r_{k+1} − r_k) ≥ 0.
"""
function solve_subproblem(x̄, ū, ūN, p̄, h_t, N_t, e_t, mode;
                          η=1.0, η_p=0.50, λ=1e3, fix_p=false,
                          p_max=Inf, monotonic_alt=false,
                          discretization=:euler)
    Nn = size(x̄, 2)      # number of nodes
    Ni = Nn - 1           # intervals
    Δτ = 1.0 / Ni

    model = Model(SOLVER)
    set_silent(model)
    if SOLVER === ECOS.Optimizer
        set_optimizer_attribute(model, "maxit", 200)
    end

    @variable(model, x[1:7, 1:Nn])
    @variable(model, u[1:3, 1:Nn])
    @variable(model, uN[1:Nn] >= 0)
    @variable(model, p >= 0)
    if fix_p
        @constraint(model, p == p̄)
    else
        # Per-iter trust region on p (paper's eq. 51 omits this, but it's
        # needed in practice to avoid artificial unboundedness in t_f).
        @constraint(model, p - p̄ <=  η_p * p̄)
        @constraint(model, p̄ - p <=  η_p * p̄)
        # Global ceiling — prevents runaway across many iterations.
        if isfinite(p_max)
            @constraint(model, p <= p_max)
        end
    end

    # Virtual controls + L1 epigraphs
    @variable(model, ν[1:7, 1:Ni])
    @variable(model, νs[1:Nn] >= 0)
    n_term = (mode == :CTC) ? 4 : 8
    @variable(model, νf[1:n_term])
    @variable(model, abs_ν[1:7, 1:Ni]  >= 0)
    @variable(model, abs_νf[1:n_term]  >= 0)
    @constraint(model, [i=1:7, k=1:Ni],     ν[i, k] <= abs_ν[i, k])
    @constraint(model, [i=1:7, k=1:Ni],    -ν[i, k] <= abs_ν[i, k])
    @constraint(model, [i=1:n_term],     νf[i] <= abs_νf[i])
    @constraint(model, [i=1:n_term],    -νf[i] <= abs_νf[i])

    # Initial state (eq. 6–8, lifted via launch_state into the reference)
    @constraint(model, x[:, 1] .== x̄[:, 1])

    # Dynamics linearization.  Both forms place ν inside a Δτ factor so
    # ν has the same per-unit-τ meaning (and penalty weight) in each.
    for k in 1:Ni
        if discretization == :foh
            # Matrix-exponential FOH (deviation from eq. 54, see header)
            A_d, Bm, Bp, F_d, r_d = foh_discretize(
                x̄[:, k], ū[:, k], ūN[k], ū[:, k+1], ūN[k+1], p̄, Δτ)
            rhs = A_d * x[:, k] + Bm * [u[:, k]; uN[k]] +
                  Bp * [u[:, k+1]; uN[k+1]] + F_d * p + r_d + Δτ .* ν[:, k]
            @constraint(model, x[:, k+1] .== rhs)
        else
            # eq. 54: explicit Euler in the state, FOH average controls
            rbar = x̄[1:3, k]; vbar = x̄[4:6, k]
            u_ref = ū[:, k];  uN_ref = ūN[k]
            Ak = p̄ .* jac_x(rbar)
            Bk = p̄ .* J_U
            # f̂_k at the reference, which equals F_k · p̄ (eq. A4, with α-sign fix)
            Fk = [vbar;
                  u_ref - rbar / norm(rbar)^3;
                  -ALPHA * uN_ref]
            # r_k = f̂(x̄,ū,p̄) − A x̄ − B ū − F p̄ = −A x̄ − B ū  (cancels with f̂)
            rk = -Ak * x̄[:, k] - Bk * [u_ref; uN_ref]
            u_avg = 0.5 * ([u[:, k]; uN[k]] + [u[:, k+1]; uN[k+1]])
            rhs = x[:, k] + Δτ * (Ak * x[:, k] + Bk * u_avg + Fk * p + rk + ν[:, k])
            @constraint(model, x[:, k+1] .== rhs)
        end
    end

    # Thrust constraints (eqs. 39, 48) at every node
    for k in 1:Nn
        @constraint(model, [uN[k]; u[:, k]] in SecondOrderCone())
        zbar = x̄[7, k]
        @constraint(model,
            uN[k] - T_MAX * exp(-zbar) * (1 - (x[7, k] - zbar)) <= νs[k])
    end

    # Mass floor (eq. 36)
    @constraint(model, x[7, Nn] >= Z_DRY)

    # Trust region (eq. 51) on control vector (u; uN)
    for k in 1:Nn
        @constraint(model,
            [η; u[:, k] .- ū[:, k]; uN[k] - ūN[k]]
            in MOI.NormInfinityCone(5))
    end

    # Optional Fig.-7 variant: altitude nondecreasing, linearized as
    # d‖r‖ ≈ r̄·dr / ‖r̄‖ ≥ 0 along the reference.
    if monotonic_alt
        for k in 1:Ni
            rbar = x̄[1:3, k]
            @constraint(model, dot(rbar, x[1:3, k+1] .- x[1:3, k]) >= 0)
        end
    end

    # Terminal linearization (eq. 49)
    H, l = terminal_linearization(x̄[:, Nn], h_t, N_t, e_t, mode)
    @constraint(model, H * x[:, Nn] + l + νf .== 0)

    # Objective — exactly eq. 55: J = −z_N + λ‖ν_f‖₁ + λ Σ(‖νᵢ‖₁ + ‖ν_si‖₁).
    # (No Δτ on the sum; eq. 55 drops it when discretizing eq. 50.)
    @objective(model, Min,
        -x[7, Nn]
        + λ * sum(abs_νf)
        + λ * (sum(abs_ν) + sum(νs)))

    optimize!(model)
    st = termination_status(model)
    ok = st in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED, MOI.ALMOST_OPTIMAL)
    if !ok
        return (ok=false, status=st)
    end
    return (ok=true, status=st,
            x = value.(x), u = value.(u), uN = value.(uN), p = value(p),
            obj = objective_value(model),
            ν_max  = maximum(abs.(value.(ν))),
            νs_max = maximum(value.(νs)),
            νf_max = maximum(abs.(value.(νf))))
end

# ─── Exact penalty function for the SCvx ratio test ──────────────────
"""
Nonlinear penalized cost matching eq. 55 with the virtual controls
replaced by the actual constraint violations:
  • dynamics defect, matched to the active discretization:
    :euler  ν_i = (x_{i+1} − x_i − Δτ p f(x_i, ū_avg)) / Δτ  (ν sits
            inside the Δτ bracket in eq. 54, hence the division);
    :foh    ν_i = (x_{i+1} − Φ(x_i, u_i, u_{i+1}, p)) / Δτ, the
            multiple-shooting defect against nonlinear propagation;
  • LCvx upper-bound violation  ν_si = max(0, u_N − T_max e^{−z});
  • terminal residual ν_f = g(x_N).
At a reference point this equals the linearized subproblem objective,
which is what makes the ratio test well-posed.
"""
function penalized_cost(x, u, uN, p, h_t, N_t, e_t, mode, λ;
                        discretization=:euler)
    Nn = size(x, 2); Ni = Nn - 1; Δτ = 1.0 / Ni
    J = -x[7, Nn]
    for k in 1:Ni
        if discretization == :foh
            xprop = propagate_interval(x[:, k], u[:, k], uN[k],
                                       u[:, k+1], uN[k+1], p, Δτ)
            J += (λ / Δτ) * sum(abs, x[:, k+1] .- xprop)
        else
            u_avg  = 0.5 .* (u[:, k] .+ u[:, k+1])
            uN_avg = 0.5  * (uN[k]   +  uN[k+1])
            rk = x[1:3, k]
            f = [x[4:6, k]; u_avg .- rk ./ norm(rk)^3; -ALPHA * uN_avg]
            J += (λ / Δτ) * sum(abs, x[:, k+1] .- x[:, k] .- (Δτ * p) .* f)
        end
    end
    for k in 1:Nn
        J += λ * max(0.0, uN[k] - T_MAX * exp(-x[7, k]))
    end
    J += λ * sum(abs, terminal_residual(x[:, Nn], h_t, N_t, e_t, mode))
    return J
end

# ─── SCvx outer loop ────────────────────────────────────────────────
"""
SCvx outer loop.  Convergence (eq. 52) additionally requires the
virtual controls to vanish — "virtual terms must be zero for a
converged solution" (Sec. III.D).

Two trust-region update rules are available via `tr_update`:

  :ratio       The SCvx rule of ref. [43] that the paper cites:
               ρ = (actual reduction of the exact penalty) /
                   (predicted reduction by the convex subproblem);
               reject the step and shrink η when ρ < ρ₀, grow η when
               ρ ≥ ρ₂.

  :stagnation  (default) Accept every step; starting from η = 1, halve
               η (floor 0.12) whenever Δx fails to decrease for three
               consecutive iterations.

DEVIATION FROM THE PAPER: the default is :stagnation, not the cited
ratio rule.  On this problem the exact-penalty ratio test measured
ρ ≈ 0.3–0.5 almost everywhere (the quadratic terminal constraints make
the linear model consistently over-optimistic), so η shrinks until
feasibility progress — rate-limited by the per-node control trust
region — no longer fits in the paper's 50-iteration budget; with an η
floor it instead deadlocks re-solving an identical rejected subproblem.
Accepting every step exploits the Gauss-Newton character of the
re-linearized terminal constraints and reproduces the paper's iteration
counts.  ρ is still computed and logged in both modes.
"""
function scvx(x̄, ū, ūN, p̄, h_t, N_t, e_t, mode;
              max_iter=50, ε=3e-3, feas_tol=1e-6, η=1.0, η_p=0.50, λ=1e3,
              fix_p=false, p_max=Inf, monotonic_alt=false,
              tr_update=:stagnation, discretization=:euler, verbose=true)
    # :ratio parameters (standard SCvx values).  No useful η floor
    # exists: a floor deadlocks on a rejected step (same reference and
    # same η reproduce the same candidate), so it is kept tiny.
    ρ_0, ρ_1, ρ_2 = 0.0, 0.25, 0.7
    β_shrink, β_grow = 2.0, 2.0
    η_min, η_max = 1e-4, 10.0
    # :stagnation parameters.
    η_floor = 0.12
    stagnation = 0
    prev_Δx = Inf

    history = NamedTuple[]
    converged = false
    last = nothing
    η_cur = η
    J_ref = penalized_cost(x̄, ū, ūN, p̄, h_t, N_t, e_t, mode, λ;
                           discretization=discretization)
    for it in 1:max_iter
        ηp_cur = η_p * η_cur / η          # p trust region scales with η
        res = solve_subproblem(x̄, ū, ūN, p̄, h_t, N_t, e_t, mode;
                               η=η_cur, η_p=ηp_cur, λ=λ,
                               fix_p=fix_p, p_max=p_max,
                               monotonic_alt=monotonic_alt,
                               discretization=discretization)
        if !res.ok
            verbose && @printf("    iter %2d  solver FAILED: %s\n", it, res.status)
            push!(history, (iter=it, ok=false, status=res.status, accepted=false,
                            ρ=NaN, η=η_cur, Δx=NaN, Δp=NaN,
                            ν=NaN, νs=NaN, νf=NaN, p=p̄))
            return (ok=false, x=x̄, u=ū, p=p̄, last=last, history=history)
        end
        J_new  = penalized_cost(res.x, res.u, res.uN, res.p,
                                h_t, N_t, e_t, mode, λ;
                                discretization=discretization)
        pred   = J_ref - res.obj          # ≥ 0: reference is subproblem-feasible
        actual = J_ref - J_new
        ρ = pred > 1e-9 ? actual / pred : 1.0
        Δx = maximum(abs.(res.x .- x̄))
        Δp = abs(res.p - p̄)
        accepted = tr_update == :ratio ? ρ >= ρ_0 : true
        verbose && @printf("    iter %2d  %s ρ=%+.3f  η=%.3f  Δx=%.3e  Δp=%.3e  νmax=%.2e  νs=%.2e  νf=%.2e  p=%.4f\n",
                           it, accepted ? "acc" : "REJ", ρ, η_cur,
                           Δx, Δp, res.ν_max, res.νs_max, res.νf_max, res.p)
        push!(history, (iter=it, ok=true, status=res.status, accepted=accepted,
                        ρ=ρ, η=η_cur, Δx=Δx, Δp=Δp,
                        ν=res.ν_max, νs=res.νs_max, νf=res.νf_max, p=res.p))
        if accepted
            last = res
            x̄ = res.x; ū = res.u; ūN = res.uN; p̄ = res.p
            J_ref = J_new
            feasible = max(res.ν_max, res.νs_max, res.νf_max) <= feas_tol
            if Δx + Δp <= ε && feasible
                converged = true
                verbose && @printf("    converged at iter %d\n", it)
                break
            end
        end
        if tr_update == :ratio
            if ρ < ρ_1
                η_cur = max(η_cur / β_shrink, η_min)
            elseif ρ >= ρ_2
                η_cur = min(η_cur * β_grow, η_max)
            end
        else
            # Shrink trust region if Δx isn't improving (limit-cycle guard).
            if Δx > 0.9 * prev_Δx
                stagnation += 1
            else
                stagnation = 0
            end
            if stagnation >= 3 && η_cur > η_floor
                η_cur *= 0.5
                stagnation = 0
                verbose && @printf("    [trust region tightened to η=%.3f  η_p=%.3f]\n",
                                   η_cur, η_p * η_cur / η)
            end
            prev_Δx = Δx
        end
    end
    return (ok=converged, x=x̄, u=ū, p=p̄, last=last, history=history)
end

# ─── Post-convergence nonlinear validation ────────────────────────────
"""
Propagate the converged FOH control profile through the full nonlinear
dynamics (eqs. 1–3 in (u, z) form) with Vern6.  The SCP state satisfies
only the linearized, Euler-discretized dynamics of eq. 54, so this
quantifies the true discretization/linearization defect.
"""
function propagate_nonlinear(x0, u, uN, p)
    Nn = size(u, 2); Ni = Nn - 1; Δτ = 1.0 / Ni
    x = copy(x0)
    for k in 1:Ni
        x = propagate_interval(x, u[:, k], uN[k], u[:, k+1], uN[k+1], p, Δτ)
    end
    return x
end

# ─── Ballistic two-body propagation (Vern6) for visualizing the orbit ─
function propagate_twobody(r0, v0; T_span_can, N_pts=400)
    function f!(dx, x, _, _t)
        dx[1:3] = x[4:6]
        dx[4:6] = -x[1:3] ./ norm(x[1:3])^3
    end
    prob = ODEProblem(f!, [r0; v0], (0.0, T_span_can))
    sol = solve(prob, Vern6(); reltol=1e-10, abstol=1e-10,
                saveat=range(0.0, T_span_can; length=N_pts+1))
    return Array(sol)[1:3, :]
end

# ─── Scenario configurations ─────────────────────────────────────────
struct Scenario
    name      :: String
    lat0_deg  :: Float64
    lon0_deg  :: Float64
    alt0_m    :: Float64
    a_t_km    :: Float64
    e_t       :: Float64
    i_t_deg   :: Float64
    Ω_t_deg   :: Float64
    ω_t_deg   :: Float64
    t_f_s     :: Float64
    N_nodes   :: Int
    t_f_fixed :: Bool      # paper fixes t_f only for scenario 1
end

const SCENARIOS = [
    Scenario("inclined_circular",   0.0, -90.0, 0.0,
             R_MOON/1000 + 200,  0.0,    60.0, 280.0,  0.0, 1200.0, 200, true),
    Scenario("equatorial_circular", 60.0, -90.0, 0.0,
             R_MOON/1000 + 200,  0.0,     0.0,   0.0,  0.0, 2400.0, 200, false),
    Scenario("elliptic",            20.0, -90.0, 0.0,
             3087.4,             0.37248, 40.0, 290.0, 90.0, 2400.0, 200, false),
]

function run_case(sc::Scenario, mode::Symbol; max_iter=50, verbose=true,
                  pitch_deg=70.0, azimuth_deg=90.0, t_f_s=sc.t_f_s,
                  monotonic_alt=false, discretization=DISCRETIZATION,
                  tr_update=TR_UPDATE, label=sc.name)
    if verbose
        @printf("\n━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\n")
        @printf("  %-22s  [%s]\n", label, mode)
        @printf("───────────────────────────────────────────────────────────────\n")
    end
    a_can = sc.a_t_km * 1000 / DU
    r_t, v_t = elements_to_rv(a_can, sc.e_t,
        deg2rad(sc.i_t_deg), deg2rad(sc.Ω_t_deg), deg2rad(sc.ω_t_deg), 0.0)
    h_t, N_t, e_t_vec = target_vectors(r_t, v_t)
    verbose && @printf("  target  a=%.1f km  e=%.5f  i=%.2f°  Ω=%.2f°  ω=%.2f°\n",
        sc.a_t_km, sc.e_t, sc.i_t_deg, sc.Ω_t_deg, sc.ω_t_deg)
    x̄, ū, p̄ = gravity_turn_reference(sc.lat0_deg, sc.lon0_deg, sc.alt0_m,
                                      t_f_s, sc.N_nodes - 1;
                                      pitch_deg=pitch_deg, azimuth_deg=azimuth_deg)
    ūN = vec(sqrt.(sum(ū.^2; dims=1)))    # LCvx slack reference = ‖ū‖
    # Cap t_f at 1.5× the initial guess to suppress artificial unboundedness.
    p_max = sc.t_f_fixed ? Inf : 1.5 * p̄
    t0 = time()
    res = scvx(x̄, ū, ūN, p̄, h_t, N_t, e_t_vec, mode;
               max_iter=max_iter, fix_p=sc.t_f_fixed, p_max=p_max,
               monotonic_alt=monotonic_alt, discretization=discretization,
               tr_update=tr_update, verbose=verbose)
    elapsed = time() - t0
    if res.ok
        elems = rv_to_elements(res.x[1:3, end], res.x[4:6, end])
        m_final = exp(res.x[7, end])
        if verbose
            @printf("  → CONVERGED   t_f = %.1f s   m_f = %.2f kg   iters = %d   wall = %.2f s\n",
                res.p * TU, m_final, length(res.history), elapsed)
            @printf("    achieved   a=%.1f km  e=%.5f  i=%.3f°  Ω=%.3f°  ω=%.3f°\n",
                elems.a*DU/1000, elems.e, rad2deg(elems.i),
                rad2deg(elems.Ω), rad2deg(elems.ω))
            # Validate against the full nonlinear dynamics (the SCP state
            # only satisfies the eq.-54 linearized Euler discretization).
            x_nl = propagate_nonlinear(res.x[:, 1], res.last.u, res.last.uN, res.p)
            g_nl = terminal_residual(x_nl, h_t, N_t, e_t_vec, mode)
            e_nl = rv_to_elements(x_nl[1:3], x_nl[4:6])
            @printf("    nonlinear-prop check  ‖g‖∞=%.2e   a=%.1f km  e=%.5f  i=%.3f°  Ω=%.3f°  Δr_f=%.2f km\n",
                maximum(abs.(g_nl)), e_nl.a*DU/1000, e_nl.e,
                rad2deg(e_nl.i), rad2deg(e_nl.Ω),
                norm(x_nl[1:3] - res.x[1:3, end])*DU/1000)
        end
    elseif verbose
        n_iter = length(res.history)
        last_status = isempty(res.history) ? :unknown : res.history[end].status
        if res.last !== nothing
            @printf("  → DID NOT CONVERGE   iters = %d   last νf = %.2e   last νs = %.2e\n",
                n_iter, res.last.νf_max, res.last.νs_max)
            # report what the linearized terminal residual achieved
            elems = rv_to_elements(res.last.x[1:3, end], res.last.x[4:6, end])
            @printf("    last-iter approx orbit  a=%.1f km  e=%.4f  i=%.2f°  Ω=%.2f°\n",
                elems.a*DU/1000, elems.e, rad2deg(elems.i), rad2deg(elems.Ω))
        else
            @printf("  → FAILED at first iteration (%s)\n", last_status)
        end
    end
    return (sc=sc, mode=mode, res=res, h_t=h_t, elapsed=elapsed, label=label)
end

# ─── Plotting ────────────────────────────────────────────────────────
function plot_trajectory(case; outdir=".")
    res = case.res; sc = case.sc; label = case.label
    res.ok || return
    # res.last is the final SOCP subproblem solution (has uN, the LCvx slack);
    # res.x/u/p are the converged reference, which equals last.x/u/p modulo ε.
    lastsol = res.last
    x = lastsol.x; p = lastsol.p
    Nn = size(x, 2)
    ts_can = collect(0:Nn-1) ./ (Nn-1) .* p
    ts_si  = ts_can .* TU

    r_si = x[1:3, :] .* DU
    v_si = x[4:6, :] .* VU
    m_si = exp.(x[7, :])
    alt  = vec(sqrt.(sum(r_si.^2; dims=1))) .- R_MOON
    vmag = vec(sqrt.(sum(v_si.^2; dims=1)))
    # Recover thrust magnitude: T = m * ‖u_SI‖, where u_SI = u_canonical * AU
    u_si_mag = lastsol.uN .* AU            # since LCvx is tight, uN ≈ ‖u‖
    T_kN = m_si .* u_si_mag ./ 1000

    # Propagate the achieved orbit one period for the figure
    elems = rv_to_elements(x[1:3, end], x[4:6, end])
    a_can = elems.a
    T_orbit = 2π * sqrt(a_can^3)   # canonical, μ=1
    rs_prop = propagate_twobody(x[1:3, end], x[4:6, end];
                                T_span_can=T_orbit, N_pts=400) .* DU

    # 3D MCI trajectory + propagated orbit + moon sphere
    u_sph = range(0, 2π; length=40); v_sph = range(0, π; length=20)
    xs = [cos(u)*sin(v)*R_MOON/1000 for u in u_sph, v in v_sph]
    ys = [sin(u)*sin(v)*R_MOON/1000 for u in u_sph, v in v_sph]
    zs = [cos(v)     *R_MOON/1000 for u in u_sph, v in v_sph]
    p3d = surface(xs, ys, zs; color=:lightgray, alpha=0.35, label="",
                  xlabel="x_MCI [km]", ylabel="y_MCI [km]", zlabel="z_MCI [km]",
                  title="$(label) [$(case.mode)] — 3D MCI", legend=:topright)
    plot!(p3d, r_si[1, :]./1000, r_si[2, :]./1000, r_si[3, :]./1000;
          color=:blue, linewidth=2.2, label="Ascent")
    plot!(p3d, rs_prop[1, :]./1000, rs_prop[2, :]./1000, rs_prop[3, :]./1000;
          color=:blue, linestyle=:dash, label="Propagated orbit")
    savefig(p3d, joinpath(outdir, "$(label)_$(case.mode)_3d.png"))

    # Ground track
    lats = Float64[]; lons = Float64[]
    for k in 1:Nn
        la, lo = mci_to_latlon(x[1:3, k], ts_can[k])
        push!(lats, la); push!(lons, lo)
    end
    # propagated ground track: time continues past the burn
    lats_p = Float64[]; lons_p = Float64[]
    Nprop = size(rs_prop, 2)
    ts_p_can = ts_can[end] .+ collect(0:Nprop-1) ./ (Nprop-1) .* T_orbit
    for k in 1:Nprop
        la, lo = mci_to_latlon(rs_prop[1:3, k] ./ DU, ts_p_can[k])
        push!(lats_p, la); push!(lons_p, lo)
    end
    p_gt = plot(lons, lats; color=:blue, linewidth=2.2, label="Ascent",
                xlabel="Longitude [°]", ylabel="Latitude [°]",
                xlims=(-180, 180), ylims=(-90, 90),
                title="$(label) [$(case.mode)] — ground track",
                legend=:topright)
    scatter!(p_gt, [sc.lon0_deg], [sc.lat0_deg];
             color=:black, marker=:star, label="Launch")
    # Drawn as dots rather than a line so the ±180° wrap doesn't streak
    scatter!(p_gt, lons_p, lats_p; color=:blue, label="Propagated",
             markersize=1, markerstrokewidth=0)
    savefig(p_gt, joinpath(outdir, "$(label)_$(case.mode)_groundtrack.png"))

    # Profiles (altitude, velocity, thrust, mass)
    p1 = plot(ts_si, alt./1000; xlabel="Time [s]", ylabel="Altitude [km]",
              color=:blue, linewidth=2, label="")
    p2 = plot(ts_si, vmag./1000; xlabel="Time [s]", ylabel="Velocity [km/s]",
              color=:blue, linewidth=2, label="")
    p3 = plot(ts_si, T_kN; xlabel="Time [s]", ylabel="Thrust [kN]",
              color=:blue, linewidth=2, label="", ylims=(-0.5, 11))
    p4 = plot(ts_si, m_si; xlabel="Time [s]", ylabel="Mass [kg]",
              color=:blue, linewidth=2, label="")
    p_prof = plot(p1, p2, p3, p4; layout=(2, 2), size=(900, 600),
                  plot_title="$(label) [$(case.mode)] — profiles")
    savefig(p_prof, joinpath(outdir, "$(label)_$(case.mode)_profiles.png"))
    @printf("    saved 3D, ground track, and profile plots\n")
end

"""Fig.-7-style overlay: lofted vs monotonically increasing altitude."""
function plot_fig7(case_lofted, case_mono; outdir=".")
    (case_lofted.res.ok && case_mono.res.ok) || return
    function series(case)
        x = case.res.last.x; p = case.res.last.p
        Nn = size(x, 2)
        ts = collect(0:Nn-1) ./ (Nn-1) .* p .* TU
        alt = vec(sqrt.(sum((x[1:3, :] .* DU).^2; dims=1))) .- R_MOON
        vmag = vec(sqrt.(sum((x[4:6, :] .* VU).^2; dims=1)))
        m = exp.(x[7, :])
        T_kN = m .* case.res.last.uN .* AU ./ 1000
        return ts, alt, vmag, m, T_kN
    end
    t1, a1, v1, m1, T1 = series(case_lofted)
    t2, a2, v2, m2, T2 = series(case_mono)
    lab1 = "Lofted"; lab2 = "Monotonically increasing"
    p1 = plot(t1, a1./1000; label=lab1, color=:blue, linewidth=2,
              xlabel="Time [s]", ylabel="Altitude [km]")
    plot!(p1, t2, a2./1000; label=lab2, color=:red, linewidth=2)
    p2 = plot(t1, v1./1000; label=lab1, color=:blue, linewidth=2,
              xlabel="Time [s]", ylabel="Velocity [km/s]")
    plot!(p2, t2, v2./1000; label=lab2, color=:red, linewidth=2)
    p3 = plot(t1, T1; label=lab1, color=:blue, linewidth=2,
              xlabel="Time [s]", ylabel="Thrust [kN]", ylims=(-0.5, 11))
    plot!(p3, t2, T2; label=lab2, color=:red, linewidth=2)
    p4 = plot(t1, m1; label=lab1, color=:blue, linewidth=2,
              xlabel="Time [s]", ylabel="Mass [kg]")
    plot!(p4, t2, m2; label=lab2, color=:red, linewidth=2)
    p_all = plot(p1, p2, p3, p4; layout=(2, 2), size=(900, 600),
                 plot_title="equatorial_circular [GTC] — lofted vs monotonic (Fig. 7)")
    savefig(p_all, joinpath(outdir, "equatorial_circular_GTC_fig7.png"))
    @printf("    saved Fig.-7 comparison plot\n")
end

# ─── Monte Carlo of Table 4/5 (equatorial circular scenario) ─────────
"""
Replicates the paper's convergence analysis: n samples with uniform
pitch ∈ [60°, 90°], launch azimuth ∈ [−180°, 180°], t_f guess
∈ [1800 s, 2600 s] (Table 4).  Paper result (Table 5, 500 samples):
GTC 48.2 % convergence (iters min 8 / mean 21 / max 50), CTC 0 %.
"""
function monte_carlo(sc::Scenario, mode::Symbol, n; rng=Random.Xoshiro(1))
    @printf("\n━━━ Monte Carlo  %-22s  [%s]  n=%d ━━━\n", sc.name, mode, n)
    iters = Int[]
    n_conv = 0
    for s in 1:n
        pitch = 60.0 + 30.0 * rand(rng)
        azim  = -180.0 + 360.0 * rand(rng)
        tf    = 1800.0 + 800.0 * rand(rng)
        c = run_case(sc, mode; verbose=false,
                     pitch_deg=pitch, azimuth_deg=azim, t_f_s=tf)
        conv = c.res.ok
        n_conv += conv
        conv && push!(iters, length(c.res.history))
        @printf("  sample %3d/%d  pitch=%5.1f°  az=%+7.1f°  t_f0=%6.1f s  →  %s\n",
                s, n, pitch, azim, tf,
                conv ? @sprintf("converged (%d iters)", length(c.res.history)) :
                       "no convergence")
    end
    @printf("  %s convergence rate: %.1f %%  (paper: %s)\n",
            mode, 100 * n_conv / n, mode == :GTC ? "48.2 %" : "0 %")
    if !isempty(iters)
        @printf("  iterations of converged runs: min %d / mean %.1f / max %d  (paper GTC: 8 / 21 / 50)\n",
                minimum(iters), sum(iters)/length(iters), maximum(iters))
    end
end

# ─── Main ────────────────────────────────────────────────────────────
function main()
    @printf("\nSüslü & Söken (2026) — lunar ascent SCvx replication\n")
    @printf("Solver: %s   Discretization: %s%s   Trust region: %s\n",
            SOLVER_NAME, DISCRETIZATION,
            DISCRETIZATION == :foh ? " (matrix-exponential FOH — deviates from eq. 54)" : " (eq. 54)",
            TR_UPDATE)
    @printf("Canonical units: DU=%.1f km, TU=%.2f s, VU=%.2f m/s, AU=%.4f m/s²\n",
            DU/1000, TU, VU, AU)

    outdir = @__DIR__
    cases = []

    # GTC for all three scenarios
    for sc in SCENARIOS
        push!(cases, run_case(sc, :GTC))
    end
    # CTC for all three scenarios — paper documents failure for 2 & 3
    for sc in SCENARIOS
        c = run_case(sc, :CTC)
        push!(cases, c)
        # Sec. IV.C: on equatorial CTC failure the paper retries with a
        # coarser N = 50 discretization (and fails again).
        if !c.res.ok && sc.name == "equatorial_circular"
            sc50 = Scenario(sc.name * "_N50", sc.lat0_deg, sc.lon0_deg,
                            sc.alt0_m, sc.a_t_km, sc.e_t, sc.i_t_deg,
                            sc.Ω_t_deg, sc.ω_t_deg, sc.t_f_s, 50, sc.t_f_fixed)
            push!(cases, run_case(sc50, :CTC; label=sc50.name))
        end
    end

    # Fig. 7: rerun the equatorial GTC case with the altitude forced to
    # be monotonically increasing; the lofted trajectory should retain
    # more final mass (paper: 673.6 kg lofted vs 645.7 kg monotonic).
    sc_eq = SCENARIOS[2]
    case_mono = run_case(sc_eq, :GTC; monotonic_alt=true,
                         label=sc_eq.name * "_mono")
    push!(cases, case_mono)
    case_lofted = cases[2]
    if case_lofted.res.ok && case_mono.res.ok
        @printf("\n  Fig. 7 comparison: lofted m_f = %.1f kg  vs  monotonic m_f = %.1f kg\n",
                exp(case_lofted.res.x[7, end]), exp(case_mono.res.x[7, end]))
    end

    @printf("\n━━━━━━━━━━━━━━━━━━━━━━━━ SUMMARY ━━━━━━━━━━━━━━━━━━━━━━━━\n")
    @printf("%-26s  %-4s  %-12s  %-10s  %-8s  %-8s  %-6s\n",
            "scenario", "mode", "status", "t_f [s]", "m_f [kg]", "wall [s]", "iters")
    for c in cases
        n_iter = length(c.res.history)
        if c.res.ok
            mf = exp(c.res.x[7, end])
            @printf("%-26s  %-4s  %-12s  %-10.1f  %-8.1f  %-8.2f  %-6d\n",
                c.label, c.mode, "converged", c.res.p*TU, mf, c.elapsed, n_iter)
        else
            @printf("%-26s  %-4s  %-12s  %-10s  %-8s  %-8.2f  %-6d\n",
                c.label, c.mode, "FAILED", "—", "—", c.elapsed, n_iter)
        end
    end

    # Plot only the successful cases
    for c in cases
        plot_trajectory(c; outdir=outdir)
    end
    plot_fig7(case_lofted, case_mono; outdir=outdir)

    @printf("\nPaper benchmarks (for reference):\n")
    @printf("  inclined_circular [GTC]: a=1937.4 km  Ω=280.0°  (CTC: Ω=270.034°)\n")
    @printf("  equatorial_circular [GTC]: converges 48.2%% over 500 MC samples;\n")
    @printf("                              CTC: 0%% convergence\n")
    @printf("  elliptic [GTC]: a=3087.4 km  e=0.37248  i=40°  Ω=290°  ω=90°\n")
    @printf("  Fig. 7: lofted m_f = 673.6 kg vs monotonic m_f = 645.7 kg\n")

    # Optional Monte Carlo (Table 4/5), e.g.  SUSLU_MC=50 julia suslu-soken-2026.jl
    n_mc = parse(Int, get(ENV, "SUSLU_MC", "0"))
    if n_mc > 0
        monte_carlo(SCENARIOS[2], :GTC, n_mc)
        monte_carlo(SCENARIOS[2], :CTC, n_mc)
    end
end

main()