import sys
IN_COLAB = 'google.colab' in sys.modules

if IN_COLAB:
    !pip install -q rasterio pyproj

import os
import glob
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import matplotlib.colors as mcolors
from matplotlib.colors import TwoSlopeNorm, BoundaryNorm
from matplotlib.patches import FancyArrowPatch
from pathlib import Path

import rasterio
from pyproj import Transformer

# Directorios
WORK_DIR   = Path('/content') if IN_COLAB else Path('.')
OUTPUT_DIR = WORK_DIR / 'output_1979'
DATA_DIR   = Path('../data') if not IN_COLAB else WORK_DIR / 'data'
DEM_PATH   = WORK_DIR / 'tumaco_dem_utm_sim.tif'

print(f"Directorio de resultados: {OUTPUT_DIR}")

# Celda 1 — Cargar DEM (auto-genera si no existe en esta sesión)
from scipy.ndimage import gaussian_filter
from rasterio.crs import CRS
from rasterio.transform import from_bounds
from rasterio.warp import calculate_default_transform, reproject, Resampling
import re, requests

def _t_from_name(path):
    m = re.search(r'h-([0-9.]+)\.np$', str(path))
    return float(m.group(1)) if m else 0.0

if not DEM_PATH.exists():
    print('DEM no encontrado — generando desde GEBCO/ETOPO...')
    LAT_MIN, LAT_MAX, LON_MIN, LON_MAX = 0.5, 3.5, -79.5, -77.0
    CRS_SRC  = CRS.from_epsg(4326)
    CRS_DEST = CRS.from_epsg(32618)
    geo_path = WORK_DIR / 'tumaco_bathy_geo.tif'
    bathy_ok = False

    # Intento 1: GEBCO 2023
    try:
        url = (f'https://download.gebco.net/a/gebco_2023'
               f'?lat1={LAT_MIN}&lat2={LAT_MAX}&lon1={LON_MIN}&lon2={LON_MAX}&format=netcdf4')
        r = requests.get(url, timeout=120)
        if r.status_code == 200 and len(r.content) > 5000:
            import xarray as xr, io
            ds      = xr.open_dataset(io.BytesIO(r.content))
            varname = [v for v in ds.data_vars if 'elev' in v.lower() or 'z' in v.lower()][0]
            da      = ds[varname]
            lat_dim = [d for d in da.dims if 'lat' in d.lower()][0]
            lon_dim = [d for d in da.dims if 'lon' in d.lower()][0]
            if da[lat_dim].values[0] < da[lat_dim].values[-1]:
                da = da.isel({lat_dim: slice(None, None, -1)})
            arr  = da.values.astype('float32')
            lats = da[lat_dim].values
            lons = da[lon_dim].values
            tf_g = from_bounds(lons.min(), lats.min(), lons.max(), lats.max(),
                               arr.shape[1], arr.shape[0])
            with rasterio.open(geo_path, 'w', driver='GTiff',
                               height=arr.shape[0], width=arr.shape[1],
                               count=1, dtype='float32', crs=CRS_SRC, transform=tf_g) as dst:
                dst.write(arr, 1)
            bathy_ok = True
            print('  ✓ GEBCO 2023')
    except Exception as e:
        print(f'  ✗ GEBCO: {e}')

    # Intento 2: ETOPO 2022
    if not bathy_ok:
        try:
            nc = int((LON_MAX - LON_MIN) / 0.004167)
            nr = int((LAT_MAX - LAT_MIN) / 0.004167)
            url2 = ('https://gis.ngdc.noaa.gov/arcgis/rest/services/DEM_mosaics/'
                    'DEM_global_mosaic/ImageServer/exportImage'
                    f'?bbox={LON_MIN},{LAT_MIN},{LON_MAX},{LAT_MAX}'
                    f'&bboxSR=4326&size={nc},{nr}&imageSR=4326&format=tiff&pixelType=F32&f=image')
            r2 = requests.get(url2, timeout=120)
            r2.raise_for_status()
            with open(geo_path, 'wb') as f:
                f.write(r2.content)
            bathy_ok = True
            print('  ✓ ETOPO 2022')
        except Exception as e:
            print(f'  ✗ ETOPO: {e}')

    # Intento 3: DEM sintético
    if not bathy_ok:
        print('  ⚠ DEM sintético (solo pruebas)')
        res    = 0.004167
        lons_s = np.arange(LON_MIN, LON_MAX, res)
        lats_s = np.arange(LAT_MIN, LAT_MAX, res)
        nr, nc = len(lats_s), len(lons_s)
        xr2 = (lons_s - LON_MIN) / (LON_MAX - LON_MIN)
        z   = np.where(xr2 < 0.70, -4000 + 3800*(xr2/0.70),
              np.where(xr2 < 0.88, -200 + 200*((xr2-0.70)/0.18), 10.0))
        dem_s  = np.tile(z, (nr, 1)).astype('float32')
        dem_s += gaussian_filter(np.random.randn(nr, nc)*50, sigma=5)
        tf_s   = from_bounds(LON_MIN, LAT_MIN, LON_MAX, LAT_MAX, nc, nr)
        with rasterio.open(geo_path, 'w', driver='GTiff', height=nr, width=nc,
                           count=1, dtype='float32', crs=CRS_SRC, transform=tf_s) as dst:
            dst.write(dem_s, 1)

    # Reproyectar a UTM 18N
    utm_path = WORK_DIR / 'tumaco_dem_utm.tif'
    with rasterio.open(geo_path) as src:
        tf_utm, w_utm, h_utm = calculate_default_transform(
            src.crs, CRS_DEST, src.width, src.height, *src.bounds)
        meta = src.meta.copy()
        meta.update({'crs': CRS_DEST, 'transform': tf_utm,
                     'width': w_utm, 'height': h_utm, 'dtype': 'float32'})
        with rasterio.open(utm_path, 'w', **meta) as dst:
            reproject(source=rasterio.band(src, 1), destination=rasterio.band(dst, 1),
                      src_transform=src.transform, src_crs=src.crs,
                      dst_transform=tf_utm, dst_crs=CRS_DEST,
                      resampling=Resampling.bilinear)

    # Invertir eje X (orientación W→E compatible con el solver)
    with rasterio.open(utm_path) as src:
        arr_u  = src.read(1)
        old_tf = src.transform
        meta2  = src.meta.copy()
    arr_flip = arr_u[:, ::-1]
    meta2['transform'] = rasterio.transform.from_origin(
        west=old_tf.c + old_tf.a*arr_u.shape[1], north=old_tf.f,
        xsize=-old_tf.a, ysize=-old_tf.e)
    with rasterio.open(DEM_PATH, 'w', **meta2) as dst:
        dst.write(arr_flip, 1)
    print(f'  ✓ DEM guardado: {DEM_PATH}')

