#!/usr/bin/env python3
"""FUNNY-TM-2026-001 — calculation of record. Every number in the memorandum comes from here.
Exact intermediates; rounding only at presentation.
Usage: python3 funny_tm_2026_001_calc.py  -> prints results as JSON."""
import json

# ---- constants [P]
MW = dict(H2=2.01588, CH4=16.0425, H2O=18.01528, CO2=44.0095, O2=31.9988)  # g/mol
OF = 3.6
H_PER_DAY = 24.0
WINDOW_D = 780.0            # synodic period
YEAR_D = 365.25

# ---- assumptions (A-numbers as in the memorandum)
A1_m_p = 1200.0             # t, band 800–1500
A2_el = 55.5                # kWh/kg H2 (60 % LHV), band 44.4–66.7
A4_tw, A4_fd = 500.0, 0.85  # d, duty
A5_s_fsp = 40.0 / 6000.0    # kWe/kg
A6_s_mw = (0.02, 0.06)      # kWe/kg, single source
A7_Pf, A7_Mf, A7_Mf_band = 850.0, 80.0, (67.0, 119.0)  # MWe net, kt — corrected 2026-09-27 (Rev F: radiator field 7–30 → 12–64 kt)
A8_tiers = (33.0, 18.6)     # kt
A9_life_yr = 10.0
fus_wear_kt_yr = (0.8, 1.0, 1.2)   # blanket/first-wall stream [E]
A11_solar_st = (0.003, 0.010)      # kWe/kg, with storage
solar_nost = (0.005, 0.009)        # kWe-avg/kg, deferrable, no storage [E]

# Fusion plant landed-mass basis (internal report 034, conditional on the Lion et al. 2025 design point) — kt
# (component, central, low, high, basis tag)
FUSION_MASS = [
    ("48 HTS coils (winding packs incl. SS316 radial plates)", 16.2, 16.2, 16.2, "P"),   # Lion Table 8 unique-coil masses x8
    ("Coil casings + inter-coil support structure",           12.0,  8.0, 16.0, "E"),
    ("Breeding blanket (WCLL, 86 % PbLi by mass)",               4.3,  4.3,  4.3, "E"),
    ("First wall + divertor targets",                           0.3,  0.3,  0.3, "E"),
    ("In-vessel neutron shield (WC-based)",                     4.0,  3.0,  4.5, "E"),
    ("Vacuum vessel",                                           2.5,  2.5,  2.5, "E"),
    ("Cryostat + thermal shield + feeders",                     3.5,  3.5,  3.5, "E"),
    ("ECRH heating system",                                     0.8,  0.8,  0.8, "E"),
    ("Cryoplant (20 K)",                                        0.8,  0.8,  0.8, "E"),
    ("Balance of plant (turbines, steam generators, pumps, electrical, tritium plant)", 10.0, 10.0, 10.0, "E"),
    ("Radiator field, Mars, radiator at 40–80 °C (corrected in Rev F)", 25.0, 12.0, 64.0, "C"),
]
# Derating chain, Mars (internal reports 030, 034) — MW
DERATE = dict(fusion_MW=2700, thermal_MW=3150, gross_earth_MWe=1054, gross_mars_MWe=(871, 1018), house_MW=90, net_mars_MWe=(780, 930), net_central_MWe=850)

def feed(m_p):
    ch4 = m_p / (1 + OF); lox = m_p - ch4
    n = ch4 / MW["CH4"]                       # Mmol-scale, units cancel (t / g/mol)
    h2 = 4 * n * MW["H2"]
    w_el = (h2 / MW["H2"]) * MW["H2O"]
    w_ret = 2 * n * MW["H2O"]
    co2 = n * MW["CO2"]
    o2 = (w_el / MW["H2O"]) / 2 * MW["O2"]
    return dict(ch4=ch4, lox=lox, h2=h2, w_el=w_el, w_ret=w_ret, w_net=w_el - w_ret, co2=co2, o2=o2, o2_surplus=o2 - lox)

def energy(m_p=A1_m_p, el=A2_el, co2_k=0.4, water_k=0.27 + 0.05, liq=(0.6, 0.9),
           zbo_1200=0.936 * 0.64, overhead=0.07, sab_credit_1200=0.0):
    """GWh per ship. water_k = 0.27 kWh_th/kg (= 2 x 0.134 melt minimum, supplied electrically) + 0.05 kWh_e/kg."""
    f = feed(m_p); s = m_p / A1_m_p
    rows = dict(
        electrolysis=f["h2"] * el / 1e3,
        co2=f["co2"] * co2_k / 1e3,
        water=f["w_net"] * water_k / 1e3,
        liquefaction=(f["lox"] * liq[0] + f["ch4"] * liq[1]) / 1e3,
        zbo=zbo_1200 * s ** (2 / 3),
        sabatier=-sab_credit_1200 * s,
    )
    sub = sum(rows.values()); rows["overhead"] = sub * overhead
    rows["total"] = sub + rows["overhead"]
    rows["e_p_kWh_kg"] = rows["total"] * 1e6 / (m_p * 1e3)
    return rows

