%% =======================================================================
%  G500 — AIRCRAFT PERFORMANCE & DESIGN (Asignación 2 + Propulsión)
%  Script integral, corregido y documentado (unidades inglesas)
% =======================================================================

clear; clc; close all;

%% =======================================================================
%% 1) AERODINÁMICA — CD0 por COMPONENT BUILD-UP (crucero, ISA)
%% =======================================================================

% ---- Condición de crucero ----
CRU.Mach  = 0.85;                    % Número de Mach
CRU.h_ft  = 45000;                   % Altitud [ft]
[CRU.rho, CRU.a, CRU.mu] = isa_atm_english(CRU.h_ft);
CRU.V     = CRU.Mach * CRU.a;        % Velocidad verdadera [ft/s]

% ---- Geometría ala (Sref de referencia) ----
W.b    = 86.333;     % envergadura [ft]
W.cR   = 19.0;       % cuerda raíz [ft]
W.cT   = 4.4;        % cuerda punta [ft]
W.tc   = 0.14;       % espesor relativo t/c
W.x_tc = 0.37;       % posición de t/c máx (x/c)
W.LE   = 30;         % flecha LE [deg]
W.cRexp= 17.6;       % cuerda raíz expuesta [ft]
W.dcut = 7.3;        % descuento por fuselaje [ft]

% ---- Estabilizador horizontal ----
H.b    = 32.7;  H.cR = 9.6;   H.cT = 4.2;
H.tc   = 0.10;  H.LE = 33;    H.d  = 1.0;

% ---- Estabilizador vertical ----
Vv.b   = 12.0;  Vv.cR=14.1;   Vv.cT = 9.0;
Vv.tc  = 0.10;  Vv.LE=40;

% ---- Fuselaje ----
F.l    = 78.8;  F.d  = 7.3;
F.S_top= 483.4; F.S_side = 528.4;

% ---- Góndolas ----
Nac.D  = 4.6;    % diámetro medio [ft]
Nac.l  = 15.25;  % longitud [ft]

% ---- Constantes / factores ----
CF.ksmooth = 2.08e-5;   % rugosidad (ft)
CF.C_lam   = 1.3828;    % Cf_lam = C/sqrt(Re)
CF.Q_wing  = 1.00;      % interferencia ala
CF.Q_tail  = 1.03;      % interferencia colas
CF.Q_nac   = 1.50;      % interferencia góndolas
CF.LP      = 1.10;      % +10% fugas/protuberancias

% ---- Áreas, AR y cuerdas medias ----
W.lambda = W.cT/W.cR;
W.Sref   = W.b/2*(W.cR+W.cT);                 % referencia
W.AR     = W.b^2/W.Sref;
W.cbar   = (2/3)*W.cR*(1+W.lambda+W.lambda^2)/(1+W.lambda);

W.Sexp   = (W.b - W.dcut)/2 * (W.cRexp + W.cT);
W.Swet   = W.Sexp*(1.977 + 0.52*W.tc);

H.lambda = H.cT/H.cR;
H.S      = (H.b - H.d)/2*(H.cR+H.cT);
H.AR     = H.b^2/H.S;
H.cbar   = (2/3)*H.cR*(1+H.lambda+H.lambda^2)/(1+H.lambda);
H.Swet   = H.S*(1.977 + 0.52*H.tc);

Vv.lambda= Vv.cT/Vv.cR;
Vv.S     = Vv.cR*Vv.b*(1+Vv.lambda)/2;
Vv.AR    = Vv.b^2/Vv.S;
Vv.cbar  = (2/3)*Vv.cR*(1+Vv.lambda+Vv.lambda^2)/(1+Vv.lambda);
Vv.Swet  = 2*Vv.S*(1.977 + 0.52*Vv.tc);

F.f      = F.l/F.d;
F.Swet   = 3.4*((F.S_top + F.S_side)/2);

Nac.Swet = 4*pi*Nac.D*Nac.l;