# Cargar DEM
with rasterio.open(str(DEM_PATH)) as src:
    dem = src.read(1)
    dem_transform = src.transform
    dem_crs = src.crs
    RESOLUTION = abs(src.transform.a)
    NROWS, NCOLS = src.height, src.width

print(f"DEM: {NROWS}×{NCOLS} px | Resolución: {RESOLUTION:.0f} m/px")

# Cargar snapshots — sort numérico para h-0 … h-1800
output_files = sorted(glob.glob(str(OUTPUT_DIR / 'h-*.np')), key=_t_from_name)

if not output_files:
    print("⚠ No se encontraron archivos de simulación en:", OUTPUT_DIR)
    print("  Ejecuta el Cuaderno 2 en esta sesión antes de continuar.")
    DEMO_MODE = True
else:
    DEMO_MODE = False
    print(f"✓ {len(output_files)} snapshots cargados | "
          f"{Path(output_files[0]).name} → {Path(output_files[-1]).name}")

# Construir diccionario tiempo → array
heightmaps = {}
if not DEMO_MODE:
    for fpath in output_files:
        t = _t_from_name(fpath)
        heightmaps[t] = np.load(fpath)
    t_range = sorted(heightmaps.keys())
    print(f"  Tiempo: {t_range[0]:.0f}s – {t_range[-1]:.0f}s")
