#!/usr/bin/env python3
"""FUNNY-TM-2026-003 — calculation of record (heat rejection: sunlit LEO vs Martian surface).
Every number in the memorandum comes from here. Exact intermediates; rounding only at presentation.
Convection uses the published module funny_tm_2026_003_conv.py, kept next to this script.
Usage: python3 funny_tm_2026_003_calc.py  -> prints results as JSON."""
import json, math, sys, io, contextlib, importlib.util, pathlib
sys.dont_write_bytecode = True

SIGMA = 5.670374419e-8          # W m-2 K-4 [P] CODATA
K0 = 273.15
# ---- assumptions (A-numbers as in the memorandum)
EPS, ALPHA = 0.88, 0.15         # A1 radiator IR emissivity, solar absorptivity [E]
H_ALT, R_E = 550.0, 6371.0      # A2 km, circular LEO, nadir-facing plate [A]
S_EARTH = 1361.0                # A3 W m-2 [M] NSSDC Earth fact sheet
T_E = (254.0, 255.0)            # A3 Earth effective temperature: NSSDC 254.0 K; white paper v001 used 255 K (240 W m-2)
ALB = (0.294, 0.30)             # A3 Earth Bond albedo: NSSDC 0.294; white paper v001 used 0.30
AREAL = (5.0, 8.0, 10.0)        # A4 radiator panel areal mass, kg m-2 (low, central, high) [E]
AREAL_LIGHT = (1.6, 3.9)        # A4 lightweight panels, kg m-2 [E]
MULT = (1.4, 1.55, 1.7)         # A4 system multiplier (loops, pumps, structure, MMOD) [E]
ATT = (15.0, 22.0)              # A5 planning area per 20 kW, m2 (favourable .. unfavourable attitude) [E]
CP_W, DT_W = 4186.0, 20.0       # A6 water, J kg-1 K-1; loop temperature rise, K [P]
HFG_NH3 = {"-33 C": 1370e3, "20 C": 1186e3, "80 C": 874e3}   # A6 J kg-1 [M] NIST saturation table
T_SKY_CLEAR = (150.0, 157.5, 165.0)   # A7 Mars clear-night sky effective temperature, K [E]
T_SKY_DUST = (200.0, 225.0)           # A7 global dust storm [E]
DAY_SHIELD_ABS = 52.0           # A8 absorbed, dusty day sky, sun shielded, W m-2 [E]
SUN_ABS, ALB_ABS = 88.0, 22.0   # A8 direct sun (alpha x 590) and ground-reflected, absorbed W m-2 [C/E]
MARS_PEAK_SOLAR = 590.0         # A8 W m-2 [M] Appelbaum & Flood
T_AIR = 210.0                   # A9 Mars near-surface air, K [M] NSSDC mean 214 K; 210 K central
K_REG, L_REG, DT_REG = 0.039, 0.10, 50.0   # A10 regolith W m-1 K-1 [M] Grott 2021; layer m; K
Q20, Q100 = 20e3, 100e3         # W, heat loads

sb = lambda T: SIGMA * T**4
def emit(tc, eps=1.0): return eps * sb(tc + K0)

def orbit(TE, alb):
    rho = math.asin(R_E / (R_E + H_ALT)); F = math.sin(rho)**2
    ir = EPS * F * sb(TE); albedo = ALPHA * F * alb * S_EARTH
    net = {t: emit(t, EPS) - ir - albedo for t in (20, 80)}
    return dict(rho_deg=math.degrees(rho), F=F, olr=sb(TE), ir=ir, alb=albedo, net20=net[20], net80=net[80],
                plate80=emit(80, EPS) + net[80])

O = orbit(T_E[0], ALB[0])          # central: NSSDC inputs
O_WP = orbit(T_E[1], ALB[1])       # white paper v001 inputs

# ---- edge-on attitude (Rev C): plate plane holds the Sun and nadir directions, both faces radiate
def f_quad(rho, lam, n=300):
    """Numerical check: view factor, one face of a plate element to a sphere of angular radius rho, face normal at lam from nadir."""
    tot, dth, dph = 0.0, (math.pi / 2) / n, (2 * math.pi) / (2 * n)
    nad = (math.sin(lam), 0.0, math.cos(lam))
    for i in range(n):
        th = (i + 0.5) * dth; st, ct = math.sin(th), math.cos(th)
        for j in range(2 * n):
            ph = (j + 0.5) * dph
            if nad[0] * st * math.cos(ph) + nad[2] * ct >= math.cos(rho):
                tot += ct * st * dth * dph
    return tot / math.pi

def edge_on(TE, alb):
    rho = math.asin(R_E / (R_E + H_ALT))
    F = (rho - math.sin(rho) * math.cos(rho)) / math.pi        # (3b) each face of a plate whose plane holds nadir [P]
    ir = EPS * F * sb(TE); albedo = ALPHA * F * alb * S_EARTH   # reflected light: full sunlit value (upper bound) .. 0 over the day-night line
    plate = {t: 2 * (emit(t, EPS) - ir) for t in (20, 80)}      # no reflected light
    plate_alb = {t: plate[t] - 2 * albedo for t in (20, 80)}    # full reflected light
    return dict(F=F, F_quad=f_quad(rho, math.pi / 2), ir=ir, alb=albedo,
                plate80=plate_alb[80], plate80_noalb=plate[80], plate20=plate_alb[20], plate20_noalb=plate[20])

