import anvil
from anvil import Q, K, Pa, m, s, kg, N, J, W, kPa, MPa

# ── Quantities with units ──────────────────────────────────────
T   = 300 * K           # Q(300, "K")
P   = 6.9 * MPa         # Q(6900000, "Pa")
rho = 1.225 * kg/m**3   # Q(1.225, "kg/m^3")

# ── Unit arithmetic ───────────────────────────────────────────
F   = Q(100, "N")
A   = Q(0.01, "m^2")
sig = F / A             # → Q(10000, "Pa"), dim auto-detected
KE  = 0.5 * Q(10,"kg") * Q(30,"m/s")**2  # → Q(4500, "J")

# ── Unit conversions ──────────────────────────────────────────
P.to("psi")             # Q(1000.7, "psi")
Q(300, "K").to("R")     # Q(540, "R")

# ── Call a built-in RSQ directly ─────────────────────────────
r = anvil.R.isentropic_ratios(M=2.0, gamma=1.4)
# → {"T0_T": 1.8, "P0_P": 7.824, "rho0_rho": 4.347}

# ── Load a pre-built System ───────────────────────────────────
nozzle = anvil.S.rocket_nozzle.copy()
nozzle.set(P0=10e6, T0=3500)
result = nozzle.solve_forward()
result.summary()

# ── Build your own System ─────────────────────────────────────
sys = anvil.system("rayleigh_duct")
sys.add("M1",    0.3)
sys.add("T01",   400.0, "K")
sys.add("P1",    200e3, "Pa")
sys.add("q_heat", 300e3, "J/kg")
sys.add("cp",    1005.0, "J/kg/K")
sys.use(my_relation)
result = sys.solve_forward()

# ── Parametric sweep ──────────────────────────────────────────
import numpy as np
sweep = sys.sweep("q_heat", np.linspace(0, 500e3, 50))
sweep.summary(outputs=["M2", "T02", "P02_P01"])

# ── Register your own RSQ ─────────────────────────────────────
proj = anvil.project("my_study", path="./work")
proj.push(my_func, domain="aero", tags=["compressible"])
proj.R.my_func(M=2.0)
proj.promote("my_func")   # → global registry

src/anvil/
├── __init__.py          top-level API: system(), relation(), push(), solve(), R, S, QDB
├── quantity.py          Quantity (Q), Dim arithmetic
├── units.py             UnitDB, UnitStub, all unit definitions
├── relation.py          Relation class, Relation.block()
├── system.py            System, Result, SweepResult, SensitivityResult, OptimizeResult
├── solvers/
│   └── __init__.py      find_root, solve_nonlinear, solve_ode, solve_ode_stiff,
│                        solve_bvp, solve_pde_heat_1d, minimize, minimize_global
├── decomp.py            POD, DMD, hankel, pod_reconstruct, dmd_reconstruct, pod_rank
├── doe.py               full_factorial, latin_hypercube, sobol, run_doe
├── registry/
│   ├── __init__.py      push, search, list, info, export, remove
│   ├── store.py         SQLite backend (Store class)
│   ├── namespace.py     R., S., QDB. dot-access namespaces
│   └── loader.py        RSQ source → live object
├── seed.py              101 built-in RSQs seeded on first import
├── project.py           Project class (isolated registry)
├── adapter.py           Adapter class (python + cli backends)
├── adapters/            real-only wrappers: XFOIL, SU2, OpenFOAM, FEniCSx,
│                        pyNastran, OpenMDAO, Cantera, NASA CEA, CoolProp,
│                        gmsh, RocketCEA, poliastro, pykep, surrogates, UQ
├── viz.py               convergence, sweep_plot, variable_trace, dependency_graph,
│                        pod_energy, dmd_spectrum
├── inspect.py           anvil.check()
├── watchdog.py          Watchdog convergence tracker
├── db/
│   ├── __init__.py      const, fluids, materials
│   └── properties.py    FluidDB, MaterialDB data tables
├── help_.py             anvil.lookup(), in-REPL help
└── cfd/                 CFD solver (mesh, BCs, flux, viz)

anvil_server/            FastAPI backend for the Web Workbench (REST + WebSocket)
anvil_web/               React + TypeScript frontend (calculator + canvas)

git clone https://github.com/c-rk/anvil-framework-s
cd anvil-framework-s
python start_anvil.py

pip install -e .                 # core (numpy + scipy)
pip install -e ".[server]"       # + web workbench (FastAPI/uvicorn)
pip install -e ".[viz]"          # + matplotlib plots
pip install -e ".[poliastro]"    # poliastro + astropy
pip install -e ".[pykep]"        # pykep
pip install -e ".[adapters]"     # all adapter deps
pip install -e ".[all]"          # viz + all adapters + server
conda install -c cantera cantera # Cantera (conda recommended)