else:
    # Modo demo
    t_range = list(range(0, 1860, 60))
    Ly = RESOLUTION * NCOLS
    for t in t_range:
        y_vec = np.linspace(0, Ly, NCOLS)
        c_wave = 200
        y_epi = Ly - 60_000
        y_wave = y_epi - c_wave * t
        lamb = 80_000
        k1, k2 = 28.416 / lamb**2, 256. / lamb**2
        eta = 2 * (8.0 * np.exp(-k1 * (y_vec - (y_wave + 0.3153*lamb))**2)
                   - (8.0/3) * np.exp(-k2 * (y_vec - y_wave)**2))
        depth = np.abs(np.minimum(dem, 0)).mean(axis=0).clip(1, None)
        amp = (depth.mean() / depth.clip(1))**0.25
        hmap = np.tile(eta * amp, (NROWS, 1)) - np.minimum(dem, 0)
        hmap = np.clip(hmap, 0.001, None)
        heightmaps[float(t)] = hmap.astype(np.float32)
    t_range = sorted(heightmaps.keys())

# Calcular η máxima en toda la simulación
# LSWE no inunda tierra: las celdas terrestres tienen h_arr=0 por construcción.
# Lo significativo es la elevación máxima de la superficie marina (η) en el océano.
# η = h_arr + min(dem, 0) = (H0+η) + b = η  para celdas oceánicas
# η = 0 + 0 = 0                              para celdas terrestres

max_eta = np.zeros_like(dem, dtype=np.float32)

for t, hmap in heightmaps.items():
    h_arr = hmap if hmap.shape == dem.shape else hmap[:dem.shape[0], :dem.shape[1]]
    eta = (h_arr + np.minimum(dem, 0)).astype(np.float32)
    max_eta = np.maximum(max_eta, eta)

ocean_mask_val = dem < 0
coast_zone_val = (dem < 0) & (dem > -200)

print(f"η máxima (todo el océano):    {max_eta[ocean_mask_val].max():.2f} m")
print(f"η máxima (zona costera <200m): {max_eta[coast_zone_val].max():.2f} m")

# Coordenadas de lugares de interés en UTM 18N
transformer = Transformer.from_crs(4326, 32618, always_xy=True)

places = {
    'Tumaco': (1.82, -78.76),
    'Bocagrande': (1.73, -78.94),
    'Epicentro\n1979': (1.598, -79.359),
}
places_utm = {
    name: transformer.transform(lon, lat)
    for name, (lat, lon) in places.items()
}

def utm_to_pixel(x_utm, y_utm, transform):
    col = int((x_utm - transform.c) / transform.a)
    row = int((y_utm - transform.f) / transform.e)
    return row, col

# Extensión geográfica del DEM en km (bordes del raster)
x_origin = dem_transform.c / 1000
y_origin = dem_transform.f / 1000
x_end = (dem_transform.c + dem_transform.a * NCOLS) / 1000
y_end = (dem_transform.f + dem_transform.e * NROWS) / 1000

# Grilla de coordenadas para ax.contour (centros de pixel, en km)
# Necesaria porque contour no acepta el parámetro extent de imshow
_x_c = (dem_transform.c + dem_transform.a * (np.arange(NCOLS) + 0.5)) / 1000
_y_c = (dem_transform.f + dem_transform.e * (np.arange(NROWS) + 0.5)) / 1000
X_ct, Y_ct = np.meshgrid(_x_c, _y_c)

ext = [x_origin, x_end, y_end, y_origin]  # extent común para imshow

eta_clim = float(np.nanpercentile(max_eta[coast_zone_val], 99)) if coast_zone_val.any() else 6.0
eta_clim = max(eta_clim, 3.0)

fig, axes = plt.subplots(1, 2, figsize=(16, 7))

