utils.py

import numpy as np
from scipy import stats
from scipy.stats import f as f_dist

def cargar_datos_montecarlo(ruta='montecarlo.dat'):
    """
    Carga el .dat de la campaña Monte Carlo y calcula las variables derivadas
    (margen estático, radio de impacto, Re y Ma en el instante del cruce).

    Retorna un dict con todos los arrays, para desempaquetar en el script principal.
    """
    data = np.loadtxt(ruta, comments='#')

    corrida    = data[:,0]
    x_imp      = data[:,1]
    y_imp      = data[:,2]
    crash_bw   = data[:,3]
    h_crash    = data[:,4]
    t_crash    = data[:,5]
    X_cg       = data[:,6]
    X_cp       = data[:,7]
    theta0     = data[:,8]
    alpha0     = data[:,9]
    M0         = data[:,10]
    Vf_crash   = data[:,11]
    WE_crash   = data[:,12]
    WN_crash   = data[:,13]
    viento_RMS = data[:,14]
    alpha_max  = data[:,15]
    Jx         = data[:,16]
    Jy         = data[:,17]
    Jz         = data[:,18]

    dc        = 0.08103
    margen    = (X_cp - X_cg) / dc
    radio_imp = np.sqrt(x_imp**2 + y_imp**2)

    T_crash   = np.maximum(288.15 - 0.0065 * h_crash, 216.65)
    rho_crash = 1.225 * (T_crash / 288.15)**4.256
    mu_crash  = 1.716e-5 * (T_crash / 273.15)**0.666
    a_crash   = np.sqrt(1.4 * 287.05 * T_crash)
    L_ref     = 0.08103

    Re_crash = rho_crash * Vf_crash * L_ref / mu_crash
    Ma_crash = Vf_crash / a_crash

    return dict(
        corrida=corrida, x_imp=x_imp, y_imp=y_imp, crash_bw=crash_bw,
        h_crash=h_crash, t_crash=t_crash, X_cg=X_cg, X_cp=X_cp,
        theta0=theta0, alpha0=alpha0, M0=M0, Vf_crash=Vf_crash,
        WE_crash=WE_crash, WN_crash=WN_crash, viento_RMS=viento_RMS,
        alpha_max=alpha_max, Jx=Jx, Jy=Jy, Jz=Jz,
        margen=margen, radio_imp=radio_imp,
        Re_crash=Re_crash, Ma_crash=Ma_crash,
    )

def elipse_confianza(X, Y, nivel=0.95, n_puntos=300):
    """
    Elipse de confianza al nivel dado usando distribución de Hotelling T².
    Estadísticamente correcta para dos dimensiones — no asume gaussianidad
    marginal, usa la distribución F exacta para el umbral multivariado.

    Retorna: (centro, semi_eje_mayor, semi_eje_menor, puntos_x, puntos_y)
    """
    data2d = np.column_stack([X, Y])
    n, p   = len(data2d), 2
    mu     = data2d.mean(axis=0)
    S      = np.cov(data2d.T)

    # Umbral de Hotelling T² con corrección exacta a distribución F
    F_crit = f_dist.ppf(nivel, p, n - p)
    c2     = (p * (n - 1)) / (n - p) * F_crit

    # Descomposición espectral: eigenvalores dan los semi-ejes
    vals, vecs = np.linalg.eigh(S)
    a_eje  = np.sqrt(vals[1] * c2)   # semi-eje mayor
    b_eje  = np.sqrt(vals[0] * c2)   # semi-eje menor

    # Parametrización de la elipse en el espacio original (sin normalizar)
    theta = np.linspace(0, 2*np.pi, n_puntos)
    pts   = (mu[:, None]
             + a_eje * np.cos(theta) * vecs[:, 1][:, None]
             + b_eje * np.sin(theta) * vecs[:, 0][:, None])

    return mu, a_eje, b_eje, pts[0], pts[1]

def extension_por_eje(ex, ey):
    """
    Extensión de la elipse proyectada sobre cada eje original (bounding box / 2),
    en las unidades propias de cada variable.
    A diferencia de a_eje/b_eje (semi-ejes del eigenespacio, que mezclan
    unidades distintas cuando las variables tienen escalas muy dispares),
    esto SÍ es interpretable directamente: "la elipse llega hasta ±X en
    esta variable, ±Y en la otra".
    """
    ext_x = (ex.max() - ex.min()) / 2.0
    ext_y = (ey.max() - ey.min()) / 2.0
    return ext_x, ext_y