anvil doctor     # report which adapters are usable on this machine, with install hints
anvil serve      # start the web workbench (same as start_anvil.py, without the venv setup)
anvil version    # print the installed version

from anvil.client import AnvilClient

srv = AnvilClient("http://127.0.0.1:8000")
srv.health()                                    # {'status': 'ok', 'rsq_count': ...}
srv.call("isentropic_ratios", M=2.0, gamma=1.4) # -> {'T0_T': 1.8, ...}

import anvil
from anvil import Q, K, Pa, m, s, kg, N, J, W, kPa, MPa

# Three equivalent styles
T = Q(300, "K")          # string unit
T = Q(300, "temperature") # category alias
T = 300 * K              # unit-stub syntax (cleanest)

P   = 101325 * Pa
v   = 340 * (m/s)
g   = 9.81 * m/s**2
rho = 1.225 * kg/m**3

F   = Q(100, "N")
A   = Q(0.01, "m^2")
sig = F / A              # → 10000.00 Pa  (dimension auto-detected)

KE  = 0.5 * Q(10,"kg") * Q(30,"m/s")**2  # → 4500.00 J

print(sig)               # 10000.00 Pa
print(KE)                # 4500.00 J
print(KE.to("kJ"))       # 4.5000 kJ
print(KE.si)             # 4500.0  (always SI float)

# Direct call, no System needed
r = anvil.R.isentropic_ratios(M=2.0, gamma=1.4)
print(r["T0_T"])    # 1.8
print(r["P0_P"])    # 7.8244
print(r["rho0_rho"]) # 4.3469

nozzle = anvil.S.rocket_nozzle.copy()   # deep copy, safe to modify
nozzle.set(P0=10e6, T0=3500)            # override values, keep units

result = nozzle.solve_forward()
result.summary()

--------------------------------------------------------
  rocket_nozzle -- results
--------------------------------------------------------
  P0                        1.0000e+07 Pa
  T0                        3500.00 K
  gamma                     1.2500
  R_gas                     320.00 J/kg/K
  A_throat                  0.01 m^2
  A_exit                    0.08 m^2
  P_amb                     101325.00 Pa
                            ---
  area_ratio                8.0000
  M_exit                    3.6423
  T0_T                      2.6469
  P0_P                      59.6093
  rho0_rho                  22.5243
  T_exit                    1322.66 K
  P_exit                    167773.97 Pa
  a_exit                    728.41 m/s
  V_exit                    2652.17 m/s
  mdot                      12.12 kg/s
  thrust                    32143.09 N
  Isp                       270.20 s
--------------------------------------------------------

def lift_coeff(W, S_ref, V, rho):
    q = 0.5 * rho * V**2
    CL = W / (q * S_ref)
    return {"q_inf": Q(q, "Pa"), "CL": CL}

sys = anvil.system("wing_loading")
sys.add("W",    50000, "N")
sys.add("S_ref", 20,   "m^2")
sys.add("V",    80,   "m/s")
sys.add("rho",  1.225, "kg/m^3")
sys.use(lift_coeff)

result = sys.solve_forward()
result.summary()
# CL = 0.7994, q_inf = 3920.00 Pa

import numpy as np

sweep = sys.sweep("V", np.linspace(50, 150, 8))
sweep.summary(outputs=["CL", "q_inf"])

----------------------------------------------------------------------
  wing_loading -- sweep over V
----------------------------------------------------------------------
               V            CL         q_inf
            [m/s]                      [Pa]
  ----------------------------------------------------------------------
              50        2.0388     1531.25
           64.29        1.2337     2534.06
           78.57        0.8249     3783.89
           92.86        0.5913     5280.73
             107        0.4453     7024.58
           121.4        0.3460     9015.45
           135.7        0.2773     11253.3
             150        0.2269     13781.3
----------------------------------------------------------------------

result = anvil.solve("isentropic_ratios", M=2.0, gamma=1.4)
print(result["P0_P"])   # 7.824

@anvil.relation(domain="aero", tags=["compressible"])
def dynamic_pressure(rho, V):
    return {"q_inf": Q(0.5 * rho * V**2, "Pa")}

# Now available globally
anvil.R.dynamic_pressure(rho=1.225, V=100)
# {"q_inf": 6125.00 Pa}

# System with a feedback loop, use solve_gauss_seidel()
result = sys.solve_gauss_seidel(
    relaxation=0.8,
    max_iter=200,
    rtol=1e-8,
    monitor=True,      # prints live residuals
)

proj = anvil.project("my_study", path="./work")

proj.push(rayleigh_ratios, domain="aero.compressible")
proj.push(rayleigh_heat,  domain="aero.compressible")