# --- Panel izquierdo: η máxima en el océano ---
ax = axes[0]
ax.imshow(np.where(dem < 0, dem, np.nan), cmap='Blues', alpha=0.4,
          extent=ext, vmin=-5000, vmax=0)
ax.imshow(np.where(dem >= 0, dem, np.nan), cmap='Greys', alpha=0.35,
          extent=ext, vmin=0, vmax=200)

eta_masked = np.ma.masked_where((max_eta <= 0.1) | (dem >= 0), max_eta)
im = ax.imshow(eta_masked, cmap='YlOrRd', vmin=0.1, vmax=eta_clim,
               extent=ext, alpha=0.85)
plt.colorbar(im, ax=ax, label='η máxima (m sobre NMM)', shrink=0.8)

ax.contour(X_ct, Y_ct, dem, levels=[0], colors='k', linewidths=0.8)

for name, (x_utm, y_utm) in places_utm.items():
    is_epi = 'Epicentro' in name
    ax.plot(x_utm/1000, y_utm/1000, 'r*' if is_epi else 'ko',
            markersize=12 if is_epi else 7, zorder=5)
    ax.annotate(name, (x_utm/1000, y_utm/1000),
                textcoords='offset points', xytext=(5, 5),
                fontsize=8, color='darkred' if is_epi else 'black')

ax.set_xlabel('Easting UTM 18N (km)')
ax.set_ylabel('Northing UTM 18N (km)')
ax.set_title('Elevación máxima de la ola — Tsunami Tumaco 1979\n(30 minutos, solver LSWE)')

# --- Panel derecho: η en zona costera (<200 m de profundidad) ---
ax2 = axes[1]
im2 = ax2.imshow(np.where(coast_zone_val, max_eta, np.nan),
                 cmap='hot_r', vmin=0, vmax=eta_clim, extent=ext)
plt.colorbar(im2, ax=ax2, label='η máxima en zona costera (m)', shrink=0.8)
ax2.contour(X_ct, Y_ct, dem, levels=[0], colors='k', linewidths=0.5)
ax2.set_xlabel('Easting UTM 18N (km)')
ax2.set_title('Amplitud máxima de la ola\nen zona costera (<200 m profundidad)')

for name, (x_utm, y_utm) in places_utm.items():
    if 'Epicentro' not in name:
        ax2.plot(x_utm/1000, y_utm/1000, 'ko', markersize=7)
        ax2.annotate(name, (x_utm/1000, y_utm/1000),
                     textcoords='offset points', xytext=(5, 3), fontsize=8)

plt.tight_layout()
plt.savefig(WORK_DIR / 'max_inundation_map.png', dpi=200, bbox_inches='tight')
plt.show()
print(f"η_max zona costera: {float(np.nanmax(np.where(coast_zone_val, max_eta, np.nan))):.2f} m")
print(f"Mapa guardado en: {WORK_DIR}/max_inundation_map.png")

# Estaciones de observación con coordenadas y run-up histórico
stations = [
    {'name': 'Tumaco (muelle)', 'lat': 1.82,  'lon': -78.76, 'obs_runup': 5.0,
     'obs_arrival': 15, 'color': 'firebrick'},
    {'name': 'Bocagrande',       'lat': 1.73,  'lon': -78.94, 'obs_runup': 3.0,
     'obs_arrival': 18, 'color': 'darkorange'},
    {'name': 'Esmeraldas (ECU)', 'lat': 0.97,  'lon': -79.65, 'obs_runup': 4.0,
     'obs_arrival': 12, 'color': 'steelblue'},
    {'name': 'El Charco',        'lat': 2.47,  'lon': -78.11, 'obs_runup': 1.5,
     'obs_arrival': 30, 'color': 'darkgreen'},
]

# Extraer series de tiempo de los heightmaps
fig, axes = plt.subplots(2, 2, figsize=(14, 8), sharex=True)
fig.suptitle('Series de tiempo — Altura de la superficie del agua\nTsunami Tumaco 1979',
             fontsize=13)