% ---- Reynolds y Re_cutoff ----
Re.W  = CRU.rho*CRU.V*W.cbar  /CRU.mu;  Re.Wc  = 38.21*(W.cbar /CF.ksmooth)^1.053;
Re.H  = CRU.rho*CRU.V*H.cbar  /CRU.mu;  Re.Hc  = 38.21*(H.cbar /CF.ksmooth)^1.053;
Re.Vv = CRU.rho*CRU.V*Vv.cbar /CRU.mu;  Re.Vvc = 38.21*(Vv.cbar/CF.ksmooth)^1.053;
Re.F  = CRU.rho*CRU.V*F.l     /CRU.mu;  Re.Fc  = 38.21*(F.l    /CF.ksmooth)^1.053;
Re.N  = CRU.rho*CRU.V*Nac.l   /CRU.mu;  Re.Nc  = 38.21*(Nac.l  /CF.ksmooth)^1.053;

% ---- Cf y FF ----
% Ala (mezcla 65/35)
Cf.W_lam  = CF.C_lam/sqrt(Re.W);
Cf.W_turb = 0.455/(log10(min(Re.W,Re.Wc))^2.58*(1+0.144*CRU.Mach^2)^0.65);
Cf.W      = 0.65*Cf.W_turb + 0.35*Cf.W_lam;
lambda_max = atan( tan(deg2rad(W.LE))*(4*W.x_tc*(1-W.lambda)/(W.AR*(1+W.lambda))) );
FF.W      = (1 + 0.6/W.x_tc*W.tc + 100*W.tc^4) * 1.34*CRU.Mach^0.18 * cos(lambda_max)^0.28;

% H-tail
Cf.H_lam  = CF.C_lam/sqrt(Re.H);
Cf.H_turb = 0.455/(log10(min(Re.H,Re.Hc))^2.58*(1+0.144*CRU.Mach^2)^0.65);
Cf.H      = 0.65*Cf.H_turb + 0.35*Cf.H_lam;
Lambda40H = atan( tand(H.LE) - 4*0.4*(1-H.lambda)/(H.AR*(1+H.lambda)) );
FF.H      = (1+0.6/0.40*H.tc+100*H.tc^4)*1.34*CRU.Mach^0.18*cos(Lambda40H)^0.28;

% V-tail
Cf.V_lam  = CF.C_lam/sqrt(Re.Vv);
Cf.V_turb = 0.455/(log10(min(Re.Vv,Re.Vvc))^2.58*(1+0.144*CRU.Mach^2)^0.65);
Cf.V      = 0.65*Cf.V_turb + 0.35*Cf.V_lam;
Lambda40V = atan( tand(Vv.LE) - 4*0.4*(1-Vv.lambda)/(2*Vv.AR*(1+Vv.lambda)) );
FF.V      = (1+0.6/0.40*Vv.tc+100*Vv.tc^4)*1.34*CRU.Mach^0.18*cos(Lambda40V)^0.28;

% Fuselaje (100% turbulento)
Cf.F_turb = 0.455/(log10(min(Re.F,Re.Fc))^2.58*(1+0.144*CRU.Mach^2)^0.65);
FF.F      = 1 + 60/F.f^3 + F.f/400;

% Góndolas (90% turb + 10% lam)
Cf.N_lam  = CF.C_lam/sqrt(Re.N);
Cf.N_turb = 0.455/(log10(min(Re.N,Re.Nc))^2.58*(1+0.144*CRU.Mach^2)^0.65);
Cf.N      = 0.9*Cf.N_turb + 0.1*Cf.N_lam;
FF.N      = 1 + 0.35/(Nac.l/Nac.D);

% ---- ΔCD0 por componente ----
dCD.W  = Cf.W*FF.W*CF.Q_wing*(W.Swet/W.Sref);
dCD.H  = Cf.H*FF.H*CF.Q_tail*(H.Swet/W.Sref);
dCD.V  = Cf.V*FF.V*CF.Q_tail*(Vv.Swet/W.Sref);
dCD.F  = Cf.F_turb*FF.F*(F.Swet/W.Sref);
dCD.N  = Cf.N*FF.N*CF.Q_nac*(Nac.Swet/W.Sref);

