





%% =================== G500 — Final Aero Constants ===================
% PURPOSE: Update the model with final Sref, CD0 and Oswald e,
%          then compute the drag polar and performance points.

% ---------------- User-provided final values ----------------
Sref   = 992.8;      % [ft^2] Wing reference area (sets the scale)
CD0    = 0.0142;     % [-] Zero-lift drag coefficient (parasite drag at CL≈0)
e      = 0.612;      % [-] Oswald span efficiency (accounts for non-elliptic lift & losses)

% ---------------- Geometry kept from earlier work -------------
AR     = 7.38;       % [-] Aspect ratio (update if your final AR differs)

% ---------------- Derived quantities --------------------------
k      = 1/(pi*e*AR);          % [-] Induced-drag factor in parabolic polar: CD = CD0 + k*CL^2
Seq    = CD0*Sref;             % [ft^2] Equivalent flat-plate area (useful magnitude check)

% -------- Best L/D point for a parabolic polar ----------------
CL_star = sqrt(CD0/k);         % [-] CL at max L/D (where d(CD/CL)/dCL = 0)
LD_max  = 1/(2*sqrt(CD0*k));   % [-] Maximum lift-to-drag ratio

% ----------------------- Reporting -----------------------------
fprintf('\n=== G500 Final Aero (Updated) ===\n');
fprintf('Sref   = %.1f ft^2\n', Sref);
fprintf('CD0    = %.5f\n', CD0);
fprintf('e      = %.3f\n', e);
fprintf('AR     = %.2f\n', AR);
fprintf('k      = %.5f  (induced-drag factor)\n', k);
fprintf('Seq    = %.2f ft^2 (equivalent area)\n', Seq);
fprintf('CL*    = %.3f at (L/D)max\n', CL_star);
fprintf('(L/D)max = %.2f\n', LD_max);

% ----------------------- Drag Polar plot -----------------------
CL = linspace(0, 1.4, 300);          % typical CL envelope
CD = CD0 + k*CL.^2;                   % parabolic polar

figure;
plot(CD, CL, 'LineWidth', 1.8); grid on; hold on;
xlabel('C_D'); ylabel('C_L');
title(sprintf('G500 Drag Polar (CD_0=%.4f, e=%.3f, AR=%.2f)', CD0, e, AR));
% mark best L/D point
CD_star = CD0 + k*CL_star^2;
plot(CD_star, CL_star, 'o', 'MarkerFaceColor','auto');
text(CD_star, CL_star, sprintf('  (L/D)_{max}=%.1f\\n  C_L^*=%.3f', LD_max, CL_star));

% (Optional) show a CL line typical of cruise if you want:
% CL_cruise = 0.45;  xline(CD0 + k*CL_cruise^2, '--');









%% ================================================================
% APD - Asignación 2: CD0 por Component Build-up (crucero)
% G500 aprox. (valores de tus láminas)
% Unidades: ft, slug, s, deg
%% ================================================================
clear; clc; close all;

%% -------------------- ENTRADAS DE CRUCERO ------------------------
M    = 0.85;          % Mach
h_ft = 45000;         % altitud [ft] ISA

% Wing (para Sref y algunas escalas)
b   = 86.333;         % envergadura [ft]
cR  = 19.0;           % cuerda raíz [ft]
cT  = 4.4;            % cuerda punta [ft]
t_c = 0.14;           % espesor relativo
x_cmaxtc = 0.37;      % x/c del espesor máx (para FF_wing)
LambdaLE = 30;        % flecha LE [deg]
cRexp = 17.6; d_cut = 7.3;

% Horizontal tail
bH=32.7; cRH=9.6; cTH=4.2; t_cH=0.10; LambdaLE_H=33; dH=1.0;

% Vertical tail
bV=12.0; cRV=14.1; cTV=9.0; t_cV=0.10; LambdaLE_V=40;

% Fuselaje
l_fus=78.8; d_fus=7.3; S_top=483.4; S_side=528.4;

