#!/usr/bin/env python3
"""
Week 1 thermal budget for the benchtop thermal-vacuum testbed.

Question this answers: at what chamber pressure does the radiative channel
between two payload pucks dominate the parasitic channels (gas conduction,
convection, standoff conduction), and what pump class does that imply?

Geometry: two parallel 50x50 mm PCB pucks, 20 mm gap, hot puck held at
T_hot by its heater trace, cold puck floating near ambient chamber wall
temperature. All parameters at the top; change them and rerun.

Models (all standard, all first-order -- this is a budget, not a paper):
- Radiation: parallel-plate grey exchange, effective emissivity
  1/(1/e1 + 1/e2 - 1).
- Gas conduction: series (continuum + free-molecular) interpolation,
  q = dT / (d/k_air + 1/(Lambda*p)). Valid across the transition regime
  to ~20-30% -- fine for a budget.
- Convection: Rayleigh number scaling Ra ~ p^2; convection assumed off
  below Ra_crit = 1708 (enclosure criterion).
- Standoffs: pure 1D conduction through N plastic posts.

Honesty notes:
- Accommodation coefficient set to 0.9 (typical engineering surfaces);
  free-molecular constant Lambda for air ~ 1.2 W m^-2 K^-1 Pa^-1 with
  that accommodation. If you care about the last 20%, measure it -- the
  pressure sweep does exactly that.
- Radiation to the chamber wall (view factor to surroundings) is a larger
  absolute flux than puck-to-puck exchange; the experiment measures the
  puck-to-puck *differential* via the occluder, so the wall term is common
  mode. Budget below is the puck-to-puck channel vs its parasitics.
"""

import numpy as np
from pathlib import Path
import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt

# ----------------------------- parameters -----------------------------
L = 0.050          # puck side length [m]
A = L * L          # facing area [m^2]
d = 0.020          # puck-to-puck gap [m]
T_hot = 350.0      # hot puck [K]  (~77 C -- safe for FR-4)
T_cold = 300.0     # cold puck / ambient [K]
eps1 = eps2 = 0.90 # soldermask emissivity (matte, typical 0.85-0.95)

k_air = 0.026      # air thermal conductivity, continuum [W/m/K]
Lambda = 1.2       # free-molecular conduction const, air, alpha~0.9 [W/m^2/K/Pa]

n_posts = 3        # standoffs supporting the floating puck
post_d = 0.003     # standoff diameter [m]
post_len = 0.015   # standoff length [m]
k_plastic = 0.20   # PLA/PETG [W/m/K]

sigma = 5.670e-8
dT = T_hot - T_cold

# ----------------------------- channels -------------------------------
p = np.logspace(-2, 5.01, 500)  # 0.01 Pa .. 1 atm

# Radiation (pressure-independent)
eps_eff = 1.0 / (1.0 / eps1 + 1.0 / eps2 - 1.0)
Q_rad = eps_eff * sigma * (T_hot**4 - T_cold**4) * A * np.ones_like(p)

# Gas conduction: series continuum + free-molecular
Q_gas = dT * A / (d / k_air + 1.0 / (Lambda * p))

# Convection threshold: Ra scales as p^2 (rho^2 in nu*alpha denominator)
g, beta = 9.81, 1.0 / ((T_hot + T_cold) / 2)
nu_atm, alpha_atm, p_atm = 1.6e-5, 2.2e-5, 101325.0
Ra_atm = g * beta * dT * d**3 / (nu_atm * alpha_atm)
p_conv_off = p_atm * np.sqrt(1708.0 / Ra_atm)

# Standoff conduction (pressure-independent)
A_post = np.pi * (post_d / 2) ** 2
Q_post = n_posts * k_plastic * A_post / post_len * dT * np.ones_like(p)

# ----------------------------- report ---------------------------------
def q_at(pressure):
    return dT * A / (d / k_air + 1.0 / (Lambda * pressure))

print(f"Puck-to-puck radiative exchange:      {Q_rad[0]*1000:7.1f} mW  (pressure-independent)")
print(f"Standoff conduction ({n_posts}x PLA posts):   {Q_post[0]*1000:7.1f} mW  (pressure-independent)")
print(f"Convection suppressed below:          {p_conv_off/1000:7.1f} kPa  (Ra_atm = {Ra_atm:.0f})")
print()
print(f"{'p [Pa]':>10} {'Q_gas [mW]':>12} {'gas/rad':>9}   pump class")
rows = [
    (2000.0, "diaphragm pump (~2 kPa ultimate)"),
    (100.0, "single-stage rotary vane"),
    (5.0, "two-stage rotary vane, typical working"),
    (1.0, "two-stage rotary vane, good"),
    (0.1, "RV + careful seals / small turbo"),
]
for pr, label in rows:
    qg = q_at(pr)
    print(f"{pr:>10.1f} {qg*1000:>12.1f} {qg/Q_rad[0]:>8.1%}   {label}")

print()
print("Decision rule: gas conduction < 10% of radiative signal wants p <~ 1 Pa;")
print("< 25% (correctable via pressure sweep) is satisfied by ~5 Pa.")
print("=> two-stage rotary vane pump. Diaphragm pump is not sufficient:")
print(f"   at 2 kPa gas conduction is {q_at(2000)/Q_rad[0]:.0%} of the radiative signal.")
print()
print("Method note: Q_rad is pressure-independent and Q_gas(p) is monotone ->")
print("sweeping pressure and extrapolating the measured exchange to p -> 0")
print("separates the two channels empirically, independent of these models.")

# ----------------------------- plot -----------------------------------
fig, ax = plt.subplots(figsize=(8, 5.2))
ax.loglog(p, Q_rad * 1000, lw=2.2, label="Radiation (puck-to-puck)")
ax.loglog(p, Q_gas * 1000, lw=2.2, label="Gas conduction (series model)")
ax.loglog(p, Q_post * 1000, lw=1.6, ls="--", label=f"Standoff conduction ({n_posts}x PLA)")
ax.axvspan(p_conv_off, p[-1], alpha=0.12, color="tab:red")
ax.text(np.sqrt(p_conv_off * p[-1]), 2600, "convection\nactive", ha="center", fontsize=9)
for pr, name in [(2000, "diaphragm"), (1.0, "2-stage RV")]:
    ax.axvline(pr, color="grey", lw=0.8, ls=":")
    ax.text(pr, 0.006, f" {name}", rotation=90, va="bottom", fontsize=8, color="grey")
ax.set_xlabel("Chamber pressure [Pa]")
ax.set_ylabel("Heat flow between pucks [mW]")
ax.set_title(
    f"Thermal channel budget - {L*1e3:.0f} mm pucks, {d*1e3:.0f} mm gap, "
    f"dT = {dT:.0f} K"
)
ax.grid(True, which="both", alpha=0.25)
ax.legend(loc="center left", fontsize=9)
fig.tight_layout()
render_dir = Path(__file__).resolve().parent / "output" / "renders"
render_dir.mkdir(parents=True, exist_ok=True)
fig.savefig(render_dir / "thermal_budget.png", dpi=160)
print("wrote thermal_budget.png")