% Upsweep fuselaje (de la lámina)
UP.u_deg = 3.4; UP.Amax = 48.0;               % ft^2
dCD.UP   = (3.83*(UP.u_deg^2.5)*UP.Amax)/W.Sref;

% ---- Suma y LP ----
CD0_build    = dCD.W + dCD.H + dCD.V + dCD.F + dCD.N + dCD.UP;
CD0_build_LP = CD0_build * CF.LP;

% ---- Cfe (método equivalente) ----
Swet_total = W.Swet + H.Swet + Vv.Swet + F.Swet + Nac.Swet;
Cfe_ref    = 0.0026;
CD0_equiv  = Cfe_ref * Swet_total / W.Sref;

% ---- Reporte ----
fprintf('\n=== CD0 por Component Build-up @ Crucero ===\n');
fprintf('Sref=%.1f ft^2 | AR=%.2f | cbar=%.2f ft | V=%.1f ft/s\n', W.Sref, W.AR, W.cbar, CRU.V);
fprintf('Re_w=%.3e  Re_H=%.3e  Re_V=%.3e\n', Re.W, Re.H, Re.Vv);
fprintf('-- ΔCD: Wing=%.5f  H=%.5f  V=%.5f  Fusel=%.5f  Nac=%.5f  Upsweep=%.5f\n', ...
    dCD.W, dCD.H, dCD.V, dCD.F, dCD.N, dCD.UP);
fprintf('CD0 (sin LP)=%.5f  |  CD0 (LP=%.2f)=%.5f\n', CD0_build, CF.LP, CD0_build_LP);
fprintf('Swet_total=%.0f ft^2  |  Cfe=%.4f  =>  CD0_equiv=%.4f\n', Swet_total, Cfe_ref, CD0_equiv);

figure; bar(categorical({'Wing','H','V','Fusel','Nac','Upsweep'}), ...
            [dCD.W dCD.H dCD.V dCD.F dCD.N dCD.UP],'LineWidth',1);
ylabel('\Delta C_D'); title('CD_0 por componente @ crucero'); grid on;

%% =======================================================================
%% 2) POLAR DE ARRASTRE — VALORES FINALES OFICIALES
%% =======================================================================

FINAL.Sref = 992.8;      % [ft^2]
FINAL.CD0  = 0.0142;     % [-]
FINAL.e    = 0.612;      % [-]
FINAL.AR   = 7.38;       % [-]

FINAL.k     = 1/(pi*FINAL.e*FINAL.AR);       % CD = CD0 + k CL^2
FINAL.Seq   = FINAL.CD0 * FINAL.Sref;        % área equivalente
FINAL.CLs   = sqrt(FINAL.CD0/FINAL.k);       % CL* en (L/D)max
FINAL.LDmax = 1/(2*sqrt(FINAL.CD0*FINAL.k)); % (L/D)max

fprintf('\n=== Polar final (oficial) ===\n');
fprintf('k=%.5f | Seq=%.2f ft^2 | CL*=%.3f | (L/D)_max=%.2f\n', ...
    FINAL.k, FINAL.Seq, FINAL.CLs, FINAL.LDmax);

CL = linspace(0,1.4,300);
CD = FINAL.CD0 + FINAL.k*CL.^2;

figure;
plot(CD,CL,'LineWidth',1.8); grid on; hold on;
xlabel('C_D'); ylabel('C_L');
title(sprintf('G500 Polar (CD_0=%.4f, e=%.3f, AR=%.2f)', FINAL.CD0, FINAL.e, FINAL.AR));

CDs = FINAL.CD0 + FINAL.k*FINAL.CLs^2;
plot(CDs, FINAL.CLs, 'ro', 'MarkerFaceColor','r');
% >>> Corrección: usar TeX con \newline para evitar el warning
txt = sprintf('(L/D)_{max}=%.1f\\newline C_L^*=%.3f', FINAL.LDmax, FINAL.CLs);
text(CDs, FINAL.CLs, txt, 'Interpreter','tex', ...
     'VerticalAlignment','bottom','FontSize',10);