% Góndolas
Dia_ave=4.6; l_nac=15.25;

% Constantes / factores
ksmooth   = 2.08e-5;  % rugosidad (ft)
CFLAM_C   = 1.3828;   % Cf_lam = C/sqrt(Re)
Q_tail    = 1.03;     % interferencia colas
Q_wing    = 1.00;     % wing
Q_nac     = 1.50;     % góndolas
LP        = 1.10;     % +10% leaks/protuberances

%% -------------------- GEOMETRÍA Y ÁREAS --------------------------
lambda  = cT/cR;
Sref    = b/2*(cR+cT);                 % área ala = referencia
AR      = b^2/Sref;
cbar    = (2/3)*cR*(1+lambda+lambda^2)/(1+lambda);

Sexp_w  = (b-d_cut)/2*(cRexp+cT);
Swet_w  = Sexp_w*(1.977+0.52*t_c);

lambdaH = cTH/cRH;
SH      = (bH-dH)/2*(cRH+cTH);
ARH     = bH^2/SH;
cbarH   = (2/3)*cRH*(1+lambdaH+lambdaH^2)/(1+lambdaH);
Swet_H  = SH*(1.977+0.52*t_cH);

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);

f_fus   = l_fus/d_fus;
Swet_F  = 3.4*((S_top+S_side)/2);

Swet_nac = 4*pi*Dia_ave*l_nac;

%% -------------------- ATMÓSFERA & REYNOLDS -----------------------
[rho,a,mu] = isa_atm_english(h_ft);
V = M*a;

Re_w = rho*V*cbar /mu;  Re_c_w = 38.21*(cbar /ksmooth)^1.053;
Re_H = rho*V*cbarH/mu;  Re_c_H = 38.21*(cbarH/ksmooth)^1.053;
Re_V = rho*V*cbarV/mu;  Re_c_V = 38.21*(cbarV/ksmooth)^1.053;

%% -------------------- COEF. FRICCIÓN Y FF POR COMPONENTE ---------
% Wing (mezcla 65/35)
Cf_lam_w  = CFLAM_C/sqrt(Re_w);
Cf_turb_w = 0.455/(log10(min(Re_w,Re_c_w))^2.58*(1+0.144*M^2)^0.65);
Cf_wing   = 0.65*Cf_turb_w + 0.35*Cf_lam_w;
lambda_max = atan( tan(deg2rad(LambdaLE))*(4*x_cmaxtc*(1-lambda)/(AR*(1+lambda))) );
FF_wing   = (1 + 0.6/x_cmaxtc*t_c + 100*t_c^4)*1.34*M^0.18*cos(lambda_max)^0.28;

