# Láminas y tablas de parámetros del paper, generadas desde las mismas constantes y datos que usa el modelo:
#   paper/tablas/parametros.tex  — constantes físicas, reología (coeficientes por muestra y valores derivados), fricción, malla, entrada
#   paper/tablas/datos.tex       — datos de entrada utilizados (fuente, resolución, uso)
#   paper/figs/ubicacion.pdf     — mapa de ubicación: cuencas, cauces, puntos de entrada, embocadura, sitios y dominio
#   paper/figs/reologia.pdf      — τ_y y η en función de C_v para las muestras, ρ_m(C_v)
#   paper/figs/forzamiento.pdf   — lluvia DMC (6 h) y hidrogramas triangulares ensayados
import json, math, numpy as np, matplotlib, rasterio
matplotlib.use('Agg'); import matplotlib.pyplot as plt
from matplotlib.colors import LightSource
from pyproj import Transformer
plt.rcParams.update({'font.family': 'serif', 'font.size': 8.5, 'axes.titlesize': 9})
SUELOS = {'Glenwood 1': (0.00283, 23.0, 0.0345, 20.1), 'Glenwood 2a': (0.0648, 6.20, 0.0765, 16.9), 'Glenwood 2b': (0.648, 6.20, 0.0765, 16.9),
          'Glenwood 2c': (0.648, 8.20, 0.0765, 16.9), 'Aspen 1': (0.0360, 22.1, 0.181, 25.7), 'Aspen 2': (0.0538, 14.5, 2.72, 10.4)}
RW, RS, G = 1000.0, 2650.0, 9.8
tau = lambda a2, b2, cv: 0.1 * a2 * math.exp(b2 * cv); eta = lambda a1, b1, cv: 0.1 * a1 * math.exp(b1 * cv)
num = lambda x, d=2: f'{x:.{d}f}'.replace('.', '{,}')
# ---------- tabla de parámetros
L = [r'\begin{table}[htbp]\centering\small', r'\caption{Parámetros físicos y numéricos aplicados en las simulaciones. Los coeficientes reológicos son los de la tabla de muestras de O\'Brien y Julien recopilada en el manual de FLO-2D \citep{flo2d2021}; $\tau_y=0{,}1\,\alpha_2e^{\beta_2C_v}$ en Pa y $\eta=0{,}1\,\alpha_1e^{\beta_1C_v}$ en Pa\,s.}\label{tab:parametros}',
     r'\begin{tabular}{lrrrrrrr}\toprule', r'Muestra & $\alpha_1$ & $\beta_1$ & $\alpha_2$ & $\beta_2$ & \multicolumn{3}{c}{$\tau_y$ (Pa) / $\eta$ (Pa\,s) con $C_v=$}\\\cmidrule(l){6-8}', r' & & & & & 0,35 & 0,40 & 0,45\\\midrule']
for k, (a1, b1, a2, b2) in SUELOS.items():
    L.append(f'{k} & {num(a1, 4)} & {num(b1)} & {num(a2, 4)} & {num(b2)} & ' + ' & '.join(f'{num(tau(a2, b2, c), 1)} / {num(eta(a1, b1, c), 2)}' for c in (0.35, 0.40, 0.45)) + r'\\')