%% =======================================================================
%% 3) PROPULSIÓN — CHORRO (F, Pt, Pexp, eta_PE)
%% =======================================================================

JET.rho   = 0.002377;          % [slug/ft^3]
JET.S_jet = 3.0;                % [ft^2]
JET.Vjet  = 1500;               % [ft/s]
JET.V0    = 800;                % [ft/s]

JET.mdot = JET.rho*JET.Vjet*JET.S_jet;
JET.F    = JET.mdot*(JET.Vjet - JET.V0);
JET.Pt   = JET.F * JET.V0;
JET.Pexp = 0.5*JET.mdot*(JET.Vjet^2 - JET.V0^2);
JET.etaPE     = JET.Pt/JET.Pexp;
JET.etaPE_alt = 2/(JET.Vjet/JET.V0 + 1);

fprintf('\n=== Jet — teoría de momento ===\n');
fprintf('mdot=%.2f slug/s | F=%.1f lbf | Pt=%.0f hp | Pexp=%.0f hp | eta_PE=%.3f\n', ...
    JET.mdot, JET.F, JET.Pt/550, JET.Pexp/550, JET.etaPE);

%% =======================================================================
%% 4) PROPULSIÓN — PARTE POWER (SFC vs T/Tmax, M)
%% =======================================================================

PP.Tfrac = linspace(0.2,1,60);
PP.M     = 0.85;
PP.cmax  = 0.50;  % [lbm/(lbf·hr)] ejemplo

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);
PP.c_partial = PP.c_ratio * PP.cmax;

figure; plot(PP.Tfrac, PP.c_ratio,'LineWidth',2); grid on;
xlabel('T/T_{max,dry}'); ylabel('c/c_{max,dry}');
title(sprintf('SFC relativo vs fracción de empuje (M=%.2f)',PP.M));

%% =======================================================================
%% 5) PROPULSIÓN — MOTOR DE PISTÓN (Power vs altitud)
%% =======================================================================

PIST.Power_SL = 200;                 % [hp]
PIST.h = linspace(0,25000,50);       % [ft]
rho0 = 0.002377;

PIST.rho = zeros(size(PIST.h));
for i=1:numel(PIST.h)
    [rho_i,~,~] = isa_atm_english(PIST.h(i));
    PIST.rho(i) = rho_i;
end
PIST.sigma = PIST.rho / rho0;
PIST.Power = PIST.Power_SL .* ( PIST.sigma - (1 - PIST.sigma)/7.55 );

figure; plot(PIST.h, PIST.Power,'LineWidth',2); grid on;
xlabel('Altitud [ft]'); ylabel('Potencia disponible [hp]');
title('Motor de pistón: Potencia vs altitud');

%% =======================================================================
%% 6) PROPULSIÓN — HÉLICE (eficiencia, coeficientes, correcciones)
%% =======================================================================

% ---- Eficiencia básica ----
PROP.T      = 2000;             % [lbf]
PROP.V      = 300;              % [ft/s]
PROP.Ps_hp  = 1200;             % [hp]
PROP.Ps     = PROP.Ps_hp*550;   % [ft·lbf/s]
PROP.eta_p  = (PROP.T*PROP.V)/PROP.Ps;

fprintf('\n=== Hélice — eficiencia básica ===\n');
fprintf('T=%.0f lbf | V=%.0f ft/s | Ps=%.0f hp | eta_p=%.3f\n', ...
    PROP.T, PROP.V, PROP.Ps_hp, PROP.eta_p);

% ---- Coeficientes adimensionales ----
rho = 0.002377;               % SL
PROP.n   = 1800/60;           % *** 1800 rpm -> 30 rev/s (corregido) ***
PROP.D   = 8;                 % [ft]
PROP.Vad = 250;               % [ft/s]
PROP.P_hp= 550;  PROP.P = PROP.P_hp*550;
PROP.eta_assumed = 0.85;