for idx, (sta, ax) in enumerate(zip(stations, axes.flatten())):
    # Convertir coordenadas a píxel
    x_utm, y_utm = transformer.transform(sta['lon'], sta['lat'])
    try:
        row, col = utm_to_pixel(x_utm, y_utm, dem_transform)
        row = np.clip(row, 0, NROWS - 1)
        col = np.clip(col, 0, NCOLS - 1)
        elev = dem[row, col]

        # Extraer serie de tiempo
        times = np.array(sorted(heightmaps.keys()))
        h_vals = np.array([heightmaps[t][row, col] + min(dem[row, col], 0)
                           for t in times])
    except (IndexError, KeyError):
        # Si la estación está fuera del dominio, mostrar señal vacía
        times = np.array(sorted(heightmaps.keys()))
        h_vals = np.zeros_like(times)
        elev = 0.0

    # Graficar
    ax.plot(times / 60, h_vals, color=sta['color'], linewidth=2, label='Simulado')
    ax.axhline(0, color='k', linewidth=0.5, linestyle=':')
    ax.axvline(sta['obs_arrival'], color='gray', linewidth=1, linestyle='--',
               label=f"Arribo obs. ({sta['obs_arrival']} min)")
    ax.axhline(sta['obs_runup'], color=sta['color'], linewidth=1.5, linestyle='--',
               alpha=0.6, label=f"Run-up obs.: {sta['obs_runup']} m")

    ax.fill_between(times / 60, h_vals, 0,
                    where=h_vals > 0, alpha=0.2, color=sta['color'])
    ax.set_title(f"{sta['name']} (z={elev:.1f}m)", fontsize=10)
    ax.set_ylabel('Η (m)', fontsize=9)
    ax.set_ylim(-3, 8)
    ax.legend(fontsize=8, loc='upper right')
    ax.grid(True, alpha=0.3)

for ax in axes[1]:
    ax.set_xlabel('Tiempo desde el sismo (minutos)')

plt.tight_layout()
plt.savefig(WORK_DIR / 'time_series.png', dpi=150, bbox_inches='tight')
plt.show()
print(f"Series de tiempo guardadas en: {WORK_DIR}/time_series.png")

# Cargar observaciones históricas
try:
    obs_path = DATA_DIR / 'observaciones_1979.csv'
    df_obs = pd.read_csv(obs_path)
    print(f"Observaciones cargadas: {len(df_obs)} estaciones")
except FileNotFoundError:
    df_obs = pd.DataFrame({
        'estacion': ['Tumaco (muelle)', 'Tumaco (playa)', 'Bocagrande',
                     'San Juan', 'El Charco', 'Esmeraldas (ECU)',
                     'Buenaventura', 'Salinas (ECU)'],
        'lat': [1.82, 1.80, 1.73, 1.23, 2.47, 0.97, 3.88, -2.20],
        'lon': [-78.76, -78.80, -78.94, -77.97, -78.11, -79.65, -77.02, -80.97],
        'runup_m': [5.0, 4.5, 3.0, 2.5, 1.5, 4.0, 0.8, 1.2],
    })

# η máxima simulada cerca de cada estación (celdas oceánicas dentro de 3×3 px)
sim_runup = []
for _, row in df_obs.iterrows():
    try:
        x_utm, y_utm = transformer.transform(row['lon'], row['lat'])
        r, c = utm_to_pixel(x_utm, y_utm, dem_transform)
        r = np.clip(r, 0, NROWS - 1)
        c = np.clip(c, 0, NCOLS - 1)
        r0, r1 = max(0, r-2), min(NROWS, r+3)
        c0, c1 = max(0, c-2), min(NCOLS, c+3)
        # Usar max_eta (sea surface η) en zona cercana
        patch = max_eta[r0:r1, c0:c1]
        ocean_patch = patch[dem[r0:r1, c0:c1] < 0]
        max_val = float(ocean_patch.max()) if ocean_patch.size > 0 else float(patch.max())
    except Exception:
        max_val = np.nan
    sim_runup.append(max_val)

