# Volumen de mezcla esperado — dos estimaciones independientes de la calibración hidráulica (revisión A1 §7: "restringir volumen").
# (1) Balance lluvia–escorrentía:  V_w = P·A·C,  V_m = V_w/(1 − C_v)
#     P: lluvia del pulso que generó el evento (DMC Tobalaba 15–21 h local: 13,8 mm; prensa: ~16–20 mm en ~1 h; factor orográfico
#        precordillera 1,0–2,0 sobre Tobalaba ~640 m s.n.m. vs cuenca media ~1.350 m) → P ∈ [12, 35] mm
#     A: cuenca H1 (cuenca.json);  C: escorrentía con suelo saturado (49 mm el 6–7 oct, año +81 %) y quemado → C ∈ [0,4, 0,9]
#     C_v total ∈ [0,30, 0,50]
# (2) Modelo empírico USGS para cuencas quemadas (Gartner, Cannon & Santi 2014, "emergency assessment"):
#     ln V = 4,22 + 0,39 √i15 + 0,36 ln B_mh + 0,13 √R      (V m³, i15 mm/h, B_mh km² quemados severidad moderada-alta, R m)
#     incendio dic-2025/ene-2026 ~1.000 ha en San Carlos de Apoquindo (prensa) → fracción quemada de la cuenca supuesta 0,4–1,0
#     i15 ∈ [25, 60] mm/h (20 mm en ~1 h con picos convectivos y granizo)
# Monte Carlo uniforme sobre los rangos; se reportan percentiles. Son rangos supuestos, no medidos: sirven de previa, no de verdad.
import json, numpy as np
rng = np.random.default_rng(20261008); N = 200000
cu = json.load(open('cuenca.json'))['H1']; A = cu['area_km2'] * 1e6; R = cu['relieve_R_m']
P = rng.uniform(0.012, 0.035, N); C = rng.uniform(0.4, 0.9, N); Cv = rng.uniform(0.30, 0.50, N)
Vm1 = P * A * C / (1 - Cv)
i15 = rng.uniform(25, 60, N); Bmh = cu['area_km2'] * rng.uniform(0.4, 1.0, N)
Vm2 = np.exp(4.22 + 0.39 * np.sqrt(i15) + 0.36 * np.log(Bmh) + 0.13 * np.sqrt(R))
q = lambda v: {f'p{p}': round(float(np.percentile(v, p))) for p in (5, 25, 50, 75, 95)}
# caudal pico empírico (Mizuyama et al. 1992): flujo de barro Qp = 0,0188 V^0,79 ; granular Qp = 0,135 V^0,78
V50 = float(np.percentile(Vm1, 50))
out = {'cuenca_H1': cu, 'V_lluvia_escorrentia_m3': q(Vm1), 'V_usgs_gartner2014_m3': q(Vm2),
       'Qp_mizuyama_barro_m3s_para_V50': round(0.0188 * V50 ** 0.79, 1), 'Qp_mizuyama_granular_m3s_para_V50': round(0.135 * V50 ** 0.78, 1),
       'rango_previo_para_calibrar': {'V_m3': [round(float(np.percentile(np.r_[Vm1, Vm2], 10))), round(float(np.percentile(np.r_[Vm1, Vm2], 90)))]}}
json.dump(out, open('res/volumen_previo.json', 'w'), indent=1, ensure_ascii=False); print(json.dumps(out, indent=1, ensure_ascii=False))