result = proj.R.rayleigh_heat(M1=0.3, T01=300.0, P1=101325.0,
                               q_heat=200e3, cp=1005.0)
print(result["M2"])     # 0.428
print(result["T02"])    # 499.00 K

proj.promote("rayleigh_heat")  # → global registry

from anvil import Q

p = Q(101325, "Pa")
T = Q(300, "K")
v = Q(340, "m/s")
mdot = Q(1.5, "g/s")       # compound unit, parsed automatically
vol  = Q(5,  "cm^3")      # ** also accepted: "cm**3"
flux = Q(200, "W/m^2")

p = Q(101325, "pressure")   # → Pa internally
T = Q(300, "temperature")   # → K internally
v = Q(10, "velocity")       # → m/s internally

# Standard (separate brackets)
d1 = Q(1, "[L][M][T-2]")   # → 1 N
# Compact (all in one bracket)
d2 = Q(1, "[LMT-2]")       # → same
# Both produce identical Dim objects
assert d1._dim == d2._dim   # True

from anvil import K, Pa, m, s, kg, mol, N, J, W
from anvil import km, cm, mm, g, kPa, MPa, GPa, bar, atm, psi
from anvil import kN, kJ, MJ, kW, BTU, ft, inch, in_, lb, lbf, hr

T   = 300 * K           # Q(300, "K")
P   = 101325 * Pa       # Q(101325, "Pa")
v   = 340 * (m/s)       # Q(340, "m/s")
g   = 9.81 * m/s**2     # Q(9.81, "m/s^2"), parens optional
rho = 1.225 * kg/m**3   # Q(1.225, "kg/m^3")

Pa*s    # → unit for dynamic viscosity
J/kg/K  # → specific heat unit stub
m**3    # → volume unit stub

import numpy as np
from anvil import Q

pressures = Q(np.array([100e3, 150e3, 200e3]), "Pa")
print(pressures)         # array(3,) Pa
print(pressures.si)      # array([100000. 150000. 200000.])

# Linspace helper
qs = Q.linspace(100, 300, 5, unit="K")  # list of 5 Quantity objects

q = Q(None, "Pa")
print(q.value)  # None
print(q.si)     # None

T = 300 * K
print(T.si)            # 300.0
print(T.value)         # 300.0 (K is SI for temperature)
print(T.unit)          # "K"
print(T.dim)           # [Θ]

p_psi = Q(14.696, "psi")
print(p_psi.si)        # 101325.0 (converted to Pa internally)
print(p_psi.value)     # 14.696  (display in psi)
print(p_psi.unit)      # "psi"

anvil.set_units("Imperial")
Q(101325, "Pa").unit   # "psi", preferred Imperial unit for pressure
anvil.set_units("SI")
Q(101325, "Pa").unit   # "Pa"

Q(100, "Pa") + Q(50, "Pa")      # → 150.00 Pa
Q(1000, "N") - Q(200, "N")      # → 800.00 N

# Dimensionless Q + scalar
Q(1.4) + 0.1                    # → 1.5 (dimensionless)

# Cross-dimension → ValueError
Q(10, "N") + Q(5, "K")          # ValueError: incompatible dimensions

# Dimensional Q + scalar → ValueError
Q(10, "N") + 5                  # ValueError: Cannot add scalar to dimensional quantity

Q(100, "N") * Q(0.1, "m")       # → Q(10, "N*m") → displayed as J (energy)
Q(100, "N") / Q(0.01, "m^2")    # → Q(10000, "Pa")
Q(10, "kg") * Q(9.81, "m/s^2")  # → Q(98.1, "N")

# Scale
Q(100, "Pa") * 2                 # → Q(200, "Pa")
Q(100, "Pa") / 4                 # → Q(25, "Pa")

Q(3, "m") ** 2                   # → Q(9, "m^2")
Q(9, "m^2") ** 0.5               # → Q(3, "m")
Q(3, "m") ** Q(2)                # → Q(9, "m^2"), Q must be dimensionless

Q(3, "m") ** Q(2, "N")          # ValueError: exponent must be dimensionless

-Q(100, "N")          # → Q(-100, "N")
abs(Q(-100, "N"))     # → Q(100, "N")

Q(100, "Pa") == Q(100, "Pa")    # True
Q(100, "Pa") < Q(200, "Pa")     # True
Q(100, "Pa") <= Q(200, "Pa")    # True  (__le__ supported)
Q(200, "Pa") >= Q(100, "Pa")    # True  (__ge__ supported)
Q(100, "Pa") == 100             # True (dimensionless check)

Q(10, "N") < Q(5, "K")
# TypeError: Cannot compare [L][M][T-2] < [Θ]: incompatible dimensions.
#            Convert to the same unit first.

