%% =======================================================================
%  G500 - AIRCRAFT PERFORMANCE & DESIGN (Asignacion 2 + Propulsion)
%  Script integral alineado con el informe LaTeX (unidades inglesas)
%  - Reporte unico en consola, por preguntas, con interpretacion
%  - Figura unica con 4 graficas (tiledlayout)
% =======================================================================

clear; clc; close all;

% ---------- Ventana de graficas ----------
fig = figure('Name','G500 Summary - Aircraft Performance & Propulsion',...
             'NumberTitle','off','Color','w');
tlo = tiledlayout(fig,2,2,'TileSpacing','compact','Padding','compact');

% ---------- Buffer de reporte ----------
report = string([]);
add = @(s) assignin('caller','report',[evalin('caller','report'); string(s)]);

%% =====================  ISA @ FL450 (valores como en LaTeX)  =====================
FL.Mach = 0.85;
FL.h_ft = 45000;
[FL.rho, FL.a, FL.mu] = isa_atm_english(FL.h_ft);   % capa isoterma > 36 kft
FL.V    = FL.Mach * FL.a;
FL.q    = 0.5*FL.rho*FL.V^2;

%% =====================  Datos finales del profesor (LaTeX)  =====================
FINAL.Sref = 992.8;      % ft^2
FINAL.CD0  = 0.0142;     % -
FINAL.e    = 0.612;      % -
FINAL.AR   = 7.51;       % -
FINAL.k    = 1/(pi*FINAL.e*FINAL.AR);   % 0.0693...

% =========  Datos de Propulsion y Mision (PDF)  =========
N_eng      = 2;              % # motores
W_cruise   = 62000;          % lbf  (peso en crucero del PDF)
TSFC_cr    = 0.60;           % lb/(lbf·h)
m_core     = 38.0;           % lb/s por motor (core airflow)
m_bleed_dx = 0.5;            % lb/s por motor (extra bleed de-icing)
Tmax_SL    = 15429;          % lbf por motor (a nivel del mar)

%% =====================  Bloque Aerodinamico: Polar & CL de crucero  =====================
CL_cruise  = W_cruise/(FL.q*FINAL.Sref);
CD_cruise  = FINAL.CD0 + FINAL.k*CL_cruise.^2;
D_total    = FL.q*FINAL.Sref*CD_cruise;       % = Thrust requerido total
Treq_each  = D_total/N_eng;

% Puntos de referencia de la polar
CL_star    = sqrt(FINAL.CD0/FINAL.k);         % en (L/D)max
LD_max     = 1/(2*sqrt(FINAL.CD0*FINAL.k));
CD_star    = FINAL.CD0 + FINAL.k*CL_star^2;

%% =====================  PREGUNTA 1 =====================
add("=== Pregunta 1 - Thrust requerido en crucero ===");
add(sprintf("Condicion: M=%.2f, h=%d ft  |  a=%.1f ft/s, V=%.1f ft/s, rho=%.2e slug/ft^3, q=%.1f lbf/ft^2", ...
    FL.Mach, FL.h_ft, FL.a, FL.V, FL.rho, FL.q));
add(sprintf("Datos finales: Sref=%.1f ft^2, CD0=%.4f, e=%.3f, AR=%.2f -> k=%.4f", ...
    FINAL.Sref, FINAL.CD0, FINAL.e, FINAL.AR, FINAL.k));
add(sprintf("W=%.0f lbf -> CL=%.3f,  CD=%.5f", W_cruise, CL_cruise, CD_cruise));
add(sprintf("Thrust requerido total D = %.0f lbf -> por motor (x%d) = %.0f lbf", ...
    D_total, N_eng, Treq_each));
add("Interpretacion: el empuje requerido en vuelo nivelado es el arrastre total.");
add("Depende de q, S y de la polar CD = CD0 + k*CL^2. Con este set coincide con el PDF.");

%% =====================  PREGUNTA 2 =====================
FF_each_hr = TSFC_cr * Treq_each;
FF_tot_hr  = N_eng    * FF_each_hr;

add(" ");
add("=== Pregunta 2 - Fuel flow con TSFC=0.60 lb/(lbf*h) ===");
add(sprintf("Fuel flow por motor = %.0f lb/h  |  total = %.0f lb/h", ...
    FF_each_hr, FF_tot_hr));
add("Interpretacion: el consumo es proporcional al thrust solicitado.");
add("Si baja CD0 o sube e*AR, baja el thrust y el consumo horario.");

%% =====================  PREGUNTA 3 =====================
frac_loss   = m_bleed_dx/m_core;       % 0.01316...
dT_each     = frac_loss * Treq_each;   % perdida por motor
Tavail_each = Treq_each - dT_each;

add(" ");
add("=== Pregunta 3 - Penalizacion de thrust por bleed adicional (de-icing) ===");
add(sprintf("Bleed adicional = %.1f/%.1f = %.3f (%.2f%%) del core", ...
    m_bleed_dx, m_core, frac_loss, 100*frac_loss));
add(sprintf("Delta T ~ %.1f lbf por motor (thrust disponible ~ %.0f lbf si el set de motor no cambia)", ...
    dT_each, Tavail_each));
add("Interpretacion: asumiendo T proporcional a m_dot_core, el bleed reduce el empuje en el mismo porcentaje.");

%% =====================  PREGUNTA 4 =====================
TSFC_new       = TSFC_cr/(1 - frac_loss);     % aumento por perdida de core flow
FF_each_hr_new = TSFC_new * Treq_each;
FF_tot_hr_new  = N_eng * FF_each_hr_new;