def regresion_estandarizada(X_raw, y_raw, nombre_y, nombres_x, verbose=True):
    """
    Regresión lineal múltiple con inputs y output estandarizados.
    Los coeficientes β* son directamente comparables entre sí —
    el mayor en valor absoluto es el parámetro más influyente.

    verbose=False permite reutilizar el cálculo (p.ej. para graficar)
    sin repetir la impresión por consola.
    """
    # Estandarizar
    X_std = (X_raw - X_raw.mean(axis=0)) / X_raw.std(axis=0)
    y_std = (y_raw - y_raw.mean()) / y_raw.std()

    # Ajuste por mínimos cuadrados
    betas, _, _, _ = np.linalg.lstsq(
        np.column_stack([np.ones(len(y_std)), X_std]),
        y_std, rcond=None)
    betas = betas[1:]   # descartar intercepto

    # R² del ajuste
    y_pred = X_std @ betas
    ss_res = np.sum((y_std - y_pred)**2)
    ss_tot = np.sum((y_std - y_std.mean())**2)
    R2     = 1 - ss_res / ss_tot

    # Correlación de Spearman (robusta a no-linealidad) para comparar
    r_spearman = np.array([
        stats.spearmanr(X_raw[:, i], y_raw).correlation
        for i in range(X_raw.shape[1])])

    if verbose:
        print(f"\n  {nombre_y}  (R² = {R2:.3f})")
        print(f"  {'Parámetro':>10}  {'β* (OLS)':>10}  {'r Spearman':>12}  {'domina':>8}")
        orden = np.argsort(np.abs(betas))[::-1]
        for k in orden:
            marca = '◄' if k == orden[0] else ''
            print(f"  {nombres_x[k]:>10}  {betas[k]:>+10.4f}  "
                  f"{r_spearman[k]:>+12.4f}  {marca}")
    return betas, r_spearman, R2


def parametro_dominante(resultado, nombres_x):
    """Nombre del parámetro con mayor |β*| en un resultado de regresion_estandarizada."""
    betas = resultado[0]
    return nombres_x[np.argmax(np.abs(betas))]


def presion_dinamica_isa(h_geopotencial_m):
    """
    Calcula densidad del aire (kg/m^3) via Atmosfera Estandar Internacional (ISA),
    valida para 0-11000 m (troposfera) y 11000-20000 m (estratosfera baja).
    h_geopotencial_m: altitud geopotencial en metros (escalar o array numpy).
    Retorna: rho (misma forma que la entrada).
    """
    h = np.asarray(h_geopotencial_m, dtype=float)
    rho = np.zeros_like(h)

    # Constantes ISA
    T0    = 288.15      # K, temperatura a nivel del mar
    P0    = 101325.0    # Pa, presion a nivel del mar
    L     = 0.0065       # K/m, gradiente termico troposfera
    R_gas = 287.05       # J/(kg*K), constante especifica del aire
    g0    = 9.80665      # m/s^2

    # --- Troposfera (0 - 11000 m) ---
    mask_tropo = h <= 11000.0
    T_tropo = T0 - L * h[mask_tropo]
    P_tropo = P0 * (T_tropo / T0) ** (g0 / (R_gas * L))
    rho[mask_tropo] = P_tropo / (R_gas * T_tropo)

    # --- Estratosfera baja (11000 - 20000 m), temperatura constante ---
    mask_estrato = h > 11000.0
    if np.any(mask_estrato):
        T_11 = T0 - L * 11000.0
        P_11 = P0 * (T_11 / T0) ** (g0 / (R_gas * L))
        h_rel = h[mask_estrato] - 11000.0
        P_estrato = P_11 * np.exp(-g0 * h_rel / (R_gas * T_11))
        rho[mask_estrato] = P_estrato / (R_gas * T_11)

    return rho

def clasificar_poblacion_crash(h_crash, Vf_crash, ht_rampa=2000.0, q_umbral=None):
    """
    Reconstruye la presion dinamica en el instante de cruce del limite de Barrowman
    y clasifica cada corrida en 'divergencia_aero' (alta q, mecanismo fisico real)
    o 'geometrico_apogeo' (Vf bajo cerca de apogeo, artefacto de alpha=atan2 con Vf->0).

    h_crash: array, altitud geodesica real en el instante de cruce (m, sobre el elipsoide)
    Vf_crash: array, velocidad total en el instante de cruce (m/s)
    ht_rampa: altitud de la rampa de lanzamiento (m), para convertir a altitud relativa/AGL
    q_umbral: umbral de presion dinamica (Pa) para separar poblaciones.
              Si None, se usa la mediana de q_dyn como umbral (punto de partida,
              ajustar tras inspeccionar el histograma bimodal si existe).

    Retorna: dict con 'q_dyn' (array, Pa), 'poblacion' (array de strings),
             'q_umbral' (float, el umbral usado)
    """
    h_agl = h_crash - ht_rampa          # altitud sobre la rampa, no sobre el elipsoide
    rho = presion_dinamica_isa(h_crash) # ISA usa altitud geopotencial real, no AGL
    q_dyn = 0.5 * rho * Vf_crash**2      # Pa

    if q_umbral is None:
        q_umbral = np.median(q_dyn)

    poblacion = np.where(q_dyn >= q_umbral, 'divergencia_aero', 'geometrico_apogeo')

    return {'q_dyn': q_dyn, 'poblacion': poblacion, 'q_umbral': q_umbral, 'h_agl': h_agl}