function atm = isa_std_atmos(h_ft, varargin)
% ISA_STD_ATMOS  Modelo simple de atmósfera estándar (ISA) centrado en 30 kft.
%
% USO:
%   atm = isa_std_atmos();                      % por defecto h_ft = 30000 [ft]
%   atm = isa_std_atmos(h_ft);                  % cálculo en h_ft [ft]
%   atm = isa_std_atmos(h_ft,'verbose',true);   % imprime resumen explicado
%
% ENTRADAS:
%   h_ft      (opcional) Altitud presión [ft]. Por defecto 30000 [ft].
%   'verbose' (opcional) true/false para imprimir resumen (false por defecto).
%
% SALIDA (struct atm):
%   atm.T_K       Temperatura [K]
%   atm.T_R       Temperatura [R] (Rankine)   T_R = (9/5) * T_K
%   atm.a_fts     Velocidad del sonido [ft/s] a = sqrt(gamma * R * T_R)
%   atm.rho       Densidad del aire [slug/ft^3]
%   atm.sigma     Relación de densidad [-] = rho / rho0
%   atm.rho0      Densidad a nivel del mar [slug/ft^3] (referencia)
%   atm.h_ft      Altitud usada [ft]
%   atm.model     Descripción breve del modelo usado
%
% CONSTANTES:
%   gamma = 1.4;  R = 1716 [ft*lbf/(slug*R)]; rho0 = 2.377e-3 [slug/ft^3]
%
% NOTAS:
% - A 30 000 ft usa un punto ISA validado (T=228.7 K; rho≈8.89e-4 slug/ft^3).
% - Fuera de ~30 kft aplica una aproximación simple (tropo/estrato) suficiente
%   para cálculos de rendimiento preliminares en unidades inglesas.
%
% Autor: Aircraft Performance & DSG

% atm = isa_std_atmos();                 % 30 kft por defecto, sin imprimir
% atm = isa_std_atmos(30000,'verbose',true);   % imprime resumen explicado
% atm = isa_std_atmos(35000,'verbose',true);   % otra altitud



    % --------------------------
    % Entradas y opciones
    % --------------------------
    if nargin < 1 || isempty(h_ft)
        h_ft = 30000;                   % valor por defecto
    end
    p = inputParser;
    p.addParameter('verbose', false, @(x)islogical(x) || isnumeric(x));
    p.parse(varargin{:});
    verbose = logical(p.Results.verbose);

    % --------------------------
    % Constantes ISA (imperial)
    % --------------------------
    gamma = 1.4;                 % razón de calores específicos [-]
    R     = 1716;                % constante de gas [ft*lbf/(slug*R)]
    rho0  = 2.377e-3;            % densidad al nivel del mar [slug/ft^3]
    T0_R  = 518.67;              % temperatura SL [R] (59°F)

    % --------------------------
    % Cálculo T, a, rho
    % --------------------------
    if abs(h_ft - 30000) < 1e-6 || (h_ft > 25000 && h_ft < 35000)
        % Punto "ancla" usado en el curso (ISA @ 30 kft):
        T_K  = 228.7;                   % [K]
        T_R  = T_K * 9/5;               % [R]
        rho  = 8.89e-4;                 % [slug/ft^3]
        mdl  = 'ISA @30kft (punto validado)';
    else
        % Aproximación simple fuera de 30 kft
        if h_ft <= 36089                % troposfera (gradiente aproximado)
            lapse_R_per_ft = -0.00356616;   % [R/ft] ≈ -6.5 K/km
            T_R = T0_R + lapse_R_per_ft * h_ft;
        else                            % estratosfera baja (isotermo aprox.)
            T_R = 389.97;               % ≈ 216.65 K en Rankine
        end
        T_K  = T_R * 5/9;               % [K]
        rho  = max(7.0e-4, rho0 * (T_R/T0_R)^4.256);  % muy aprox con cota
        mdl  = 'ISA simplificada (tropo/estrato)';
    end

    a_fts = sqrt(gamma * R * T_R);      % velocidad del sonido [ft/s]
    sigma = rho / rho0;                 % densidad adimensional

    % --------------------------
    % Empaquetar salida
    % --------------------------
    atm = struct();
    atm.T_K    = T_K;                   % [K]
    atm.T_R    = T_R;                   % [R]
    atm.a_fts  = a_fts;                 % [ft/s]
    atm.rho    = rho;                   % [slug/ft^3]
    atm.sigma  = sigma;                 % [-]
    atm.rho0   = rho0;                  % [slug/ft^3]
    atm.h_ft   = h_ft;                  % [ft]
    atm.model  = mdl;

    % --------------------------
    % Modo "verbose": imprimir resumen explicado
    % --------------------------
    if verbose
        fprintf('\n=== Atmósfera ISA (h = %.0f ft) ===\n', h_ft);
        explain('T_K',   atm.T_K,   'K',        'Temperatura absoluta (Kelvin)');
        explain('T_R',   atm.T_R,   'R',        'Temperatura absoluta (Rankine)');
        explain('a',     atm.a_fts, 'ft/s',     'Velocidad del sonido, a = sqrt(γ R T)');
        explain('rho',   atm.rho,   'slug/ft^3','Densidad del aire');
        explain('sigma', atm.sigma, '-',        'Relación de densidad σ = ρ/ρ₀');
        explain('rho0',  atm.rho0,  'slug/ft^3','Densidad a nivel del mar (referencia)');
        fprintf('Modelo usado: %s\n\n', atm.model);
    end
