"""R1 — Perch vs hover energy model.
Momentum-theory hover power with a figure of merit, versus sensor-only power when perched.
Outputs research/out/perch_vs_hover.json. Deterministic; no random inputs."""
import json, math
import numpy as np

G = 9.81

def air_density(alt_m, temp_c):
    # ISA pressure with altitude, ideal gas with actual temperature
    p = 101325 * (1 - 2.25577e-5 * alt_m) ** 5.25588
    return p / (287.05 * (temp_c + 273.15))

def hover_power(mass_kg, disk_area_m2, rho=1.225, fm=0.65, eta=0.85):
    T = mass_kg * G
    p_ideal = T ** 1.5 / math.sqrt(2 * rho * disk_area_m2)
    return p_ideal / (fm * eta)

out = {}

# 1) Reference airframe (Kestrel-class, 41 kg, 4 x 0.80 m rotors)
m, A = 41.0, 4 * math.pi * 0.40 ** 2
P_h = hover_power(m, A)
batt_wh = P_h * (38 / 60)            # energy that yields the published 38 min hover
out["reference"] = {"mass_kg": m, "disk_area_m2": round(A, 3), "hover_power_w": round(P_h), "usable_energy_wh": round(batt_wh),
                    "disk_loading_n_m2": round(m * G / A, 1)}
perch_loads = [15, 30, 45, 60, 90, 120]
out["reference"]["perch"] = [{"sensor_w": s, "perch_hours": round(batt_wh / s, 1), "ratio_vs_hover": round(P_h / s)} for s in perch_loads]

# 2) Parametric family: mass 2..45 kg, constant disk loading 200 N/m2, battery 35% of mass, 180 Wh/kg usable
masses = np.linspace(2, 45, 44)
fam = []
for mk in masses:
    area = mk * G / 200.0
    ph = hover_power(mk, area)
    e = mk * 0.35 * 180
    fam.append({"mass": round(float(mk), 1), "hover_min": round(e / ph * 60, 1),
                "perch_h_15w": round(e / 15, 1), "perch_h_45w": round(e / 45, 1), "perch_h_90w": round(e / 90, 1)})
out["family"] = fam

# 3) Density sensitivity: hover endurance at altitude / temperature, relative to ISA sea level
cases = [(0, 15), (0, 35), (1000, 25), (2000, 15), (2000, 35), (3000, 10)]
base = hover_power(m, A, rho=air_density(0, 15))
out["density"] = [{"alt_m": a, "temp_c": t, "rho": round(air_density(a, t), 3),
                   "endurance_rel": round(base / hover_power(m, A, rho=air_density(a, t)), 3)} for a, t in cases]

# 4) Duty cycle of a perch-and-shift cycle: time-on-watch fraction vs relocation interval
#    each relocation: 9 s lift + transit (distance/22 m/s) + 30 s perch settle; watch is lost only while airborne
#    (in practice EO/IR can watch in transit, which we conservatively ignore)
def tow(reloc_interval_h, dist_m):
    blind = 9 + dist_m / 22 + 30
    return 1 - blind / (reloc_interval_h * 3600)
out["duty"] = [{"interval_h": h, "tow_500m": round(tow(h, 500), 5), "tow_3km": round(tow(h, 3000), 5)} for h in [0.25, 0.5, 1, 2, 4, 8]]

json.dump(out, open("out/perch_vs_hover.json", "w"), indent=1)
print(json.dumps(out["reference"], indent=1)); print(out["density"]); print(out["duty"]); print(fam[0], fam[-1])