L += [r'\midrule', r'\multicolumn{8}{l}{\textbf{Mezcla y constantes}}\\',
      rf'\multicolumn{{5}}{{l}}{{Densidad del agua $\rho_w$ / de los sólidos $\rho_s$}} & \multicolumn{{3}}{{r}}{{\SI{{1000}}{{kg/m^3}} / \SI{{2650}}{{kg/m^3}}}}\\',
      r'\multicolumn{5}{l}{Densidad de la mezcla $\rho_m=\rho_w+C_v(\rho_s-\rho_w)$, $C_v=0{,}35/0{,}40/0{,}45$} & \multicolumn{3}{r}{' + ' / '.join(f'{RW + c * (RS - RW):.0f}' for c in (0.35, 0.40, 0.45)) + r' kg/m$^3$}\\',
      r'\multicolumn{5}{l}{Gravedad $g$} & \multicolumn{3}{r}{\SI{9,8}{m/s^2}}\\',
      r'\multicolumn{5}{l}{Resistencia laminar $K$ (calles / resto)} & \multicolumn{3}{r}{100 / 400}\\',
      r'\multicolumn{5}{l}{Manning turbulento-dispersivo $n_{td}$ (calles / resto)} & \multicolumn{3}{r}{0,025 / 0,050 s\,m$^{-1/3}$}\\',
      r'\multicolumn{8}{l}{\textbf{Entrada, malla y cálculo}}\\',
      r'\multicolumn{5}{l}{Volumen de mezcla $V$ / caudal pico $Q_p$ / tiempo al pico $t_p$} & \multicolumn{3}{r}{80--300 mil m$^3$ / 100--250 m$^3$/s / 600 s}\\',
      r'\multicolumn{5}{l}{Velocidad de entrada / radio de cada fuente} & \multicolumn{3}{r}{2 m/s / 14 m}\\',
      r'\multicolumn{5}{l}{Área máxima de triángulo (corredor / resto); triángulos} & \multicolumn{3}{r}{40 / 600 m$^2$; $\approx$430\,000}\\',
      r'\multicolumn{5}{l}{Grilla del terreno; rebaje de calles; altura supuesta de edificios} & \multicolumn{3}{r}{4 m; 0,15 m; 6 m}\\',
      r'\multicolumn{5}{l}{Número de Courant (cribado / lluvia / final)} & \multicolumn{3}{r}{1 / 0,5 / $\le$0,25}\\',
      r'\multicolumn{5}{l}{Umbral de mojado para la evaluación $h^*$} & \multicolumn{3}{r}{5 cm (sensibilidad 2 y 10 cm)}\\',
      r'\bottomrule\end{tabular}\end{table}']
open('paper/tablas/parametros.tex', 'w').write('\n'.join(L) + '\n')
# ---------- tabla de datos
D = [(r'Terreno desnudo', r'FABDEM V1-2 \citep{fabdem2023}', r'1 segundo de arco ($\approx$\SI{30}{m}), remuestreado a \SI{4}{m}', 'Cotas del suelo y cuencas'),
     (r'Superficie', r'Copernicus GLO-30 \citep{copernicus_glo30}', r'\SI{30}{m}', 'Altura de edificios grandes'),
     (r'Terreno alternativo', r'Fotogrametría SPOT 2014 \citep{hussain2021}', r'0,07 puntos/m$^2$ ($\approx$\SI{8}{}--\SI{10}{m})', 'Prueba de sensibilidad'),
     (r'Edificios, calles, cauces', r'OpenStreetMap \citep{osm2026}', 'vectorial', '1.103 bloques, calles, trazado del cauce'),
     (r'Lluvia', r'DMC, estación 330019 \citep{dmc2026}', r'horaria y sinóptica', 'Volumen previo'),
     (r'Observaciones', r'Prensa, municipio, Senapred, Sernageomin', r'puntual', 'Evaluación (no entra al modelo)')]
T = [r'\begin{table}[htbp]\centering\small', r'\caption{Datos de entrada utilizados.}\label{tab:datos}', r'\begin{tabular}{p{.17\textwidth}p{.3\textwidth}p{.22\textwidth}p{.22\textwidth}}\toprule', r'Dato & Fuente & Resolución & Uso\\\midrule']
T += [' & '.join(r) + r'\\' for r in D] + [r'\bottomrule\end{tabular}\end{table}']
open('paper/tablas/datos.tex', 'w').write('\n'.join(T) + '\n')
# ---------- lámina de ubicación
tf = Transformer.from_crs('EPSG:4326', 'EPSG:32719', always_xy=True)
from rasterio.warp import reproject, Resampling
T0_ = np.load('terreno.npz'); _x0 = float(T0_['X0']) - 600; _y0 = float(T0_['Y0']) - 600; _x1, _y1 = 366500.0, 6305000.0; _R = 20.0
_W, _H = int((_x1 - _x0) / _R), int((_y1 - _y0) / _R); _tr = rasterio.transform.from_origin(_x0, _y1, _R, _R); Z = np.zeros((_H, _W), 'float32')
with rasterio.open('/home/linuxoficina/Proyectos/dem/fabdem/S34W071_FABDEM_V1-2.tif') as _d:
    reproject(rasterio.band(_d, 1), Z, dst_transform=_tr, dst_crs='EPSG:32719', resampling=Resampling.bilinear)