PROP.J  = PROP.Vad/(PROP.n*PROP.D);
PROP.cP = PROP.P /(rho*PROP.n^3*PROP.D^5);
PROP.Tf = PROP.P*PROP.eta_assumed/PROP.Vad;
PROP.cT = PROP.Tf /(rho*PROP.n^2*PROP.D^4);
PROP.cS = PROP.Vad^5 * sqrt(rho/(PROP.P*PROP.n^2));
PROP.T_static = (PROP.cT/PROP.cP) * (PROP.P/(PROP.n*PROP.D));

fprintf('\n=== Hélice — coeficientes ===\n');
fprintf('J=%.3f  cP=%.3f  cT=%.3f  cS=%.3f  T_fwd=%.0f lbf  T_static=%.0f lbf\n', ...
    PROP.J, PROP.cP, PROP.cT, PROP.cS, PROP.Tf, PROP.T_static);

% ---- Correcciones: góndola y Mach punta ----
Vcorr = 250;                   % ft/s (para las correcciones)
ncorr = 1800/60;               % *** 30 rev/s (corregido) ***
Dcorr = 8;  Sc = 5;  t_c = 0.12;  eta_p_nom = 0.85;

Jcorr_nom = Vcorr/(ncorr*Dcorr);
Jcorr     = Jcorr_nom*(1 - 0.329*Sc/Dcorr^2);

a0 = 1116;                       % ft/s
Vtip = pi*Dcorr*ncorr;           % ft/s
Mtip = Vtip/a0;

if Mtip > 0.89
    eta_p_corr = eta_p_nom - (Mtip - 0.89)*(0.16/(0.48 - 3*t_c));
else
    eta_p_corr = eta_p_nom;
end

fprintf('\n=== Hélice — correcciones ===\n');
fprintf('J_nom=%.3f  J_corr=%.3f  M_tip=%.3f  eta_nom=%.3f  eta_corr=%.3f\n', ...
        Jcorr_nom, Jcorr, Mtip, eta_p_nom, eta_p_corr);

% ---- Scrubbing, Cooling y Misc ----
PROP.rho0 = 0.002377;  PROP.rho = 0.0020;  PROP.sigma = PROP.rho/PROP.rho0;
PROP.CfSwet_washed = 30;                        % ft^2
eta_p_eff = eta_p_nom * (1 - (1.558/Dcorr^2)*(PROP.rho/PROP.rho0)*PROP.CfSwet_washed);

P_hp_cool = 550;  T_engineR = 500;  Vcool = 250;
Dq_cooling1 = (4.9e-7)*P_hp_cool*T_engineR^2/(PROP.sigma*Vcool);
Dq_cooling2 = (6e-8)  *P_hp_cool*T_engineR^2/(PROP.sigma*Vcool);

P_hp_misc = 550; 
Dq_misc1 = 2e-4  * P_hp_misc;
Dq_misc2 = 2.5e-5* P_hp_misc;

fprintf('\n=== Hélice — scrubbing, cooling y misc ===\n');
fprintf('eta_eff=%.3f | (D/q)_cool=%.5f..%.5f | (D/q)_misc=%.5f..%.5f\n', ...
    eta_p_eff, Dq_cooling1, Dq_cooling2, Dq_misc1, Dq_misc2);

%% =======================================================================
% --------------------------- FUNCIÓN ISA -------------------------------
%% =======================================================================

function [rho, a, mu] = isa_atm_english(h_ft)
% Propiedades ISA estándar (unidades inglesas)
% IN:  h_ft [ft]
% 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;
        P36 = P0*(T/T0)^(-g0/(L*R));
        P = P36*exp(-g0/(R*T)*(h_ft-36089));
    end

    rho = P/(R*T);
    a   = sqrt(gamma*R*T);

    % Sutherland a 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;  % -> slug/(ft·s)
end