% H-tail
Cf_lam_H  = CFLAM_C/sqrt(Re_H);
Cf_turb_H = 0.455/(log10(min(Re_H,Re_c_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;

% V-tail
Cf_lam_V  = CFLAM_C/sqrt(Re_V);
Cf_turb_V = 0.455/(log10(min(Re_V,Re_c_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;

% Fuselaje (100% turbulento en la hoja)
Re_fus    = Re_V*(l_fus/cbarV);      Re_c_fus = 38.21*(l_fus/ksmooth)^1.053;
Cf_turb_F = 0.455/(log10(min(Re_fus,Re_c_fus))^2.58*(1+0.144*M^2)^0.65);
FF_F      = 1 + 60/f_fus^3 + f_fus/400;

% Nacelles (90% turb + 10% lam)
Re_nac    = Re_w*(l_nac/cbar);       Re_c_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,Re_c_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);

%% -------------------- ΔCD0 DE CADA COMPONENTE --------------------
DeltaCD_wing = Cf_wing*FF_wing*Q_wing*(Swet_w/Sref);
DeltaCD_H    = Cf_H   *FF_H   *Q_tail*(Swet_H/Sref);
DeltaCD_V    = Cf_V   *FF_V   *Q_tail*(Swet_V/Sref);
DeltaCD_fus  = Cf_turb_F*FF_F*(Swet_F/Sref);
DeltaCD_nac  = Cf_N   *FF_N   *Q_nac*(Swet_nac/Sref);

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

% Suma y corrección LP
CD0_sum   = DeltaCD_wing + DeltaCD_H + DeltaCD_V + DeltaCD_fus + DeltaCD_nac + DeltaCD_upswp;
CD0_total = CD0_sum * LP;

%% -------------------- COMPARACIÓN CON Cfe ------------------------
Swet_total = Swet_w + Swet_H + Swet_V + Swet_F + Swet_nac;
Cfe = 0.0026;                   % "Civil Transport"
CD0_equiv = Cfe * Swet_total / Sref;

%% -------------------- PRINT ORDENADO ------------------------------
fprintf('\n========== Asignación 2 — CD0 @ Cruise ==========\n');
fprintf('Sref=%.1f ft^2 | AR=%.2f | cbar=%.2f ft | V=%.1f ft/s\n', Sref, AR, cbar, V);
fprintf('Re_w=%.3e  Re_H=%.3e  Re_V=%.3e\n', Re_w, Re_H, Re_V);

fprintf('\n-- Componentes ΔCD --\n');
fprintf('Wing      : %.5f\n', DeltaCD_wing);
fprintf('H-tail    : %.5f\n', DeltaCD_H);
fprintf('V-tail    : %.5f\n', DeltaCD_V);
fprintf('Fuselage  : %.5f\n', DeltaCD_fus);
fprintf('Nacelles  : %.5f\n', DeltaCD_nac);
fprintf('Upsweep   : %.5f\n', DeltaCD_upswp);
fprintf('--------------------------------\n');
fprintf('CD0 (sin LP)  = %.5f\n', CD0_sum);
fprintf('CD0 (LP=%.2f) = %.5f   <-- RESULTADO ASIG. 2\n', LP, CD0_total);

fprintf('\n-- Swet total y Cfe --\n');
fprintf('Swet_total = %.1f ft^2 | CD0_equiv(Cfe=%.4f) = %.4f\n', Swet_total, Cfe, CD0_equiv);

%% -------------------- GRÁFICA OPCIONAL ---------------------------
labels = {'Wing','H-tail','V-tail','Fuselage','Nacelles','Upsweep'};
vals   = [DeltaCD_wing,DeltaCD_H,DeltaCD_V,DeltaCD_fus,DeltaCD_nac,DeltaCD_upswp];
figure; bar(categorical(labels), vals,'LineWidth',1.0);
ylabel('\DeltaC_D'); title('CD_0 por componente @ cruise'); grid on;
text(1:length(vals), vals, compose('%.4f',vals),'HorizontalAlignment','center','VerticalAlignment','bottom');



%% ================= FUNCION ATMOSFERA ISA ======================
function [rho, a, mu] = isa_atm_english(h_ft)
% isa_atm_english: calcula propiedades ISA estándar en unidades inglesas
% Entrada:
%   h_ft - altitud [ft]
% Salida:
%   rho - densidad [slug/ft^3]
%   a   - velocidad del sonido [ft/s]
%   mu  - viscosidad dinámica [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);

    % 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



%% ================================================================
%% ================================================================
%% ================================================================
%% ================================================================



%% ================================================================
% Propulsión por chorro simple — fórmulas de empuje y potencia
% (Fuente: Aircraft Performance)
% Sirve para ilustrar la teoría de momento en la asignación
%% ================================================================

% --- Parámetros de ejemplo (puedes editarlos según caso de estudio) ---
rho   = 0.002377;   % densidad [slug/ft^3] ~ nivel del mar ISA
S_jet = 3.0;        % área de sección de chorro [ft^2]
V_jet = 1500;       % velocidad de chorro [ft/s]
V0    = 800;        % velocidad de vuelo [ft/s]

% --- Masa de aire que atraviesa el chorro por segundo ---
mdot = rho * V_jet * S_jet;   % [slug/s]

% --- 1) Empuje neto (F) ---
F = mdot*(V_jet - V0);        % [lbf] si trabajas en sistema inglés
% >> Este es el empuje producido por acelerar el flujo del aire.

% --- 2) Potencia de tracción (Pt) ---
Pt = F * V0;                  % [ft·lbf/s]
% >> Es la potencia útil asociada al empuje a la velocidad de vuelo.

% --- 3) Potencia expandida (Pexp) ---
Pexp = 0.5*mdot*(V_jet^2 - V0^2);   % [ft·lbf/s]
% >> Es la tasa de incremento de energía cinética del flujo.

% --- Mostrar resultados ---
fprintf('\n======= Propulsión (teoría de momento) =======\n');
fprintf('mdot   = %.2f slug/s\n', mdot);
fprintf('F      = %.2f lbf (Empuje)\n', F);
fprintf('Pt     = %.2f ft·lbf/s (Potencia útil)\n', Pt);
fprintf('Pexp   = %.2f ft·lbf/s (Potencia expandida)\n', Pexp);



%% -------------------- Eficiencia propulsiva ---------------------
% Definición:
% eta_PE = Pt / Pexp
% También: eta_PE = 2 / (V/V0 + 1)

eta_PE = Pt / Pexp;
eta_PE_alt = 2 / (V_jet/V0 + 1);

% --- Mostrar resultados ---
fprintf('\n======= Eficiencia propulsiva =======\n');
fprintf('eta_PE (definición) = %.4f\n', eta_PE);
fprintf('eta_PE (fórmula cerrada) = %.4f\n', eta_PE_alt);


%% ================================================================
% Propulsión — Parte Power Operation (Specific Fuel Consumption)
% Ecuación de aproximación para c/c_max_dry según T/Tmax_dry y Mach
% Fuente: ME4932 Aircraft Performance & Design
%% ================================================================

% --- Parámetros de entrada ---
T_frac = linspace(0.2,1,50);    % razón T/Tmax_dry (20% a 100%)
M = 0.85;                       % número de Mach
c_max_dry = 0.5;                % valor ejemplo de SFC max dry [lbm/(lbf·hr)]

% --- Ecuación de aproximación ---
c_ratio = (0.1 ./ (T_frac)) ...
        + (0.24 ./ (T_frac.^0.8)) ...
        + (0.66 .* (T_frac.^0.8)) ...
        + (0.1*M).*(1./T_frac - T_frac);

% --- SFC en condiciones parciales ---
c_partial = c_ratio * c_max_dry;

% --- Mostrar resultados numéricos para algunos puntos ---
fprintf('\n======= Propulsión — Parte Power Operation =======\n');
for i = [1 10 20 30 40 50]
    fprintf('T/Tmax=%.2f  |  c/cmax=%.3f  |  c=%.3f [lbm/(lbf·hr)]\n', ...
        T_frac(i), c_ratio(i), c_partial(i));
end

% --- Gráfica ---
figure;
plot(T_frac, c_ratio,'LineWidth',2); grid on;
xlabel('T/T_{max,dry}'); ylabel('c/c_{max,dry}');
title(sprintf('SFC ratio vs Thrust Fraction (M=%.2f)',M));



%% ================================================================
% Propulsión — Motores de Pistón
% Potencia disponible en función de densidad relativa
% Fuente: ME4932 Aircraft Performance & Design
%% ================================================================

% --- Parámetros ---
Power_SL = 200;    % [hp] ejemplo de potencia al nivel del mar
h_vec = linspace(0,25000,50);   % altitud de 0 a 25,000 ft
rho0 = 0.002377;   % densidad nivel del mar [slug/ft^3]

% --- Calcular densidad atmosférica (usando función isa_atm_english) ---
rho_vec = zeros(size(h_vec));
for i=1:length(h_vec)
    [rho,~,mu] = isa_atm_english(h_vec(i));
    rho_vec(i) = rho;
end

% --- Relación de densidad ---
sigma = rho_vec / rho0;   % adimensional

% --- Potencia de pistón a cada altitud ---
Power = Power_SL .* ( sigma - (1 - sigma)/7.55 );

% --- Resultados ---
fprintf('\n======= Propulsión — Piston Engines =======\n');
for i = 1:10:length(h_vec)
    fprintf('h=%.0f ft | sigma=%.3f | Power=%.1f hp\n', ...
        h_vec(i), sigma(i), Power(i));
end

% --- Gráfica ---
figure;
plot(h_vec, Power,'LineWidth',2); grid on;
xlabel('Altitud [ft]'); ylabel('Potencia disponible [hp]');
title('Motor de Pistón: Potencia disponible vs Altitud');




%% ================================================================
% Propulsión — Análisis de hélice
% Eficiencia = Thrust Power / Shaft Power
% Fuente: ME4932 Aircraft Performance & Design
%% ================================================================

% --- Parámetros de ejemplo ---
T = 2000;        % [lbf] empuje de la hélice
V = 300;         % [ft/s] velocidad de vuelo (~180 kt)
Ps_hp = 1200;    % potencia al eje en HP
Ps = Ps_hp * 550; % convertir HP a ft·lbf/s

% --- Cálculo de eficiencia ---
ThrustPower = T * V;           % potencia útil (ft·lbf/s)
eta_p = ThrustPower / Ps;      % eficiencia propulsiva de la hélice

% --- Mostrar resultados ---
fprintf('\n======= Propulsión — Hélice =======\n');
fprintf('Empuje T = %.1f lbf\n', T);
fprintf('Velocidad V = %.1f ft/s\n', V);
fprintf('Shaft Power = %.1f hp\n', Ps_hp);
fprintf('Thrust Power = %.1f hp\n', ThrustPower/550);
fprintf('Eficiencia de hélice eta_p = %.3f (%.1f%%)\n', eta_p, eta_p*100);

% --- Exploración de eficiencia vs velocidad ---
V_vec = linspace(100,500,50);       % rango de velocidades [ft/s]
eta_curve = (T .* V_vec) ./ Ps;     % eficiencia para cada V

figure;
plot(V_vec, eta_curve,'LineWidth',2); grid on;
xlabel('Velocidad [ft/s]'); ylabel('\eta_p');
title('Eficiencia de hélice vs Velocidad');
ylim([0 1]);






%% ================================================================
% Propulsión — Hélices: Coeficientes adimensionales
% Fuente: ME4932 Aircraft Performance & Design
%% ================================================================

% --- Parámetros de ejemplo ---
rho = 0.002377;    % slug/ft^3 (densidad a nivel del mar)
V = 250;           % ft/s (velocidad de vuelo)
n = 30/60;         % rev/s (ejemplo: 1800 RPM)
D = 8;             % ft (diámetro hélice)
P_hp = 550;        % potencia entregada [hp]
P = P_hp * 550;    % convertir a ft·lbf/s
eta_p = 0.85;      % eficiencia propulsiva supuesta

% --- Coeficientes ---
J  = V/(n*D);                          % advance ratio
cP = P/(rho*n^3*D^5);                  % power coefficient
T   = P*eta_p/V;                       % thrust en vuelo
cT = T/(rho*n^2*D^4);                  % thrust coefficient
cS = V^5*sqrt(rho/(P*n^2));            % speed-power coefficient

% --- Thrust estático (usando cT/cP) ---
T_static = (cT/cP) * (P/(n*D));

% --- Mostrar resultados ---
fprintf('\n======= Propeller Coefficients =======\n');
fprintf('Advance ratio J = %.3f\n', J);
fprintf('Power coefficient cP = %.3f\n', cP);
fprintf('Thrust coefficient cT = %.3f\n', cT);
fprintf('Speed-power coefficient cS = %.3f\n', cS);
fprintf('Thrust (forward) = %.1f lbf\n', T);
fprintf('Thrust (static)  = %.1f lbf\n', T_static);









%% ================================================================
% Propulsión — Hélices: Correcciones por góndolas y Mach en la punta
%% ================================================================

% --- Parámetros de ejemplo ---
V = 250;        % ft/s velocidad del avión
n = 30/60;      % rev/s (ejemplo: 1800 RPM)
D = 8;          % diámetro hélice [ft]
Sc = 5;         % área de góndola bloqueando flujo [ft^2]
t_c = 0.12;     % espesor relativo de pala
eta_p_nom = 0.85; % eficiencia nominal sin corrección

% --- Advance ratio nominal ---
J = V/(n*D);

% --- Corrección por góndola ---
J_corr = J*(1 - 0.329*Sc/D^2);

% --- Mach en la punta ---
a = 1116;          % ft/s, vel. sonido a nivel del mar
Vtip = pi*D*n;     % vel. punta (circunferencia*rev/s)
Mtip = Vtip/a;

% --- Corrección por Mach en la punta ---
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

% --- Mostrar resultados ---
fprintf('\n======= Propeller Corrections =======\n');
fprintf('J nominal   = %.3f\n', J);
fprintf('J corrected = %.3f\n', J_corr);
fprintf('M_tip       = %.3f\n', Mtip);
fprintf('eta_p (nominal)  = %.3f\n', eta_p_nom);
fprintf('eta_p (corrected)= %.3f\n', eta_p_corr);




%% ================================================================
% Propulsión — Hélices: Scrubbing drag y Cooling drag
%% ================================================================

% --- Parámetros de ejemplo ---
eta_p = 0.85;        % eficiencia nominal
D = 8;               % diámetro de hélice [ft]
rho0 = 0.002377;     % slug/ft^3 densidad nivel del mar
rho = 0.0020;        % slug/ft^3 densidad a altitud
CfSwet_washed = 30;  % [ft^2] contribución de superficies bañadas
P_hp = 550;          % potencia [hp]
P = P_hp;            % se usa directamente en bhp
T_engine = 500;      % temperatura característica [R]
sigma = rho/rho0;
V = 250;             % velocidad [ft/s]

% --- 1) Eficiencia efectiva con scrubbing drag ---
eta_p_eff = eta_p * (1 - (1.558/D^2)*(rho/rho0)*CfSwet_washed);

% --- 2) Cooling drag (dos formas equivalentes) ---
Dq_cooling1 = (4.9e-7)*P*T_engine^2/(sigma*V);
Dq_cooling2 = (6e-8)*P*T_engine^2/(sigma*V);

% --- Resultados ---
fprintf('\n======= Scrubbing & Cooling Drag =======\n');
fprintf('Eficiencia nominal eta_p = %.3f\n', eta_p);
fprintf('Eficiencia efectiva eta_p_eff = %.3f\n', eta_p_eff);
fprintf('Cooling drag (form 1) D/q = %.5f\n', Dq_cooling1);
fprintf('Cooling drag (form 2) D/q = %.5f\n', Dq_cooling2);



%% ================================================================
% Propulsión — Hélices: Miscellaneous Engine Drag
%% ================================================================

% --- Parámetros de ejemplo ---
P_hp = 550;   % potencia en hp (ejemplo de motor)
P = P_hp;     % asumimos P en bhp directamente

% --- Cálculo del drag misc ---
Dq_misc1 = 2e-4 * P;   % forma en función de bhp
Dq_misc2 = 2.5e-5 * P; % forma alternativa equivalente

% --- Resultados ---
fprintf('\n======= Miscellaneous Engine Drag =======\n');
fprintf('(D/q)_misc (form 1) = %.5f\n', Dq_misc1);
fprintf('(D/q)_misc (form 2) = %.5f\n', Dq_misc2);