class _B: left, right, bottom, top = _x0, _x1, _y0, _y1
b = _B()
T0 = np.load('terreno.npz'); X0, Y0, DX, NX, NY = [float(T0[k]) for k in ('X0', 'Y0', 'DX', 'NX', 'NY')]
fig, ax = plt.subplots(figsize=(7.2, 5.0))
ax.imshow(LightSource(315, 40).hillshade(Z, vert_exag=2, dx=20, dy=20), cmap='gray', extent=(b.left, b.right, b.bottom, b.top), alpha=0.9)
cs = ax.contour(np.flipud(Z), levels=range(800, 3000, 200), colors='0.35', linewidths=0.25, extent=(b.left, b.right, b.bottom, b.top))
o = json.load(open('osm.json'))
for c in o['calles']:
    if c['k'] in ('primary', 'secondary', 'tertiary'):
        x, y = tf.transform([q[1] for q in c['pts']], [q[0] for q in c['pts']]); ax.plot(x, y, c='0.15', lw=0.35)
for c in o['cauces']:
    x, y = tf.transform([q[1] for q in c['pts']], [q[0] for q in c['pts']]); ax.plot(x, y, c='#1f6fb5', lw=0.7 if not c.get('tunnel') else 1.2, ls='-' if not c.get('tunnel') else ':')
COL = {'H1': '#d95f02', 'S1': '#c51b7d', 'S2': '#7570b3'}
for k, col in COL.items():
    for f in json.load(open(f'cuenca_{k}.geojson'))['features']:
        for ring in f['geometry']['coordinates']: r = np.array(ring); ax.plot(r[:, 0], r[:, 1], c=col, lw=1.3)
    cu = json.load(open('cuenca.json'))[k]; ax.plot([], [], c=col, lw=1.3, label=f'Cuenca {k} ({cu["area_km2"]:.2f} km²)'.replace('.', ','))
ax.add_patch(plt.Rectangle((X0, Y0), NX * DX, NY * DX, fill=False, ec='k', lw=1, ls='--'))
ax.text(X0 + 60, Y0 + NY * DX + 60, 'dominio del modelo', fontsize=7)
P = {'H1': (-33.40030, -70.49960, '^', '#d95f02'), 'S1': (-33.40703, -70.49092, '^', '#c51b7d'), 'S2': (-33.40870, -70.49990, '^', '#7570b3'),
     'Embocadura Bulnes Correa': (-33.40586, -70.51609, 'o', 'red'), 'GB × S. C. de Apoquindo': (-33.4018, -70.5099, 'o', 'red'),
     'UC Christus': (-33.40014, -70.50874, 's', 'w'), 'Clínica U. Andes': (-33.40727, -70.50373, 's', 'w'), 'Duoc UC': (-33.40063, -70.50540, 's', 'w'),
     'Campus U. Andes': (-33.40424, -70.50660, 's', 'w'), 'Padre Hurtado × GB': (-33.40918, -70.54224, 'D', 'y'), 'Rotonda Atenas': (-33.41768, -70.55272, 'D', 'y')}
for k, (la, lo, mk, col) in P.items():
    x, y = tf.transform(lo, la); ax.plot(x, y, mk, ms=5, mfc=col, mec='k', mew=0.6, zorder=5)
    if k not in ('UC Christus', 'Duoc UC'): ax.annotate(k, (x, y), (4, 4), textcoords='offset points', fontsize=6.3, zorder=6, bbox=dict(fc='w', ec='none', alpha=0.7, pad=0.6))
