%% ================================================================
%  G500 — Aircraft Performance & Design (Study Script)
%  Autor: Antonio Perez
%  Unidades: ft, slug, s, deg (lbf implícito)
%  Qué aprenderás:
%    1) Cómo pasar de geometría y condiciones de crucero a CL(α)
%    2) Cómo construir CD0 por componentes (build-up)
%    3) Cómo comparar con Cfe y estimar la eficiencia de Oswald
%    4) Cómo interpretar Swet total, área equivalente Seq y la polar CD(CL)
% ================================================================

clear; clc; close all;

%% ================== (1) ENTRADAS DE DISEÑO =======================
% --- Condiciones de vuelo (crucero) ---
M          = 0.85;         % Mach: relación V/a (adimensional)
h_ft       = 45000;        % altitud [ft] (ISA estándar)

% --- Ala trapecial (wing) ---
b          = 86.333;       % envergadura [ft]
cR         = 19.0;         % cuerda raíz [ft]
cT         = 4.4;          % cuerda punta [ft]
Lambda_LE  = 30;           % flecha del borde de ataque [deg]
t_c        = 0.14;         % espesor relativo del ala (t/c)

% Datos extra usados en correlaciones de pendiente/forma (láminas):
x_cmaxtc   = 0.37;         % posición x/c del espesor máximo (~0.37)
cRexp      = 17.6;         % cuerda expuesta en raíz (para Sexp) [ft]
d_cut      = 7.3;          % recorte / fuselage cutout [ft]

% --- Estabilizador horizontal (H-tail) ---
bH         = 32.7;         % envergadura [ft]
cRH        = 9.6;          % cuerda raíz [ft]
cTH        = 4.2;          % cuerda punta [ft]
LambdaLE_H = 33;           % flecha LE [deg]
t_cH       = 0.10;         % espesor relativo

% --- Estabilizador vertical (V-tail) ---
bV         = 12.0;         % altura [ft]
cRV        = 14.1;         % cuerda raíz [ft]
cTV        = 9.0;          % cuerda punta [ft]
LambdaLE_V = 40;           % flecha LE [deg]
t_cV       = 0.10;         % espesor relativo

% --- Fuselaje ---
l_fus      = 78.8;         % longitud fuselaje [ft]
d_fus      = 7.3;          % diámetro fuselaje [ft]
S_top      = 483.4;        % área top proyectada usada en lámina [ft^2]
S_side     = 528.4;        % área side proyectada usada en lámina [ft^2]

% --- Góndolas (nacelles) ---
Dia_ave    = 4.6;          % diámetro promedio [ft]
l_nac      = 15.25;        % longitud góndola [ft]

% --- Parámetros "hoja" de la curva CL(α) ---
alpha_OL   = -4.25;        % ángulo de cero sustentación [deg]
Delta_alpha_CLmax = 2.7;   % corrimiento ángulo a CLmax [deg]
CLmax_airfoil = 2.23;      % cl_max del perfil a Re≈18e6
CLMAX_clmax   = 0.78;      % factor de conversión perfil→ala
CLMAX_CLMAX05 = 0.85;      % factor adicional
Lambda_tbl_deg  = [20 40]; % tabla ΔCLmax vs flecha
DeltaCL_tbl     = [-0.33 -0.15];

% --- Constantes de fricción/forma ---
ksmooth    = 2.08e-5;      % rugosidad equivalente (ft)
CFLAM_C    = 1.3828;       % constante para Cf_lam = C/sqrt(Re)
LP         = 1.10;         % 10% por "leakages & protuberances"
Q_tail     = 1.03;         % factor de interferencia de colas (H/V)
Q_wing     = 1.00;         % wing interference (≈1)
Q_nac      = 1.50;         % góndolas (alta interferencia)

%% ================== (2) GEOMETRÍA BÁSICA =========================
% Wing planform
lambda   = cT/cR;                               % razón de estrechamiento
Sref     = b/2 * (cR + cT);                     % área de referencia [ft^2]
AR       = b^2 / Sref;                          % alargamiento
cbar     = (2/3)*cR*(1 + lambda + lambda^2)/(1 + lambda);  % MAC [ft]
% Expuesta y mojada del ala (lámina)
Sexp_w   = (b - d_cut)/2 * (cRexp + cT);        % área expuesta [ft^2]
Swet_w   = Sexp_w * (1.977 + 0.52*t_c);         % área mojada ala [ft^2]