q.to("unit_string")   # returns new Quantity; original unchanged

Q(101325, "Pa").to("atm")      # 1.0000 atm
Q(101325, "Pa").to("psi")      # 14.6959 psi
Q(101325, "Pa").to("kPa")      # 101.3250 kPa
Q(300,   "K").to("R")         # 540.00 R
Q(1000,  "N").to("lbf")       # 224.81 lbf
Q(10,    "m/s").to("cm/s")    # 1000.00 cm/s
Q(10,    "m/s").to("mph")     # 22.37 mph
Q(0.001, "kg/s").to("g/s")    # 1.0000 g/s
Q(1e-6,  "m^3").to("cm^3")    # 1.0000 cm^3
Q(100,   "J").to("BTU")       # 0.0948 BTU

Q(100, "N").to("K")
# ValueError: Cannot convert [L][M][T-2] to 'K' ([Θ]): incompatible dimensions.

T = Q(350, "K", bounds=(200, 500))
T.in_bounds()       # True

T_hot = Q(600, "K", bounds=(200, 500))
T_hot.in_bounds()   # False

# System validation checks bounds automatically
sys.add("T", 350, "K", bounds=(200, 500))
sys.validate()      # warns if any quantity is outside bounds

T = 300 * K
print(f"Temperature: {T:.2f}")   # Temperature: 300.00 (formats display value)
print(f"Temperature: {T:.4e}")   # Temperature: 3.0000e+02

Q(340, "m/s")     # → styled HTML: 340.0000 m/s

temps = Q.linspace(200, 500, 4, unit="K")
# [200.00 K, 300.00 K, 400.00 K, 500.00 K]

ps = Q.array([100e3, 200e3, 300e3], unit="Pa")
# array(3,) Pa
print(ps.si)   # [100000. 200000. 300000.]

q1 = Q(100, "Pa")
q2 = Q(q1, "Pa")     # works: q2.si == 100.0

Q(10, "N") + 5     # ValueError, non-dimensionless Q can't add to scalar
Q(1.4) + 0.1       # OK, dimensionless

q = Q(5, "flarps")
r = q * Q(3, "s")
print(r)            # 15.0000 [T][flarps], custom dim preserved

Q(10, "N") < Q(5, "K")
# TypeError: Cannot compare [L][M][T-2] < [Θ]: incompatible dimensions.

Q(10, "N") >= Q(5, "K")
# TypeError: Cannot compare [L][M][T-2] >= [Θ]: incompatible dimensions.

Q(3, "m") ** Q(2, "N")    # ValueError: Exponent must be dimensionless.

from anvil.units import Dim

force     = Dim(L=1, M=1, T=-2)     # Newton: [L][M][T-2]
pressure  = Dim(L=-1, M=1, T=-2)    # Pascal: [L-1][M][T-2]
dimless   = Dim()                    # dimensionless (empty)
dimless   = Dim.dimensionless()      # singleton version

# Parse from string
d = Dim.parse("[L][M][T-2]")        # standard
d = Dim.parse("[LMT-2]")            # compact, same result
d = Dim.parse("[L2][T-2]")          # with exponents

velocity  = Dim(L=1) / Dim(T=1)         # [L][T-1]
kinetic_e = Dim(M=1) * Dim(L=1, T=-1)**2 # [L2][M][T-2] = J

force = Dim(L=1, M=1, T=-2)
area  = Dim(L=2)
pressure = force / area                 # [L-1][M][T-2] = Pa

Dim(L=1, M=1, T=-2) == Dim(L=1, M=1, T=-2)   # True
Dim(L=1) == Dim(T=1)                            # False
hash(d)                                          # hashable, usable as dict key

from anvil.units import db

scale, dim = db.lookup("Pa")
# scale = 1.0, dim = Dim(L=-1, M=1, T=-2)

scale, dim = db.lookup("psi")
# scale = 6894.757, dim = Dim(L=-1, M=1, T=-2)

scale, dim = db.lookup("cm/s")
# scale = 0.01, dim = Dim(L=1, T=-1), parsed as compound

scale, dim = db.lookup("flarps")
# scale = 1.0, dim = Dim(flarps=1)  , custom dimension created

from anvil.units import Dim
db.find_unit(Dim(L=1, M=1, T=-2), "SI")        # ("N", 1.0)
db.find_unit(Dim(L=1, M=1, T=-2), "Imperial")   # ("lbf", 4.448...)
db.find_unit(Dim(L=-1, M=1, T=-2), "SI")        # ("Pa", 1.0)

db.conversion_factor("Pa", "psi")   # 1/6894.757 ≈ 0.000145
db.conversion_factor("m/s", "ft/s") # 1/0.3048  ≈ 3.28084