df_obs['sim_runup_m'] = sim_runup
df_val = df_obs.dropna(subset=['sim_runup_m'])
df_val = df_val[df_val['sim_runup_m'] > 0.01]   # excluir puntos fuera del dominio

if len(df_val) > 0:
    rmse = np.sqrt(((df_val['runup_m'] - df_val['sim_runup_m'])**2).mean())
    bias = (df_val['sim_runup_m'] - df_val['runup_m']).mean()
    r_sq = np.corrcoef(df_val['runup_m'], df_val['sim_runup_m'])[0,1]**2
    print(f"Métricas de validación (n={len(df_val)})")
    print(f"  RMSE: {rmse:.2f} m  |  Sesgo: {bias:+.2f} m  |  R²: {r_sq:.3f}")

print("\nComparación por estación:")
print(df_obs[['estacion', 'runup_m', 'sim_runup_m']].to_string(index=False))

# Gráfica de validación: observado vs. simulado
fig, axes = plt.subplots(1, 2, figsize=(13, 5))

# Panel izquierdo: scatter observado vs. simulado
ax = axes[0]
v_max = max(df_obs['runup_m'].max(), df_obs['sim_runup_m'].fillna(0).max()) * 1.15
ax.plot([0, v_max], [0, v_max], 'k--', linewidth=1, label='1:1')
ax.plot([0, v_max], [0, 2*v_max], 'gray', linewidth=0.5, linestyle=':', alpha=0.5)
ax.plot([0, v_max], [0, v_max/2], 'gray', linewidth=0.5, linestyle=':', alpha=0.5)

sc = ax.scatter(df_obs['runup_m'], df_obs['sim_runup_m'],
                c=np.abs(df_obs['lat']), cmap='viridis', s=80, zorder=5)
plt.colorbar(sc, ax=ax, label='|Latitud| (°N)')

for _, row in df_obs.iterrows():
    ax.annotate(row['estacion'], (row['runup_m'], row['sim_runup_m'] or 0),
                fontsize=7, textcoords='offset points', xytext=(4, 3))

ax.set_xlim(0, v_max)
ax.set_ylim(0, v_max)
ax.set_xlabel('Run-up observado (m) — Herd et al. 1981')
ax.set_ylabel('Run-up simulado (m)')
ax.set_title('Validación del modelo\nTsunami Tumaco 1979')
ax.grid(True, alpha=0.3)

if len(df_val) > 0:
    ax.text(0.05, 0.92, f'RMSE = {rmse:.2f} m\nR² = {r_sq:.2f}\nSesgo = {bias:+.2f} m',
            transform=ax.transAxes, fontsize=9, verticalalignment='top',
            bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.8))

# Panel derecho: barras comparativas
ax2 = axes[1]
x_pos = np.arange(len(df_obs))
width = 0.38
ax2.bar(x_pos - width/2, df_obs['runup_m'], width,
        label='Observado (Herd et al. 1981)', color='steelblue', alpha=0.85)
ax2.bar(x_pos + width/2, df_obs['sim_runup_m'].fillna(0), width,
        label='Simulado (N-wave Mw 8.2)', color='firebrick', alpha=0.85)
ax2.set_xticks(x_pos)
ax2.set_xticklabels(df_obs['estacion'], rotation=35, ha='right', fontsize=8)
ax2.set_ylabel('Run-up (m)')
ax2.set_title('Comparación run-up por estación')
ax2.legend(fontsize=9)
ax2.grid(True, axis='y', alpha=0.3)

plt.tight_layout()
plt.savefig(WORK_DIR / 'validation.png', dpi=150, bbox_inches='tight')
plt.show()
print(f"Gráfica de validación guardada en: {WORK_DIR}/validation.png")

from matplotlib.animation import FuncAnimation

fig, ax = plt.subplots(figsize=(10, 7))