% H-tail
lambdaH  = cTH/cRH;
SH       = (bH - 1.0)/2 * (cRH + cTH);          % 1.0 ft de cutout (lámina)
ARH      = bH^2 / SH;
cbarH    = (2/3)*cRH*(1 + lambdaH + lambdaH^2)/(1 + lambdaH);
Swet_H   = SH * (1.977 + 0.52*t_cH);

% V-tail
lambdaV  = cTV/cRV;
SV       = cRV*bV*(1 + lambdaV)/2;
ARV      = bV^2 / SV;
cbarV    = (2/3)*cRV*(1 + lambdaV + lambdaV^2)/(1 + lambdaV);
Swet_V   = 2 * SV * (1.977 + 0.52*t_cV);        % *2 caras

% Fuselaje y góndolas
f_fus    = l_fus/d_fus;                         % fineza fuselaje
Swet_F   = 3.4 * ((S_top + S_side)/2);          % correlación de lámina
Swet_nac = 4*pi*Dia_ave*l_nac;                  % cilindro 4πDL

%% ================== (3) ATMÓSFERA, VELOCIDAD, REYNOLDS ===========
% Propiedades ISA a h_ft y velocidad de crucero V = M*a
[rho, a, mu] = isa_atm_english(h_ft);
V      = M * a;                                  % velocidad verdadera [ft/s]
Re_w   = rho*V*cbar / mu;                        % Re en MAC del ala
Re_H   = rho*V*cbarH/ mu;
Re_V   = rho*V*cbarV/ mu;

% "Cutoff" para rugosidad (Raymer)
Rec_w  = 38.21*(cbar / ksmooth)^1.053;
Rec_H  = 38.21*(cbarH/ ksmooth)^1.053;
Rec_V  = 38.21*(cbarV/ ksmooth)^1.053;

%% ================== (4) CURVA CL(α) A CRUCERO ====================
% 4.1  CL_max del ala (de la lámina con ΔCLmax tabulado)
DeltaCLmax  = interp1(Lambda_tbl_deg, DeltaCL_tbl, Lambda_LE, 'linear', 'extrap');
CLmax_wing  = CLMAX_CLMAX05 * (CLmax_airfoil*CLMAX_clmax + DeltaCLmax);

% 4.2  Pendiente CLα (Polhamus modificado, lámina)
beta   = sqrt(1 - M^2);                          % compresibilidad
% η: factor empírico (lámina con α_OL)
eta    = 57.3 * sqrt(1 - 0.15^2) / ((4 - alpha_OL)/(2*pi));
% Flecha a x/c_max
lambda_max = atan( tan(deg2rad(Lambda_LE)) * ( 4*x_cmaxtc*(1 - lambda)/(AR*(1 + lambda)) ) );
% FF (wing) dentro de CLα
FF_for_slope = 0.98;                             % "Use: FF = 0.98" (lámina)
CLalpha_rad = (2*pi*AR*FF_for_slope) / ...
              ( 2 + sqrt( 4 + (AR^2*beta^2/eta^2) * (1 + (tan(lambda_max)^2)/beta^2) ) );
CLalpha_deg = CLalpha_rad * (pi/180);            % [1/deg]

% 4.3  Ángulo a CLmax y curva
alpha_CLmax = CLmax_wing/CLalpha_deg + alpha_OL + Delta_alpha_CLmax;
alpha_vec   = linspace(alpha_OL-2, alpha_CLmax, 200);
CL_vec      = CLalpha_deg * (alpha_vec - alpha_OL);

%% ============= (5) CD0 POR COMPONENTES (BUILD-UP) ================
% Interpretación: CD0 suma pérdidas de fricción y forma a CL≈0.  
% Fórmulas tipo Raymer/Hoerner + ajustes de tus láminas.

% 5.1 Wing skin friction + form factor
Cf_lam_w  = CFLAM_C/sqrt(Re_w);
Cf_turb_w = 0.455/(log10(min(Re_w,Rec_w))^2.58 * (1 + 0.144*M^2)^0.65);
Cf_wing   = 0.65*Cf_turb_w + 0.35*Cf_lam_w;                         % mezcla 65/35
lambda_max_geom = lambda_max;                                       % reuso
FF_wing   = (1 + 0.6/x_cmaxtc*t_c + 100*t_c^4) * 1.34*M^0.18 * cos(lambda_max_geom)^0.28;
DeltaCD_wing = Cf_wing * FF_wing * Q_wing * (Swet_w/Sref);