def p_cap(E, tw=A4_tw, fd=A4_fd): return E * 1e3 / (tw * H_PER_DAY * fd)       # MWe capacity
def p_avg(E, tw=A4_tw): return E * 1e3 / (tw * H_PER_DAY)                      # MW time average
def m_per_ship(P, s): return P * 1e3 / s / 1e3                                   # t

E_c = energy()
E_lo = energy(m_p=800, el=44.4, co2_k=0.2, water_k=0.134 + 0.05, liq=(0.4, 0.6), zbo_1200=0.3, overhead=0.05, sab_credit_1200=0.33)
E_hi = energy(m_p=1500, el=66.7, co2_k=0.8, water_k=4 * 0.134 + 0.05, liq=(1.2, 1.6), zbo_1200=1.2, overhead=0.10)
P_c = p_cap(E_c["total"])
r = WINDOW_D / (A9_life_yr * YEAR_D)            # fraction of fission fleet re-landed per window
wear_w = [w * WINDOW_D / YEAR_D for w in fus_wear_kt_yr]   # kt per window

def nstar(Mf_kt, s, P=P_c): return Mf_kt * 1e3 / m_per_ship(P, s)
def nstar5(s, W=wear_w[1], P=P_c, w=5):
    """Eq. (5): M_fus(w) = M_f0 + W(w-1);  M_fis(w) = N m_fis (1 + r(w-1));  N* = M_fus(w) / (m_fis (1 + r(w-1)))"""
    return (A7_Mf + W * (w - 1)) * 1e3 / (m_per_ship(P, s) * (1 + r * (w - 1)))

# Moon
lox_leg = (130 * OF / (1 + OF), 290 * OF / (1 + OF))    # t LOX per leg (NRHO light ... direct-Earth heavy)
ch4_leg = (130 / (1 + OF), 290 / (1 + OF))
e_lox = 10.8                                             # kWh/kg all-electric
def moon_P(legs, lox_t): return legs * lox_t * e_lox / 1e3 * 1e3 / (YEAR_D * H_PER_DAY) / A4_fd   # MWe cap

R = dict(
    fusion_mass=FUSION_MASS, derate=DERATE,
    fusion_mass_sum=(sum(r[1] for r in FUSION_MASS), sum(r[2] for r in FUSION_MASS), sum(r[3] for r in FUSION_MASS)),
    feed=feed(A1_m_p), E_c=E_c, E_lo=E_lo["total"], E_hi=E_hi["total"],
    P_cap=P_c, P_avg=p_avg(E_c["total"]),
    table4={f"{E:.1f}": {str(int(t)): p_cap(E, t) for t in (300, 500, 600, 780)} for E in (E_lo["total"], E_c["total"], E_hi["total"])},
    fis_per_ship={"fsp": m_per_ship(P_c, A5_s_fsp), "mw02": m_per_ship(P_c, A6_s_mw[0]), "mw06": m_per_ship(P_c, A6_s_mw[1])},
    ships_per_plant={f"{E:.1f}": {str(int(t)): A7_Pf / p_cap(E, t) for t in (300, 500, 600, 780)} for E in (E_c["total"], E_hi["total"], E_lo["total"])},
    nstar={b: {k: nstar(M, s) for k, s in (("fsp", A5_s_fsp), ("mw02", A6_s_mw[0]), ("mw06", A6_s_mw[1]))}
           for b, M in (("land", A7_Mf), ("tier1", A8_tiers[0]), ("tier2", A8_tiers[1]))},
    nstar_land_range=(A7_Mf_band[0] * 1e3 / m_per_ship(p_cap(E_c["total"], 300), A5_s_fsp),
                      A7_Mf_band[1] * 1e3 / m_per_ship(p_cap(E_c["total"], 600), A5_s_fsp)),
    r=r, wear_w=wear_w,
    nstar5={"fsp": [nstar5(A5_s_fsp, W) for W in wear_w], "mw02": [nstar5(A6_s_mw[0], W) for W in wear_w]},
    solar={"st_t": [m_per_ship(P_c, s) for s in A11_solar_st[::-1]], "nost_t": [m_per_ship(P_c, s) for s in solar_nost[::-1]]},
    moon={"lox": lox_leg, "ch4": ch4_leg, "E_leg": [x * e_lox / 1e3 for x in lox_leg],
          "P": {n: [moon_P(n, x) for x in lox_leg] for n in (10, 100, 1000)}},
    dust_loss=(120 / A4_tw) * (1 - 0.35),
    sens={"m800": energy(m_p=800)["total"], "m1500": energy(m_p=1500)["total"],
          "el50": energy(el=66.7)["total"], "el75": energy(el=44.4)["total"],
          "sab_full": energy(sab_credit_1200=0.33)["total"]},
)
R["solar"]["nstar_st"] = [A7_Mf * 1e3 / t for t in R["solar"]["st_t"]][::-1]
R["solar"]["nstar_nost"] = [A7_Mf * 1e3 / t for t in R["solar"]["nost_t"]][::-1]
R["sens"]["nstar_m"] = [nstar(A7_Mf, A5_s_fsp, p_cap(R["sens"]["m800"])), nstar(A7_Mf, A5_s_fsp, p_cap(R["sens"]["m1500"]))]
if __name__ == "__main__":
    print(json.dumps(R, indent=1, default=lambda x: round(x, 4)))