# X_ct, Y_ct definidos en cell-6 (coordenadas geográficas en km para ax.contour)
from matplotlib.colors import LightSource
ls = LightSource(azdeg=315, altdeg=35)
dem_vis_hs = np.where(dem > 0, dem, 0).astype(float)
hillshade = ls.hillshade(dem_vis_hs, vert_exag=0.001)
ax.imshow(hillshade, cmap='gray', alpha=0.5, extent=ext)

ocean_mask_anim = np.ma.masked_where(dem >= 0, dem)
ax.imshow(ocean_mask_anim, cmap='Blues', vmin=-5000, vmax=0, alpha=0.4, extent=ext)

ax.contour(X_ct, Y_ct, dem, levels=[0], colors='k', linewidths=0.8)

epi_x_utm, epi_y_utm = transformer.transform(-79.359, 1.598)
ax.plot(epi_x_utm/1000, epi_y_utm/1000, 'r*', markersize=14, zorder=10,
        label='Epicentro 1979')

ANIM_CLIM = 5.0
t_sorted = sorted(heightmaps.keys())

def _wave_eta(hmap):
    return hmap + np.minimum(dem, 0)

first_eta = _wave_eta(heightmaps[t_sorted[0]])
water_layer = ax.imshow(
    np.where(np.abs(first_eta) > 0.1, first_eta, np.nan),
    cmap='RdBu_r', vmin=-ANIM_CLIM, vmax=ANIM_CLIM, alpha=0.75, extent=ext
)
plt.colorbar(water_layer, ax=ax, label='η (m sobre NMM)', shrink=0.7)

title_obj = ax.set_title('t = 0 s (0 min)', fontsize=12)
ax.set_xlabel('Easting UTM 18N (km)')
ax.set_ylabel('Northing UTM 18N (km)')
ax.legend(fontsize=9, loc='upper right')

def update_frame(frame_idx):
    t = t_sorted[frame_idx]
    eta = _wave_eta(heightmaps[t])
    water_layer.set_array(np.where(np.abs(eta) > 0.1, eta, np.nan))
    title_obj.set_text(f't = {t:.0f} s ({t/60:.0f} min)\nSimulación Tsunami Tumaco 1979')
    return [water_layer, title_obj]

anim = FuncAnimation(fig, update_frame, frames=len(t_sorted), interval=250, blit=True)
plt.tight_layout()

gif_path = WORK_DIR / 'tsunami_propagation.gif'
anim.save(str(gif_path), writer='pillow', fps=4, dpi=100)
print(f"Animación guardada: {gif_path}")
plt.close()

# Mostrar la animación en Colab
from IPython.display import Image, display
if gif_path.exists():
    display(Image(str(gif_path)))

print("=" * 60)
print("RESUMEN DE RESULTADOS — SIMULACIÓN TSUNAMI TUMACO 1979")
print("=" * 60)
print(f"  Modelo:            LSWE NumPy / Lax-Friedrichs 2D")
print(f"  Evento:            12-dic-1979, Mw 8.2")
print(f"  Resolución DEM:    {RESOLUTION:.0f} m (GEBCO 2023 / ETOPO 2022)")
print(f"  Duración sim.:     {max(t_sorted)/60:.0f} minutos")
print(f"  Snapshots:         {len(t_sorted)}")
print()
print("  --- Resultados ---")
print(f"  η máxima (océano): {max_eta[dem < 0].max():.2f} m")
print(f"  η máxima (costera <200m): {max_eta[coast_zone_val].max():.2f} m")
if len(df_val) > 0:
    print(f"  RMSE validación:   {rmse:.2f} m")
    print(f"  R²:                {r_sq:.3f}")
    print(f"  Sesgo:             {bias:+.2f} m")
print()
print("  --- Archivos generados ---")
for fname in ['max_inundation_map.png', 'time_series.png',
              'validation.png', 'tsunami_propagation.gif']:
    fp = WORK_DIR / fname
    if fp.exists():
        print(f"  ✓ {fname} ({fp.stat().st_size/1e3:.0f} KB)")
    else:
        print(f"  ✗ {fname} (no generado)")