add(" ");
add("=== Pregunta 4 - TSFC de crucero tras instalar el sistema de bleed ===");
add(sprintf("TSFC_old=%.3f  ->  TSFC_new=%.3f  (Delta=%.2f%%)", ...
    TSFC_cr, TSFC_new, 100*(TSFC_new/TSFC_cr - 1)));
add(sprintf("Fuel flow por motor para sostener el mismo thrust: %.0f -> %.0f lb/h", ...
    FF_each_hr, FF_each_hr_new));
add("Interpretacion: al perder parte del flujo del compresor, sube el combustible");
add("necesario para un mismo thrust efectivo, por eso TSFC crece.");

%% =====================  PREGUNTA 5 =====================
[rho_SL,~,~] = isa_atm_english(0);
sigma        = FL.rho/rho_SL;
Tmax_each_h  = Tmax_SL * sigma;
Tmax_tot_h   = N_eng   * Tmax_each_h;

add(" ");
add("=== Pregunta 5 - Thrust maximo a altitud (aprox proporcional a rho) ===");
add(sprintf("sigma = rho(h)/rho_SL = %.4f  ->  T_max por motor = %.0f lbf  |  Total = %.0f lbf", ...
    sigma, Tmax_each_h, Tmax_tot_h));
add(sprintf("Margen frente a T_req total: %.1f %%", 100*(Tmax_tot_h/D_total - 1)));
add("Interpretacion: con la hipotesis de T_max ~ rho, el tope de thrust decrece a FL450.");
add("Da el margen disponible frente al thrust requerido calculado en P1.");

%% =====================  Curvas para estudio  =====================
% Curva polar
CL = linspace(0,1.4,300);
CD = FINAL.CD0 + FINAL.k*CL.^2;

nexttile(tlo,1); % polar
plot(CD,CL,'LineWidth',1.6); grid on; hold on;
plot(CD_star,CL_star,'o','MarkerFaceColor',[0.85 0.1 0.1],'Color',[0.85 0.1 0.1]);
xlabel('C_D'); ylabel('C_L');
title(sprintf('Polar final (CD_0=%.4f, e=%.3f, AR=%.2f)',FINAL.CD0,FINAL.e,FINAL.AR));
% Etiqueta del punto de (L/D)max sin interprete TeX (evita errores con \newline)
txt = sprintf('(L/D)_max=%.1f\nC_L^*=%.3f', LD_max, CL_star);
text(CD_star, CL_star, txt, 'VerticalAlignment','bottom','FontSize',9, 'Interpreter','none');

% SFC parcial (tabla del curso, para sensibilidad)
PP.Tfrac = linspace(0.2,1,60); PP.M=FL.Mach; PP.cmax=0.50;
PP.c_ratio = (0.1./PP.Tfrac) + (0.24./PP.Tfrac.^0.8) + 0.66*(PP.Tfrac.^0.8) ...
           + (0.1*PP.M).*(1./PP.Tfrac - PP.Tfrac);
nexttile(tlo,2);
plot(PP.Tfrac,PP.c_ratio,'LineWidth',1.6); grid on;
xlabel('T/T_{max,dry}'); ylabel('c/c_{max,dry}');
title('SFC relativo vs fraccion de empuje');

% Motor de piston (para estudio)
PIST.Power_SL=200; PIST.h=linspace(0,25000,50); rho0=0.002377;
PIST.rho=zeros(size(PIST.h));
for i=1:numel(PIST.h)
    [ri,~,~]=isa_atm_english(PIST.h(i));
    PIST.rho(i)=ri;
end
PIST.sigma=PIST.rho/rho0; PIST.Power=PIST.Power_SL.*(PIST.sigma-(1-PIST.sigma)/7.55);
nexttile(tlo,3);
plot(PIST.h,PIST.Power,'LineWidth',1.6); grid on;
xlabel('Altitud [ft]'); ylabel('Potencia [hp]'); title('Piston: Potencia vs Altitud');

% (Espacio para otra grafica si deseas)

%% =====================  Reporte unico =====================
disp(newline + "================= REPORTE UNICO - G500 (Asignacion 2) ================");
disp(join(report,newline));
disp("=======================================================================");

%% =====================  Funcion ISA coherente con FL450  =====================
function [rho, a, mu] = isa_atm_english(h_ft)
% ISA estandar (unidades inglesas) coherente con FL450 (isotermia >36kft)
% OUT: rho [slug/ft^3], a [ft/s], mu [slug/(ft·s)]
    g0 = 32.174;  R = 1716.59;  gamma = 1.4;
    T0 = 518.67;  P0 = 2116.22; L = -0.00356616;  % 0–36 kft
    if h_ft <= 36089
        T = T0 + L*h_ft;
        P = P0*(T/T0)^(-g0/(L*R));
    else
        T = 389.97; % 216.65 K
        T36 = T0 + L*36089;
        P36 = P0*(T36/T0)^(-g0/(L*R));
        P = P36*exp(-g0*(h_ft-36089)/(R*T));
    end
    rho = P/(R*T);
    a   = sqrt(gamma*R*T);
    % Viscosidad (Sutherland) -> slug/(ft*s)
    T_K  = (T-459.67)/1.8 + 273.15;
    mu_SI= 1.716e-5*(T_K/273.15)^(3/2)*(273.15+110.4)/(T_K+110.4); % Pa·s
    mu   = mu_SI * 1.4881639;
end