db.get_offset("degC")   # 273.15
db.get_offset("°C")     # 273.15
db.get_offset("degF")   # 255.3722...  (= 459.67 × 5/9)
db.get_offset("K")      # 0.0
db.get_offset("Pa")     # 0.0

db.compatible("Pa", "psi")  # True, same dimension
db.compatible("Pa", "K")    # False
db.compatible("degC", "K")  # True, same [Θ] dimension
db.compatible("degF", "R")  # True

from anvil import Q

Q(1.5, "g/s")       # → 0.0015 kg/s     [M][T-1]
Q(10, "cm/s")      # → 0.1 m/s         [L][T-1]
Q(5,  "cm^3")      # → 5e-6 m^3        [L3]
Q(5,  "cm**3")     # ** also accepted
Q(200, "W/m^2")     # → W/m^2           [M][T-3]
Q(1e-3,"Pa*s")      # → Pa·s            [L-1][M][T-1]
Q(1,  "kg*m/s^2")  # → 1 N             [L][M][T-2]
Q(1,  "J/kg/K")    # → J/kg/K          [L2][T-2][Θ-1]
Q(1,  "m^0.5")     # fractional exponent OK

from anvil import Q

Q(25, "degC").si          # 298.15  (stored as K)
Q(25, "degC").value       # 25.0    (display in °C)
Q(25, "degC").to("K")     # 298.15 K
Q(25, "degC").to("degF")  # 77.00 degF
Q(100, "degC").to("degF") # 212.00 degF
Q(32, "degF").to("K")    # 273.15 K
Q(373.15, "K").to("degC") # 100.00 degC

# Unicode forms
Q(0, "°C").si             # 273.15
Q(32, "°F").to("°C")      # 0.00 °C

from anvil import m, s, kg, Pa

v_stub   = m/s          # UnitStub("m/s")
g_stub   = m/s**2       # UnitStub("m/s^2")
rho_stub = kg/m**3      # UnitStub("kg/m^3")

v   = 340 * v_stub      # Q(340, "m/s")
g   = 9.81 * g_stub     # Q(9.81, "m/s^2")
rho = 1.225 * rho_stub  # Q(1.225, "kg/m^3")

K, Pa, m, s, kg, mol, A, N, J, W, rad, deg
km, cm, mm, um, g, tonne, ms, us, hr
kPa, MPa, GPa, bar, atm, psi
kN, kJ, MJ, kW, BTU
ft, inch, in_, lb, lbf
kmol, g_mol, kg_mol

anvil.set_units("SI")        # Pa, N, m/s, K, kg, J, W  (default)
anvil.set_units("Imperial")  # psi, lbf, ft/s, R, slug, BTU, hp

from anvil import Q
q1 = Q(5, "flarps")     # custom dim: Dim(flarps=1)
q2 = Q(3, "widgets")    # custom dim: Dim(widgets=1)

r = q1 * q2             # Dim(flarps=1, widgets=1) → displayed as "[flarps][widgets]"
r2 = q1 / Q(2, "s")    # Dim(T=-1, flarps=1) → "[T-1][flarps]"

# Both are equivalent:
anvil.R.j2_precession(a=6878e3, e=0.001, i_deg=97.4)
anvil.R.j2_precession(a=6878e3, e=0.001, i_deg=Q(97.4, "deg"))

# Also works from a System:
sys.add("i_deg", 97.4, "deg")   # stored as Q(97.4, "deg")
sys.use("j2_precession")         # RSQ receives the Q and converts correctly

from anvil import units

units.resolve("Pa")                  # (1.0, Dim(L=-1, M=1, T=-2))
units.find_unit(Dim(L=1, M=1, T=-2)) # ("N", 1.0)
units.compatible("Pa", "psi")        # True
units.conversion_factor("m", "ft")   # 3.28084
units.list_units()                   # sorted list of all registered unit strings
units.list_categories()              # sorted list of category names
units.set_system("Imperial")
units.get_system()                   # "Imperial"

from anvil import Q

def nozzle_thrust(mdot, V_exit, P_exit, P_amb, A_exit):
    F = mdot * V_exit + (P_exit - P_amb) * A_exit
    return {"thrust": Q(F, "N")}

def isentropic(M, gamma=1.4):
    T_ratio = 1 + ((gamma - 1) / 2) * M**2
    P_ratio = T_ratio ** (gamma / (gamma - 1))
    return {"T0_T": T_ratio, "P0_P": P_ratio}

import anvil

def lift(rho, V, S_ref, CL):
    return {"lift": Q(0.5 * rho * V**2 * S_ref * CL, "N")}

sys = anvil.system("wing")
sys.add("rho", 1.225, "kg/m^3")
sys.add("V", 80, "m/s")
sys.add("S_ref", 20, "m^2")
sys.add("CL", 0.5)
sys.use(lift)