end

% -------- función local para impresión formateada --------
function explain(symb, val, unit, txt)
    if isnumeric(val), sval = num2str(val,'%.5g'); else, sval = char(val); end
    fprintf('  %-6s = %-12s [%s]  %s\n', symb, sval, unit, txt);
end



% ============================================================
%   AIRCRAFT PERFORMANCE & DESIGN – PIPELINE ASIGNACIONES 1–4
% ============================================================
%
% 🔹 Propósito
% Este conjunto de funciones automatiza los cálculos de 
% performance (Asignaciones 1–4), permitiendo trabajar tanto
% con CD0 base constante como con la curva ΔCD0(M) digitalizada.
% Los módulos están diseñados para ejecutarse en cadena y 
% generar tablas, gráficas y resúmenes con explicación.
%
% 🔹 Flujo de ejecución principal (pipeline)
%
% run_airperf.m
%   └─> isa_std_atmos     → atm
%   └─> base_quantities   → base
%   └─> build_CD0_total   → CD0tot
%   └─> build_drag_table  → T
%   └─> check_mach_cap    → cap
%   └─> breguet_ER        → breg
%   └─> describe_airperf_outputs(atm, base, CD0tot, T, cap, breg)
%
% 🔹 Descripción de cada módulo
%
% isa_std_atmos(h_ft)
%   ENTRADA: altitud [ft]
%   SALIDA : atm struct con:
%       T_K, T_R, a_fts, rho, sigma
%   → Modelo de atmósfera estándar (ISA).
%
% base_quantities(atm, inputs)
%   ENTRADA: atm, struct inputs con {S, W, AR, e, TSFC, T_SL_tot}
%   SALIDA : base struct con:
%       q0, q0S, k, TSFC corregido
%   → Calcula magnitudes base aerodinámicas.
%
% build_CD0_total(M_grid, CD0_base, dCD0)
%   ENTRADA: malla Mach, CD0_base, curva ΔCD0(M)
%   SALIDA : CD0tot vector
%   → Construye la curva de arrastre parasitario.
%     (dos caminos: CD0_base constante o digitalización ΔCD0).
%
% build_drag_table(M_grid, CD0tot, base, inputs)
%   ENTRADA: Mach grid, CD0tot, base, inputs
%   SALIDA : tabla T con columnas:
%       M, C_D0, C_D, C_L, D_par, D_ind, D_tot
%   → Calcula arrastres parasitario, inducido y total.
%
% check_mach_cap(M_cap, CD0tot, base, inputs)
%   ENTRADA: Mach operativo
%   SALIDA : cap struct con:
%       C_D_cap, D_cap, T_avail
%   → Compara empuje disponible vs arrastre.
%     Determina si gobierna Mach límite o empuje.
%
% breguet_ER(inputs, base, atm)
%   ENTRADA: TSFC, (L/D), Wi/Wf, Mach operativo
%   SALIDA : breg struct con:
%       E (endurance), R (range)
%   → Implementa ecuaciones de Breguet.
%
% describe_airperf_outputs(atm, base, CD0tot, T, cap, breg)
%   ENTRADA: todos los structs anteriores
%   SALIDA : impresión en consola con unidades y explicación
%   → Resume y documenta resultados.
%
% 🔹 Notas de uso
% - Siempre arrancar desde run_airperf.m
% - Cada módulo puede probarse aislado llamando con sus argumentos.
% - Diseñado para Matlab Online con entrada/salida clara y
%   documentada para reportes técnicos en LaTeX.
%
% ============================================================