EDGE = edge_on(T_E[0], ALB[0])
NZ_NOALB = {t: emit(t, EPS) - O["ir"] + emit(t, EPS) for t in (20, 80)}   # nadir/zenith plate without reflected light
NZ_20 = O["net20"] + emit(20, EPS)                                       # nadir/zenith plate at 20 C, full reflected light
EDGE["gain80"] = (EDGE["plate80"] / O["plate80"] - 1, EDGE["plate80_noalb"] / NZ_NOALB[80] - 1)   # same reflected-light case
EDGE["gain20"] = (EDGE["plate20"] / NZ_20 - 1, EDGE["plate20_noalb"] / NZ_NOALB[20] - 1)
EDGE["nz_plate80_noalb"], EDGE["nz_plate20"], EDGE["nz_plate20_noalb"] = NZ_NOALB[80], NZ_20, NZ_NOALB[20]
EDGE["A_MW"] = (1e6 / EDGE["plate80_noalb"], 1e6 / EDGE["plate80"])
A20_ideal = Q20 / (2 * emit(80, EPS)); A20 = Q20 / O["plate80"]; A100 = Q100 / O["plate80"]
A_MW = 1e6 / O["plate80"]; A_MW_bad = ATT[1] / 20e3 * 1e6
mass = lambda area, areal, mult: area * areal * mult / 1e3            # t per MW
M_MW = dict(central=mass(A_MW, AREAL[1], MULT[1]), formula_lo=mass(A_MW, 8, MULT[0]), formula_hi=mass(A_MW, 8, MULT[2]),
            band_lo=mass(A_MW, AREAL[0], MULT[0]), band_hi=mass(A_MW_bad, AREAL[2], MULT[2]),
            light_lo=mass(A_MW, AREAL_LIGHT[0], MULT[0]), light_hi=mass(A_MW_bad, AREAL_LIGHT[1], MULT[2]))
water = Q100 / (CP_W * DT_W); nh3 = {k: Q100 / v for k, v in HFG_NH3.items()}
ratio = {k: water / v for k, v in nh3.items()}

# ---- convection: the module next to this script (the build script hands the same module in)
conv = sys.modules.get("conv")
if conv is None:
    _here = pathlib.Path(__file__).resolve().parent
    _conv = next((p for p in (_here / "funny_tm_2026_003_conv.py", _here / "conv.py") if p.exists()), None)
    if _conv is None:
        sys.exit("funny_tm_2026_003_conv.py not found next to this script")
    spec = importlib.util.spec_from_file_location("conv", _conv)
    conv = importlib.util.module_from_spec(spec)
    with contextlib.redirect_stdout(io.StringIO()):
        spec.loader.exec_module(conv)
CONV = {"natural_1m": conv.natural(1.0, 60.0)["h"]}
for U in (0.5, 4.0, 10.0, 25.0):
    CONV[f"1m_{U:g}"] = conv.forced(1.0, U)["h_mix"]
for L in (0.004, 0.02, 0.1):
    CONV[f"{L:g}m_4"] = conv.forced(L, 4.0)["h_mix"]
H_CENTRAL = 0.5   # A11 field-scale central value used for the Mars budget [C]
q_conv = lambda h, Tr: h * ((Tr + K0) - T_AIR)

# ---- Mars radiative cases, face-up radiator at 80 C
E80 = emit(80, EPS)
sky = lambda T: EPS * sb(T)
MARS = {"night_clear": [E80 - sky(T) for T in T_SKY_CLEAR], "night_dust": [E80 - sky(T) for T in T_SKY_DUST],
        "day_shielded": E80 - DAY_SHIELD_ABS, "day_exposed": E80 - DAY_SHIELD_ABS - SUN_ABS - ALB_ABS,
        "abs_clear": [sky(T) for T in T_SKY_CLEAR], "abs_dust": [sky(T) for T in T_SKY_DUST]}
MARS["conv80_central"] = q_conv(H_CENTRAL, 80); MARS["conv80_band"] = (q_conv(0.2, 80), q_conv(1.0, 80))
mars_face = MARS["day_shielded"] + MARS["conv80_central"]
PLATE_RATIO = O["plate80"] / mars_face            # orbital two-sided plate vs Mars face-up (day shielded + convection)
EDGE["M_MW"] = tuple(mass(a, AREAL[1], MULT[1]) for a in EDGE["A_MW"])   # central areal mass and multiplier (A4)
EDGE["ratio_mars"] = (EDGE["plate80"] / mars_face, EDGE["plate80_noalb"] / mars_face,      # vs Mars face-up with convection
                      EDGE["plate80_noalb"] / MARS["day_shielded"])                          # .. and radiation only
CONV_SHARE = MARS["conv80_central"] / mars_face
q_reg = K_REG * DT_REG / L_REG; A_reg_kW = 1e3 / q_reg
BLOWN = dict(q_footprint=930.0, fan_pct=(1.0, 4.0))   # finned HX, per footprint, fan power fraction [C]

R = dict(emit={t: (emit(t), emit(t, EPS)) for t in (20, 80, 97)}, orbit=O, orbit_wp=O_WP,
         A20_ideal=A20_ideal, A20=A20, A20_band=ATT, A100=A100, A100_band=(ATT[0] * 5, ATT[1] * 5),
         A_MW=A_MW, A_MW_bad=A_MW_bad, M_MW=M_MW, water=water, nh3=nh3, ratio=ratio,
         conv=CONV, mars=MARS, mars_face=mars_face, plate_ratio=PLATE_RATIO, conv_share=CONV_SHARE,
         q_reg=q_reg, A_reg_kW=A_reg_kW, blown=BLOWN, edge=EDGE,
         sky_ratio=(MARS["night_clear"][1] / O["net80"], MARS["night_dust"][1] / O["net80"]))
if __name__ == "__main__":
    print(json.dumps(R, indent=1, default=lambda x: round(x, 4)))