from anvil.relation import Relation

thrust_rel = Relation(
    nozzle_thrust,
    name="nozzle_thrust",
    tags=["propulsion", "nozzle"],
    desc="Rocket thrust from momentum and pressure.",
)

print(thrust_rel._inputs)   # ['mdot', 'V_exit', 'P_exit', 'P_amb', 'A_exit']
print(thrust_rel._outputs)  # [], discovered on first call
thrust_rel(mdot=10, V_exit=3000, P_exit=50000, P_amb=0, A_exit=0.5)
print(thrust_rel._outputs)  # ['thrust']

@anvil.relation(domain="propulsion", tags=["nozzle"])
def rocket_thrust(mdot, V_exit, P_exit, P_amb, A_exit):
    """Rocket thrust from momentum and pressure terms."""
    F = mdot * V_exit + (P_exit - P_amb) * A_exit
    return {"thrust": Q(F, "N")}

# Now callable via anvil.R.rocket_thrust and usable in any system
r = anvil.R.rocket_thrust(mdot=10, V_exit=3000, P_exit=50000, P_amb=0, A_exit=0.5)
print(r["thrust"])  # 30000.00 N

@anvil.relation
def speed_of_sound(gamma, R_gas, T):
    return {"a": Q((gamma * R_gas * T) ** 0.5, "m/s")}

@anvil.relation(domain="test", register=False)
def draft_relation(UA, C_min):
    return {"NTU": UA / C_min}

from anvil.relation import Relation
from anvil import Q

def area_ratio(A_exit, A_throat):
    return {"area_ratio": A_exit / A_throat}

def mach_from_area(area_ratio, gamma=1.4):
    from anvil import solvers
    def residual(M):
        t = (2/(gamma+1))*(1+(gamma-1)/2*M**2)
        return (1/M)*t**((gamma+1)/(2*(gamma-1))) - area_ratio
    M = solvers.find_root(residual, bracket=(1.001, 30.0))
    return {"M_exit": M}

def exit_velocity(M_exit, T0, gamma=1.4, R_gas=287.0):
    T0_T = 1 + (gamma-1)/2 * M_exit**2
    T_e = T0 / T0_T
    a_e = (gamma * R_gas * T_e) ** 0.5
    V = M_exit * a_e
    return {"V_exit": Q(V, "m/s"), "T_exit": Q(T_e, "K")}

nozzle_chain = Relation.block("nozzle_physics",
    steps=[area_ratio, mach_from_area, exit_velocity],
    tags=["propulsion", "nozzle"],
)

print(nozzle_chain._inputs)   # ['A_exit', 'A_throat', 'T0', 'R_gas', 'gamma']
print(nozzle_chain._outputs)  # ['M_exit', 'T_exit', 'V_exit', 'area_ratio']

r = nozzle_chain(A_exit=0.08, A_throat=0.01, T0=3500, gamma=1.25, R_gas=320)
print(r["M_exit"])    # ~3.3
print(r["V_exit"])    # ~2600 m/s

import inspect
sig = inspect.signature(func)
inputs = list(sig.parameters.keys())
defaults = {k: p.default for k, p in sig.parameters.items()
            if p.default is not inspect.Parameter.empty}

import ast, inspect
source = inspect.getsource(func)
tree = ast.parse(source)
# Walk return statements, collect string keys from dict literals

rel = Relation(func)

rel.name        # str, function name or explicitly set
rel.func        # the callable
rel._inputs     # list of input parameter names
rel._outputs    # list of output dict keys (populated lazily)
rel._defaults   # dict of default values
rel.tags        # list of search tags
rel.desc        # description string
rel.info()      # formatted print of inputs/outputs/tags

r = anvil.R.isentropic_ratios(M=2.0, gamma=1.4)
# Returns the raw dict, same as calling the function

r["T0_T"]   # 1.8
r["P0_P"]   # 7.8244

r = anvil.R.isentropic_ratios(M=2.0)  # gamma=1.4 from default

def good(F, A):
    return {"sigma": F / A}    # Q/Q → Quantity with Pa dim

def bad_numpy(x):
    import numpy as np
    return {"y": np.sqrt(x)}  # np.sqrt(Q) fails → falls back to float

# isentropic_ratios expects 'M' but workspace has 'M_exit'
sys.use("isentropic_ratios", map={"M": "M_exit"})

# Map is {relation_param: workspace_name}

def bad(x):
    return x * 2   # not a dict

sys.use(bad)
sys.solve_forward()
# RuntimeError: Relation 'bad' returned 'Quantity' instead of a dict.
#   Relations must return a dict mapping output names to values.
#   Example: return {"result": Q(F, "N")}