% 5.2 Horizontal tail
Cf_lam_H  = CFLAM_C/sqrt(Re_H);
Cf_turb_H = 0.455/(log10(min(Re_H,Rec_H))^2.58 * (1 + 0.144*M^2)^0.65);
Cf_H      = 0.65*Cf_turb_H + 0.35*Cf_lam_H;
Lambda40H = atan( tand(LambdaLE_H) - 4*0.4*(1 - lambdaH)/(ARH*(1 + lambdaH)) );
FF_H      = (1 + 0.6/0.40*t_cH + 100*t_cH^4) * 1.34*M^0.18 * cos(Lambda40H)^0.28;
DeltaCD_H = Cf_H * FF_H * Q_tail * (Swet_H/Sref);

% 5.3 Vertical tail
Cf_lam_V  = CFLAM_C/sqrt(Re_V);
Cf_turb_V = 0.455/(log10(min(Re_V,Rec_V))^2.58 * (1 + 0.144*M^2)^0.65);
Cf_V      = 0.65*Cf_turb_V + 0.35*Cf_lam_V;
Lambda40V = atan( tand(LambdaLE_V) - 4*0.4*(1 - lambdaV)/(2*ARV*(1 + lambdaV)) );
FF_V      = (1 + 0.6/0.40*t_cV + 100*t_cV^4) * 1.34*M^0.18 * cos(Lambda40V)^0.28;
DeltaCD_V = Cf_V * FF_V * Q_tail * (Swet_V/Sref);

% 5.4 Fuselaje
Re_fus    = Re_V*(l_fus/cbarV);                                     % escala a longitud
Rec_fus   = 38.21*(l_fus/ksmooth)^1.053;
Cf_lam_F  = CFLAM_C/sqrt(Re_fus);
Cf_turb_F = 0.455/(log10(min(Re_fus,Rec_fus))^2.58 * (1 + 0.144*M^2)^0.65);
Cf_F      = Cf_turb_F;                                              % 100% turbulento (lámina)
FF_F      = 1 + 60/( (l_fus/d_fus)^3 ) + (l_fus/d_fus)/400;
DeltaCD_fus = Cf_F * FF_F * (Swet_F/Sref);

% 5.5 Góndolas (nacelles) — modelo 90% turb / 10% lam + Q_nac
Re_nac    = Re_w * (l_nac/cbar);
Rec_nac   = 38.21*(l_nac/ksmooth)^1.053;
Cf_lam_N  = CFLAM_C/sqrt(Re_nac);
Cf_turb_N = 0.455/(log10(min(Re_nac,Rec_nac))^2.58 * (1 + 0.144*M^2)^0.65);
Cf_N      = 0.9*Cf_turb_N + 0.1*Cf_lam_N;
FF_N      = 1 + 0.35/(l_nac/Dia_ave);
DeltaCD_nac = Cf_N * FF_N * Q_nac * (Swet_nac/Sref);

% 5.6 Upsweep del fuselaje (pérdida por base)
u_deg     = 3.4;
Amax      = 48.0;                                                   % [ft^2] de la lámina
DeltaCD_upswp = (3.83 * (u_deg^2.5) * Amax) / Sref;

% 5.7 CD0 total con "leakages & protuberances"
CD0_sum   = DeltaCD_wing + DeltaCD_H + DeltaCD_V + DeltaCD_fus + DeltaCD_nac + DeltaCD_upswp;
CD0_total = CD0_sum * LP;

%% ================ (6) ÁREAS: Swet total y Seq ====================
Swet_total = Swet_w + Swet_H + Swet_V + Swet_F + Swet_nac;
Seq        = CD0_total * Sref;                                      % área equivalente

%% ====== (7) COMPARACIÓN: MÉTODO Cfe Y EFICIENCIA DE OSWALD =======
Cfe           = 0.0026;                                             % "civil transport"
CD0_equivCfe  = Cfe * Swet_total / Sref;                            % CD0 alternativo