ax.plot([], [], c='#1f6fb5', lw=0.7, label='cauce abierto (OSM)'); ax.plot([], [], c='#1f6fb5', lw=1.2, ls=':', label='tramo entubado (OSM)')
ax.plot([], [], '^', mfc='w', mec='k', label='punto de entrada'); ax.plot([], [], 's', mfc='w', mec='k', label='sitio afectado'); ax.plot([], [], 'o', mfc='red', mec='k', label='nodo del drenaje')
ax.set_xlim(X0 - 400, 366500); ax.set_ylim(Y0 - 300, 6304700); ax.set_aspect('equal')
ax.set_xlabel('Este UTM 19S (m)'); ax.set_ylabel('Norte UTM 19S (m)'); ax.ticklabel_format(style='plain'); ax.tick_params(labelsize=6.5)
ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.13), ncol=4, fontsize=6.2, frameon=False); fig.tight_layout(); fig.savefig('paper/figs/ubicacion.pdf', dpi=250); plt.close(fig)
# ---------- lámina de reología
cv = np.linspace(0.25, 0.55, 100)
fig, axs = plt.subplots(1, 3, figsize=(7.2, 2.5))
for k, (a1, b1, a2, b2) in SUELOS.items():
    if k in ('Glenwood 2a', 'Glenwood 2c'): continue
    axs[0].semilogy(cv, [tau(a2, b2, c) for c in cv], label=k); axs[1].semilogy(cv, [eta(a1, b1, c) for c in cv], label=k)
for ax_ in axs[:2]: ax_.axvspan(0.35, 0.45, color='0.85', zorder=0); ax_.set_xlabel('$C_v$')
axs[0].set_ylabel(r'$\tau_y$ (Pa)'); axs[1].set_ylabel(r'$\eta$ (Pa·s)'); axs[0].legend(fontsize=6.5, frameon=False)
axs[2].plot(cv, RW + cv * (RS - RW), c='k'); axs[2].axvspan(0.35, 0.45, color='0.85'); axs[2].set_xlabel('$C_v$'); axs[2].set_ylabel(r'$\rho_m$ (kg/m³)')
fig.tight_layout(); fig.savefig('paper/figs/reologia.pdf'); plt.close(fig)
# ---------- lámina de forzamiento
fig, axs = plt.subplots(1, 2, figsize=(7.2, 2.5))
fechas = ['6-oct 18', '7-oct 00', '7-oct 06', '7-oct 12', '7-oct 18', '8-oct 00', '8-oct 06', '8-oct 12', '8-oct 18', '9-oct 00', '9-oct 06']
pp = [np.nan, np.nan, 49.2, 12.6, 5.2, 0.2, 0.0, 0.0, 3.2, 13.8, 17.8]
axs[0].bar(range(len(pp)), pp, color='#1f6fb5'); axs[0].set_xticks(range(len(pp))); axs[0].set_xticklabels(fechas, rotation=60, fontsize=6)
axs[0].set_ylabel('lluvia en 6 h (mm)'); axs[0].set_title('Estación Tobalaba (DMC), totales sinópticos UTC', fontsize=7.5)
axs[0].annotate('aluvión ~20 UTC', (8.5, 14), (5.2, 30), fontsize=6.5, arrowprops=dict(arrowstyle='->', lw=0.6))
for V, Qp in ((80000, 100), (150000, 250), (300000, 250)):
    TB = 2 * V / Qp; t = np.array([0, 600, TB]); axs[1].plot(t / 60, [0, Qp, 0], label=f'{V // 1000} mil m³, {Qp} m³/s')
axs[1].set_xlabel('tiempo (min)'); axs[1].set_ylabel('caudal de mezcla (m³/s)'); axs[1].legend(fontsize=6.5, frameon=False); axs[1].set_title('Hidrogramas triangulares ensayados', fontsize=7.5)
fig.tight_layout(); fig.savefig('paper/figs/forzamiento.pdf'); plt.close(fig)
print('láminas y tablas listas')