def bad_keys(x):
    return {Q(1, "Pa"): x * 2}   # Q as key

# Q is not hashable → TypeError: cannot use Quantity as dict key
# System solve crashes with TypeError

def bad_positional(x, y, /):   # positional-only (Python 3.8+)
    return {"z": x + y}

# Relation wraps it, but calling with **kwargs → TypeError
# inspect.signature() still captures the names, but calling fails

import anvil
sys = anvil.system("my_problem")

# Classic style: name, value, unit
sys.add("P0",   6.9e6, "Pa")
sys.add("T0",   3500,  "K")
sys.add("gamma", 1.25)           # dimensionless
sys.add("T0",   Q(3500, "K"))   # Quantity object

# kwargs style (unit inferred from Q)
sys.add(T0=3500*K, P0=6.9*MPa, gamma=1.25)

# With metadata
sys.add("T0", 3500, "K",
        desc="Chamber stagnation temperature",
        bounds=(200, 5000))

# Chain-able
sys.add("M", 2.0).add("gamma", 1.4)

sys.set(P0=8e6)              # bare number: keeps existing unit (Pa)
sys.set(P0=Q(1000, "psi"))   # Q object: overrides value AND unit/dim
sys.set(P0=8e6, T0=3200, gamma=1.3)  # multiple at once

sys.add("T", 300, "K")
sys.set(T=400)     # OK, 400 K
sys.set(T=Q(200, "R"))  # OK, 111.11 K internally, displays as 200 R

sys.set(new_var=5)  # KeyError: 'new_var' not in system

sys.use(my_function)              # plain function
sys.use(my_relation)             # Relation object
sys.use(other_system)            # System, inherits its defaults, wraps as Relation
sys.use("isentropic_ratios")     # registry name (string lookup)
sys.use(my_adapter)              # Adapter object

# With name mapping
sys.use("isentropic_ratios", map={"M": "M_exit"})
# map = {relation_param_name: workspace_variable_name}

nozzle_a = anvil.S.rocket_nozzle.copy()
nozzle_b = anvil.S.rocket_nozzle.copy()
nozzle_a.set(P0=10e6)    # does not affect nozzle_b

warnings = sys.validate()   # raises ValidationError on errors; returns warning list

# Example: missing input
sys = anvil.system("bad")
sys.add("M", 2.0)
sys.use("isentropic_ratios")  # needs M and gamma
try:
    sys.validate()
except ValidationError as e:
    print(e)
# ValidationError: Validation failed:
#   * 'isentropic_ratios' needs 'gamma' -- not provided.

result = sys.solve_forward()

result = sys.solve_gauss_seidel(
    relaxation=0.8,    # under-relaxation: 0 < ω ≤ 1. Lower = more stable, slower.
    max_iter=200,      # iteration limit
    rtol=1e-6,         # relative tolerance on max workspace change
    monitor=True,      # print live residuals per iteration
    verbose=False,     # brief residual print (less output than monitor)
)

  iter    0  |  residual = 2.0000e+00  |  t = 0.000s
  iter    1  |  residual = 2.3426e-01  |  t = 0.001s
  ...
  iter   32  |  residual = 6.9048e-11  |  t = 0.005s
  converged in 33 iterations

result = sys.solve_newton(
    max_iter=50,
    rtol=1e-10,
    verbose=True,
)

result = sys.solve()                   # auto: forward or gauss_seidel
result = sys.solve(method="gauss_seidel", relaxation=0.7)
result = sys.solve(method="newton")
result = sys.solve(method="forward")

hist = sys.history()
# [{"iteration": 0, "residual": 0.5, "wallclock": 0.001, "variables": {...}}, ...]

for h in hist:
    print(h["iteration"], h["residual"], h["wallclock"])

result = sys.solve_forward()

result["thrust"]              # Quantity(32143.09, "N")
result["thrust"].value        # 32143.09  (display value)
result["thrust"].si           # 32143.09  (SI float)
result["thrust"].to("kN")     # Quantity(32.143, "kN")
result["thrust"].to("lbf").value  # 7224.1

result.summary()              # formatted table
result.summary(keys=["thrust", "Isp"])  # subset

result.keys()                 # all variable names
result.to_dict()              # {name: display_value}
result.to_dict(si=True)       # {name: SI_float}
result.to_csv("out.csv")      # variable,value,unit per row
result.to_json("out.json")    # {name: {value, unit}}

nozzle_rel = nozzle.as_relation(
    inputs=["P0", "T0", "gamma", "R_gas", "A_throat", "A_exit", "P_amb"],
    outputs=["thrust", "Isp", "mdot"],
)