% Oswald e (dos fórmulas de lámina para estudiar su efecto)
e_straight    = 1.78*(1 - 0.045*AR^0.68) - 0.64;
e_swept       = 4.61*(1 - 0.045*AR^0.68) * (cosd(Lambda_LE))^0.15 - 3.1;
e_use         = e_swept;                                            % elige el que usarás

%% ============== (8) POLAR DE ARRASTRE PARA ESTUDIO ===============
% Interpretación: CD = CD0 + k CL^2, con k = 1/(pi e AR). Útil para L/D.
k_induced = 1/(pi*e_use*AR);
CL_grid   = linspace(0, 1.4, 200);
CD_grid   = CD0_total + k_induced*CL_grid.^2;

% Punto de máximo L/D (mínimo CD/CL): CL* = sqrt(CD0/k)
CL_star   = sqrt(CD0_total/k_induced);
CD_star   = CD0_total + k_induced*CL_star^2;
LD_max    = CL_star/CD_star;

%% ============================ SALIDAS =============================
fprintf('\n================= LIFT CURVE (cruise) =================\n');
fprintf('AR=%.2f, cbar=%.2f ft, Re_w=%.3e, beta=%.4f\n', AR, cbar, Re_w, beta);
fprintf('CLalpha = %.2f 1/rad (%.4f 1/deg)\n', CLalpha_rad, CLalpha_deg);
fprintf('CLmax_wing ≈ %.2f   alpha(CLmax) ≈ %.2f deg\n', CLmax_wing, alpha_CLmax);

fprintf('\n================= CD0 BUILD-UP (cruise) ===============\n');
fprintf('ΔCD_wing = %.5f | ΔCD_H = %.5f | ΔCD_V = %.5f\n', DeltaCD_wing, DeltaCD_H, DeltaCD_V);
fprintf('ΔCD_fus  = %.5f | ΔCD_nac = %.5f | ΔCD_upswp = %.5f\n', DeltaCD_fus, DeltaCD_nac, DeltaCD_upswp);
fprintf('CD0 sum (sin LP) = %.5f   ->  CD0 total (LP=%.2f) = %.5f\n', CD0_sum, LP, CD0_total);

fprintf('\n================= AREAS & Cfe COMPARISON ===============\n');
fprintf('Swet_total = %.1f ft^2   |   Seq = %.1f ft^2\n', Swet_total, Seq);
fprintf('CD0 (equiv, Cfe=%.4f) = %.4f   |  (build-up) = %.4f\n', Cfe, CD0_equivCfe, CD0_total);

fprintf('\n================= OSWALD & DRAG POLAR ==================\n');
fprintf('e_straight=%.3f, e_swept=%.3f -> e_use=%.3f\n', e_straight, e_swept, e_use);
fprintf('k_induced = %.4f | (L/D)_max ≈ %.1f at CL*=%.2f\n', k_induced, LD_max, CL_star);

%% ============================== PLOTS =============================
% 1) Lift curve
figure; plot(alpha_vec, CL_vec, 'LineWidth', 1.8); grid on;
xlabel('\alpha [deg]'); ylabel('C_L'); title('Lift Curve at Cruise');
xline(alpha_OL,'--'); yline(0,'-'); hold on;
plot(alpha_CLmax, CLmax_wing,'o','MarkerFaceColor','auto');
text(alpha_CLmax, CLmax_wing, sprintf('  C_{Lmax}=%.2f',CLmax_wing));

% 2) Drag polar (parabólica) y punto de (L/D)_max
figure; plot(CD_grid, CL_grid, 'LineWidth', 1.8); grid on;
xlabel('C_D'); ylabel('C_L'); title(sprintf('Drag Polar (e=%.3f, AR=%.2f)', e_use, AR));
hold on; plot(CD_star, CL_star, 'o', 'MarkerFaceColor','auto');
text(CD_star, CL_star, sprintf('  (L/D)_{max}=%.1f', LD_max));

%% ======================== FUNCIONES LOCALES =======================
function [rho, a, mu] = isa_atm_english(h_ft)
% isa_atm_english: atmósfera estándar "ligera" en unidades inglesas.
%  ENTRADA:  h_ft  altitud [ft]
%  SALIDA:   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;                             % capa con gradiente
        P = P0*(T/T0)^(-g0/(L*R));
    else
        T = 389.97;                                  % 36–65 kft isotérmica
        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);
    % Viscosidad dinámica (Sutherland) en 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