rocket_stage = anvil.system("full_rocket")
rocket_stage.add("P0", 10e6, "Pa")
rocket_stage.add("T0", 3500, "K")
# ... other inputs ...
rocket_stage.use(nozzle_rel)
rocket_stage.use(tsiolkovsky_rel)  # uses mdot from nozzle
result = rocket_stage.solve_forward()

opt = sys.optimize(
    objective,         # str, output quantity name to optimize
    design_vars,       # dict, {var_name: (lo, hi)}
    minimize=True,     # False = maximize
    method="differential_evolution",
    seed=None,         # int, random seed for reproducibility
    tol=1e-6,
    maxiter=1000,
    verbose=False,     # print progress every 25 evals
    **solver_kwargs,   # passed to System.solve() for every evaluation
)

opt.x           # dict: {var_name: optimal_value_in_display_units}
opt.fun         # objective value at optimum
opt.success     # bool
opt.message     # str, optimizer message
opt.nit         # optimizer iterations
opt.nfev        # total System.solve() calls performed

opt["thrust"]   # Quantity at optimum, same as result["thrust"]
"Isp" in opt    # bool, check if key exists in optimum result

opt.summary()   # print design vars + full result table

import numpy as np
import anvil

# --- Maximize nozzle thrust ---
nozzle = anvil.S.rocket_nozzle.copy()
nozzle.set(P0=8e6, T0=3200, gamma=1.25, R_gas=400, P_amb=101325)

opt = nozzle.optimize(
    objective="thrust",
    design_vars={
        "A_throat": (0.002, 0.030),  # bounds in m² (declared unit)
        "A_exit":   (0.010, 0.300),
    },
    minimize=False,
    method="differential_evolution",
    seed=42,
    maxiter=400,
    verbose=True,
)

opt.summary()
print(f"Thrust: {opt.fun/1000:.2f} kN")
print(f"A_throat: {opt.x['A_throat']*1e4:.1f} cm²")
print(f"Isp: {float(opt['Isp'].value):.1f} s")   # subscript access

# --- Gradient method for smooth landscapes ---
opt_fast = nozzle.optimize(
    objective="Isp",
    design_vars={"A_throat": (0.005, 0.020), "A_exit": (0.05, 0.25)},
    minimize=False,
    method="L-BFGS-B",
    maxiter=200,
)

# --- Coupled system: pass solver_kwargs ---
hx = anvil.system("heat_exchanger")
# ... setup ...
opt_hx = hx.optimize(
    objective="Q_duty",
    design_vars={"mdot_h": (0.1, 2.0), "mdot_c": (0.1, 2.0)},
    minimize=False,
    method="differential_evolution",
    # These go to System.solve():
    method_solve="gauss_seidel",  # ← note: solver_kwargs cannot reuse 'method'
    relaxation=0.7,
    max_iter=200,
)

print(sys.info())
# System: my_nozzle
#   Inputs:
#     P0                    6.9e+06 Pa
#     T0                    3500 K
#   Relations:
#     nozzle_area_ratio
#     area_mach_supersonic
#     ...

repr(sys)
# System('rocket_nozzle', 7 vars, 8 relations)

result = anvil.solve(isentropic, M=2.0, gamma=1.4)
result["T0_T"]   # 1.8

result = anvil.solve("normal_shock", M1=2.5)
result["M2"]     # ~0.513

sys.add("T0_T", 999.0)       # declared as input
sys.use("isentropic_ratios")  # also produces T0_T

sys.validate()
# WARNING: variable(s) declared via .add() are also produced by a relation, 
#   declared value will be overwritten after solve: ['T0_T']
#   (This is intentional for iterative initial guesses; for forward
#   systems it may indicate a naming mismatch.)

sys.add("extra_var", 42.0)  # not an input to any relation

sys.validate()
# WARNING: variable(s) declared via .add() are not used by any relation: ['extra_var']
#   (Possible typo or unused parameter.)

from anvil import solvers
# or
import anvil
anvil.solvers.find_root(...)

x = solvers.find_root(func, x0=None, bracket=None, method="auto", tol=1e-12, maxiter=100)

import numpy as np

# Brent's method, needs sign change in bracket
f = lambda x: x**3 - 2*x - 5
root = solvers.find_root(f, bracket=(2.0, 3.0))
# root = 2.0945514815  (exact to 1e-12)

# Newton's method, needs initial guess
root = solvers.find_root(f, x0=2.5, method="newton")
# root = 2.0945514815

# Used in built-in RSQs (e.g., area_mach_supersonic)
def residual(M, area_ratio=8.0, gamma=1.4):
    t = (2/(gamma+1))*(1+(gamma-1)/2*M**2)
    return (1/M)*t**((gamma+1)/(2*(gamma-1))) - area_ratio
M = solvers.find_root(residual, bracket=(1.001, 30.0))
# M = 3.9155...