% =========================================================================
% Master_Optimization_Suite.m
% Combined MDO Optimization Tool & Presentation Dashboard
% =========================================================================
clear; clc; close all;

%% ========================================================================
% PART 1: THE MASTER BATCH OPTIMIZATION LOOP
% =========================================================================
aircraftList = {'737', 'C130', 'RQ4'}; 

fprintf('\n==================================================\n');
fprintf('       STARTING MASTER BATCH OPTIMIZATION         \n');
fprintf('==================================================\n');

opt_Span  = zeros(1, 3);
opt_Sweep = zeros(1, 3);
opt_AR    = zeros(1, 3);
opt_Taper = zeros(1, 3);
opt_Sref  = zeros(1, 3);
opt_LD    = zeros(1, 3);
target_Sref = zeros(1, 3);

labels_Span  = cell(1, 3);
labels_Sweep = cell(1, 3);
labels_AR    = cell(1, 3);

all_history_X = cell(1, 3);
all_history_F = cell(1, 3);

for i = 1:length(aircraftList)
    TargetAircraft = aircraftList{i};
    clear runWingSimulation; % Flushes the persistent physics memory 
    [~, ~, params] = objective_constraints();
    
    params.lb = [0, 0.28, 25, -5]; 
    params.ub = [45, 1.0, 150, 0]; 
    params.k_span = 0;    
    params.k_sweep = 0;   
    params.useVSP = true; 
    params.constrainFuelCapacity = true;
    params.constrainCL = false;
    params.badConstraint = []; 

    switch TargetAircraft
        case '737'
            params.TargetAircraftName = '737';
            params.V_cruise_kt = 450; 
            params.alt_ft = 35000;
            params.W_ZFW_lbf = 115000;
            params.fuel_total_lbf = 46000;
            params.target_Span = 112.6;
            params.target_Sweep = 25.0;
            params.target_AR = 9.4;
            params.c_root = 19.1;   
            params.kappa = 0.87;    
            params.tc_root = 0.15;
            params.tc_tip  = 0.10;
            params.N_z = 3.75;
        case 'C130'
            params.TargetAircraftName = 'C130';
            params.V_cruise_kt = 320; 
            params.alt_ft = 28000;
            params.W_ZFW_lbf = 115000; 
            params.fuel_total_lbf = 40000;
            params.target_Span = 132.6;
            params.target_Sweep = 0.0;
            params.target_AR = 10.1;
            params.c_root = 15.2;   
            params.kappa = 0.80;    
            params.tc_root = 0.18;
            params.tc_tip  = 0.12;
            params.N_z = 3.0;
        case 'RQ4'
            params.TargetAircraftName = 'RQ4';
            params.V_cruise_kt = 340; 
            params.alt_ft = 60000;
            params.W_ZFW_lbf = 25000;
            params.fuel_total_lbf = 15000;
            params.target_Span = 130.9;
            params.target_Sweep = 0.0;
            params.target_AR = 25.0;
            params.c_root = 7.1;    
            params.kappa = 0.98;    
            params.tc_root = 0.14;
            params.tc_tip  = 0.12;
            params.N_z = 2.5;
    end
    [objFun, conFun, params] = objective_constraints(params);
    
    popN = 150;
    eliteN = ceil(0.03*popN);
    options = optimoptions('ga', 'Display','off', 'PopulationSize', popN, ...
        'EliteCount', eliteN, 'CrossoverFraction', 0.75, ...
        'MutationFcn', {@mutationadaptfeasible, 0.15}, ...
        'MaxGenerations', 200, 'MaxStallGenerations', 500, ...
        'UseParallel', true, 'OutputFcn', @saveGAHistory, ...
        'NonlinearConstraintAlgorithm', 'penalty');
    fprintf('Optimizing %s...\n', TargetAircraft);
    [x_best, ~, ~, ~] = ga(objFun, 4, [], [], [], [], params.lb, params.ub, conFun, options);
    sim_opt = runWingSimulation(x_best, params);
    
    error_Sweep = abs(x_best(1) - params.target_Sweep);
    error_AR    = abs(sim_opt.AR - params.target_AR) / params.target_AR * 100;
    error_Span  = abs(x_best(3) - params.target_Span) / params.target_Span * 100;

    opt_Span(i)  = x_best(3);
    opt_Sweep(i) = x_best(1);
    opt_AR(i)    = sim_opt.AR;
    opt_Taper(i) = x_best(2);
    opt_Sref(i)  = sim_opt.S_ref_ft2;
    opt_LD(i)    = sim_opt.LD_cruise;
    target_Sref(i) = params.target_Span^2 / params.target_AR;
    
    if abs(x_best(3) - params.ub(3)) < 0.1
        labels_Span{i} = 'MAX LIMIT';
    else
        labels_Span{i} = sprintf('%.1f%% Diff', error_Span);
    end
    labels_Sweep{i} = sprintf('%.1f deg Diff', error_Sweep);
    labels_AR{i}    = sprintf('%.1f%% Diff', error_AR);

    all_history_X{i} = evalin('base', 'ga_history_designs');
    all_history_F{i} = evalin('base', 'ga_history_scores');
    
    fprintf('--------------------------------------------------\n');
    fprintf('RESULTS: %s\n', TargetAircraft);
    fprintf('Taper Optimized: %.2f\n', x_best(2));
    
    if abs(x_best(3) - params.ub(3)) < 0.1
        fprintf('Span  Baseline: %5.1f ft   | Optimized: %5.1f ft  | [MAX LIMIT REACHED]\n', params.target_Span, x_best(3));
    else
        fprintf('Span  Baseline: %5.1f ft   | Optimized: %5.1f ft  | Diff: %5.1f%%\n', params.target_Span, x_best(3), error_Span);
    end
    
    fprintf('Sweep Baseline: %5.1f deg  | Optimized: %5.1f deg | Diff: %5.1f deg\n', params.target_Sweep, x_best(1), error_Sweep);
    fprintf('AR    Baseline: %5.1f      | Optimized: %5.1f     | Diff: %5.1f%%\n', params.target_AR, sim_opt.AR, error_AR);
    fprintf('--------------------------------------------------\n\n');
end
fprintf('\nMULTI-AIRCRAFT OPTIMIZATION COMPLETE. Rendering Dashboard...\n');

%% ========================================================================
% PART 2: THE DASHBOARD RENDERER
% =========================================================================
aircraft = {'Boeing 737', 'C-130', 'RQ-4'};
X = categorical(aircraft);
X = reordercats(X, aircraft);

target_Span  = [112.6, 132.6, 130.9];
target_Sweep = [25.0, 0.0, 0.0];
target_AR    = [9.4, 10.1, 25.0];

% --- FIGURE 1: BAR CHARTS ---
fig1 = figure('Name', 'Optimization Bar Charts', 'Color', 'w');
fig1.Position(3) = 1300; fig1.Position(4) = 450;  
sgtitle('Aerodynamic Optimization vs Historical Baselines', 'FontSize', 16, 'FontWeight', 'bold');
col_T = [0.1, 0.2, 0.4]; col_O = [0.2, 0.6, 0.6];

subplot(1,3,1); b1 = bar(X, [target_Span', opt_Span'], 'grouped');
set(b1(1), 'FaceColor', col_T); set(b1(2), 'FaceColor', col_O);
title('Wingspan Optimization'); ylabel('Span (ft)'); grid on; 
legend('Baseline', 'Optimized', 'Location', 'northwest');
for i = 1:3
    h_max = max(target_Span(i), opt_Span(i));
    text(b1(2).XEndPoints(i), h_max+8, labels_Span{i}, 'HorizontalAlignment', 'center', 'Color', [0.6 0 0]);
end
ax1 = gca; ax1.YLimMode = 'manual'; ax1.YLim(1) = 0; ax1.YLim(2) = 200;

subplot(1,3,2); b2 = bar(X, [target_Sweep', opt_Sweep'], 'grouped');
set(b2(1), 'FaceColor', col_T); set(b2(2), 'FaceColor', col_O);
title('Sweep Angle Optimization'); ylabel('Sweep (deg)'); grid on; 
for i = 1:3
    h_max = max(target_Sweep(i), opt_Sweep(i));
    text(b2(2).XEndPoints(i), max(h_max+3, 6), labels_Sweep{i}, 'HorizontalAlignment', 'center', 'Color', [0.6 0 0]);
end
ax2 = gca; ax2.YLimMode = 'manual'; ax2.YLim(1) = 0; ax2.YLim(2) = 50;

subplot(1,3,3); b3 = bar(X, [target_AR', opt_AR'], 'grouped');
set(b3(1), 'FaceColor', col_T); set(b3(2), 'FaceColor', col_O);
title('Aspect Ratio Optimization'); ylabel('AR'); grid on; 
for i = 1:3
    h_max = max(target_AR(i), opt_AR(i));
    text(b3(2).XEndPoints(i), h_max+2, labels_AR{i}, 'HorizontalAlignment', 'center', 'Color', [0.6 0 0]);
end
ax3 = gca; ax3.YLimMode = 'manual'; ax3.YLim(1) = 0; ax3.YLim(2) = 45;

% --- FIGURE 2: DATA TABLE ---
fig2 = figure('Name', 'Optimization Data Table', 'Color', 'w');
fig2.Position(3) = 1150; fig2.Position(4) = 180;
arrow = char(8594); 

tableData = cell(3, 6);
for i = 1:3
    tableData{i,1} = sprintf('%.1f ft %s %.1f ft', target_Span(i), arrow, opt_Span(i));
    tableData{i,2} = sprintf('%.1f deg %s %.1f deg', target_Sweep(i), arrow, opt_Sweep(i));
    tableData{i,3} = sprintf('%.1f %s %.1f', target_AR(i), arrow, opt_AR(i));
    tableData{i,4} = sprintf('%.0f sq ft %s %.0f sq ft', target_Sref(i), arrow, opt_Sref(i));
    tableData{i,5} = sprintf('%.2f', opt_Taper(i));
    tableData{i,6} = sprintf('%.2f', opt_LD(i));
end
colNames = {'Wingspan (Baseline → Optimized)', 'Sweep Angle (Baseline → Optimized)', ...
            'Aspect Ratio (Baseline → Optimized)', 'Wing Area (Baseline → Optimized)', ...
            'Final Taper Ratio', 'Final Lift-to-Drag (L/D)'};
t = uitable(fig2, 'Data', tableData, 'ColumnName', colNames, 'RowName', aircraft);
t.Units = 'normalized'; t.Position = [0.02 0.1 0.96 0.8]; t.FontSize = 12;

% --- FIGURE 3: 3D SCATTER CONVERGENCE ---
fig3 = figure('Name', 'GA Convergence Profiles', 'Color', 'w');
fig3.Position(3) = 1400; fig3.Position(4) = 450;
sgtitle('Genetic Algorithm Convergence Profiles', 'FontSize', 16, 'FontWeight', 'bold');

for i = 1:3
    subplot(1, 3, i);
    hist_X = all_history_X{i};
    hist_F = all_history_F{i};
    
    valid_idx = hist_F < 0 & hist_F > -1e8; 
    X_plot = hist_X(valid_idx, :);
    Range_plot = -hist_F(valid_idx); 

    scatter3(X_plot(:,1), X_plot(:,3), Range_plot, 15, Range_plot, 'filled', 'MarkerEdgeColor', 'none'); hold on;
    [best_Range, best_idx] = max(Range_plot);
    plot3(X_plot(best_idx, 1), X_plot(best_idx, 3), best_Range, 'p', 'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'k', 'MarkerSize', 14);
    
    grid on; xlabel('Sweep (deg)', 'FontWeight', 'bold'); ylabel('Span (ft)', 'FontWeight', 'bold'); zlabel('Range (NM)', 'FontWeight', 'bold');
    title(sprintf('%s Optimization', aircraft{i}), 'FontSize', 12);
    colormap(turbo); view(45, 25);
end

% --- FIGURE 4: PLANFORM OVERLAYS ---
fig4 = figure('Color', 'w');
set(fig4, 'Name', 'Planform Optimization Overlays');
set(fig4, 'Position');
sgtitle('Historical Baseline vs. GA Optimized Wing Planforms', 'FontSize', 16, 'FontWeight', 'bold');

% Root Chords and historical tapers to anchor the drawings
c_root_list = [19.1, 15.2, 7.1];
taper_base_list = [0.254, 0.727, 0.475];

for i = 1:3
    ax = subplot(1, 3, i); hold on; grid on;
    
    % Package parameters for the drawing function
    data_base = [target_Span(i), target_Sweep(i), c_root_list(i), taper_base_list(i)];
    data_opt  = [opt_Span(i), opt_Sweep(i), c_root_list(i), opt_Taper(i)];
    
    h_base = drawWing(ax, data_base, 'baseline');
    h_opt  = drawWing(ax, data_opt, 'optimized');
    
    title(aircraft{i});
    formatAxis(ax);
    
    % Attach the detached global legend to the 3rd subplot
    if i == 3
        lgd = legend(ax, [h_base, h_opt], {'Historical Baseline', 'GA Optimized'}, 'Orientation', 'horizontal');
        lgd.Units = 'normalized';
        lgd.Position(1) = 0.5 - (lgd.Position(3) / 2); % Center horizontally
        lgd.Position(2) = 0.08; % Float it 8% off the bottom edge
    end
end

%% ========================================================================
% PART 3: LOCAL HELPER FUNCTIONS (The Engine Room)
% =========================================================================
function [objFun, conFun, params] = objective_constraints(customParams)
    if nargin < 1, params = defaultParams(); else, params = customParams; end
    objFun = @obj; conFun = @nonlcon;

    function f = obj(x)
        sim = runWingSimulation(x, params);
        if ~sim.ok, f = params.badFitness; return; end
        RangeNM = missionRangeJet(sim, params);
        lb_val = params.lb(:)'; ub_val = params.ub(:)';
        x_norm = (x - lb_val) ./ (ub_val - lb_val + eps);
        span_pen  = x_norm(3); sweep_pen = 1 - x_norm(1);
        J_pen = params.k_span * span_pen^2 + params.k_sweep * sweep_pen^2;
        f = -RangeNM + J_pen;
    end

    function [c, ceq] = nonlcon(x)
        c = []; ceq = [];
        sim = runWingSimulation(x, params);
        if ~sim.ok
            if isfield(params, 'constrainFuelCapacity') && params.constrainFuelCapacity
                c(end+1,1) = 1e6; 
            end
            c(end+1,1) = 1e6; return;
        end
        if isfield(params, 'constrainFuelCapacity') && params.constrainFuelCapacity
            c(end+1,1) = params.fuel_total_lbf - sim.W_fuel_capacity; 
        end
        max_stress = 15000; 
        if isfield(params, 'TargetAircraftName')
            switch params.TargetAircraftName
                case 'RQ4', max_stress = 25000; 
                case 'C130', max_stress = 15000; 
                case '737', max_stress = 15000; 
            end
        else
            max_stress = 15000; 
        end
        c(end+1,1) = sim.root_stress_psi - max_stress;
    end
end

function sim = runWingSimulation(x, params)
    SweepLE_deg = x(1); lambda = x(2); b_new_ft = x(3); twist_deg = x(4);
    persistent base const lastTarget
    if isempty(base) || isempty(lastTarget) || ~strcmp(lastTarget, params.TargetAircraftName)
        base = []; const = []; lastTarget = params.TargetAircraftName;
    end
    if isempty(base) || isempty(const)
        base.S_ref_ft2 = 1341; 
        if isfield(params, 'W_ZFW_lbf'), base.W_ZFW_lbf = params.W_ZFW_lbf; else, base.W_ZFW_lbf = 115000; end
        const.rho_sl   = 0.002377; 
        if isfield(params, 'alt_ft'), const.alt_ft = params.alt_ft; else, const.alt_ft = 35000; end
        [~, ~, ~, const.rho] = strato(const.alt_ft);
        const.V_kt     = params.V_cruise_kt; const.V_fps    = const.V_kt * 1.68781;
        const.q_psf    = 0.5 * const.rho * const.V_fps^2;
        if isfield(params, 'tc_root'), const.tc_root = params.tc_root; else, const.tc_root = 0.15; end
        if isfield(params, 'tc_tip'),  const.tc_tip  = params.tc_tip;  else, const.tc_tip  = 0.10; end
        const.cl_root  = 0.4; const.cl_tip   = 0.2;
        const.CD0_other = 0.0150; const.e = 0.85; 
        if isfield(params, 'N_z'), const.N_z = params.N_z; else, const.N_z = 3.75; end
    end
    SweepLE_rad = deg2rad(SweepLE_deg);
    S_ref_new_ft2 = b_new_ft * params.c_root * (1 + lambda) / 2;
    AR_new = b_new_ft^2 / S_ref_new_ft2;
    [tc_avg, cl_avg, vol_integ] = getIntegratedProps(b_new_ft, params.c_root, lambda, const.tc_root, const.tc_tip, const.cl_root, const.cl_tip);
    W_guess = base.W_ZFW_lbf + (params.fuel_total_lbf * 0.6); 
    CL_req = W_guess / (const.q_psf * S_ref_new_ft2);
    const.a = sqrt(1.4 * 1716 * (518.6 - 0.00356 * const.alt_ft)); 
    Mach = const.V_fps / const.a;
    if isfield(params, 'kappa'), kappa = params.kappa; else, kappa = 0.87; end
    M_dd = (kappa / cos(SweepLE_rad)) - (tc_avg / cos(SweepLE_rad)^2) - (CL_req / (10 * cos(SweepLE_rad)^3));
    CD_wave = 0;
    if Mach > M_dd, CD_wave = 50 * (Mach - M_dd)^2; end
    CD_i = CL_req^2 / (pi * AR_new * const.e);
    CD0_wing = 0.0080 * (1 + 2*tc_avg + 100*tc_avg^4); 
    CD_total = const.CD0_other + CD0_wing + CD_i + CD_wave;
    L_D_new = CL_req / CD_total;
    W_wing_new = 0.0051 * (W_guess * const.N_z)^0.557 * S_ref_new_ft2^0.649 * AR_new^0.5 * tc_avg^-0.4 * (1+lambda)^0.1 * (cos(SweepLE_rad))^-0.5 * b_new_ft^0.1;
    delta_W_wing = W_wing_new; 
    c_root_ft = (2 * S_ref_new_ft2) / (b_new_ft * (1 + lambda));
    c_root_in = c_root_ft * 12; t_root_in = c_root_in * tc_avg; b_new_in  = b_new_ft * 12;
    Lift_half_lbs = (W_guess * const.N_z) / 2; Moment_arm_in = (b_new_in / 2) * 0.45; 
    M_root_lbin = Lift_half_lbs * Moment_arm_in;
    I_root = 0.05 * c_root_in * (t_root_in^3); y_max = t_root_in / 2;
    sim.root_stress_psi = (M_root_lbin * y_max) / I_root; 
    sim.ok = true;
    if CL_req > 1.5 || CL_req < 0, sim.ok = false; end 
    sim.W_fuel_capacity = vol_integ * 50.0 * 0.60;
    sim.W_ZFW_lbf = base.W_ZFW_lbf + delta_W_wing;
    sim.fuel_total_lbf = params.fuel_total_lbf;
    sim.fuel_cruise_lbf = params.fuel_total_lbf * 0.9; 
    sim.W_start_lbf = sim.W_ZFW_lbf + params.fuel_total_lbf;
    sim.LD_cruise = L_D_new; sim.AR = AR_new; sim.span_ft = b_new_ft; sim.S_ref_ft2 = S_ref_new_ft2;
    sim.CL_req = CL_req; sim.CD_req = CD_total; sim.W_cruise_lbf = W_guess;

    function [T, P, m, rho] = strato(alt)
        if alt < 36089
            T = 518.6 - 0.00356 * alt; P = 2116 * (T / 518.6)^5.256;
        else
            T = 389.98; P = 472.68 * exp(1.73 - 0.000048 * alt);
        end
        rho = P / (1716 * T); m = 0; 
    end
end

function [tc_avg, cl_avg, vol_integ] = getIntegratedProps(b, c_root, lambda, tc_r, tc_t, cl_r, cl_t)
    n_steps = 20; dy = (b / 2) / n_steps; 
    y_mid = linspace(dy/2, (b/2) - dy/2, n_steps);
    y_frac = y_mid / (b/2); 
    c_y = c_root * (1 - (1 - lambda) * y_frac);
    tc_y = tc_r + (tc_t - tc_r) * y_frac;
    cl_y = cl_r + (cl_t - cl_r) * y_frac;
    dA = c_y * dy; S_ref_half = sum(dA);
    tc_avg = sum(c_y .* tc_y .* dy) / S_ref_half;
    cl_avg = sum(c_y .* cl_y .* dy) / S_ref_half;
    vol_integ = 2 * sum( (c_y.^2) .* tc_y .* dy );
end

function RangeNM = missionRangeJet(sim, params)
    if ~sim.ok, RangeNM = NaN; return; end
    W_start_lbf = sim.W_start_lbf; fuel_avail_lbf = sim.fuel_cruise_lbf;  
    V_fun = @(~) params.V_cruise_kt; TSFC_fun = @(~) params.TSFC_hr; LoverD_fun = @(~) sim.LD_cruise;
    
    R_nm = 0; W = W_start_lbf; fuel_left_lbf = fuel_avail_lbf;
    min_step = max(0.0005 * fuel_avail_lbf, 0.1); base_step = params.stepFrac * fuel_avail_lbf;
    while fuel_left_lbf > 0
        dW = min(max(base_step, min_step), fuel_left_lbf);
        if dW <= 0, break; end
        Wf_lbf = W - dW;
        if Wf_lbf <= 0, break; end
        R_nm = R_nm + (V_fun(W) / TSFC_fun(W)) * LoverD_fun(W) * log(W / Wf_lbf);
        W = Wf_lbf; fuel_left_lbf = fuel_left_lbf - dW;
    end
    RangeNM = R_nm;
    if isnan(RangeNM) || isinf(RangeNM) || ~isreal(RangeNM), RangeNM = NaN; end
end

function [state, options, optchanged] = saveGAHistory(options, state, flag)
    persistent all_designs all_scores
    optchanged = false;
    switch flag
        case 'init'
            total_individuals = options.PopulationSize * options.MaxGenerations;
            all_designs = zeros(total_individuals, size(state.Population,2));
            all_scores  = zeros(total_individuals, 1);
        case 'iter'
            start_idx = state.Generation * options.PopulationSize + 1;
            end_idx   = (state.Generation + 1) * options.PopulationSize;
            all_designs(start_idx:end_idx, :) = state.Population;
            all_scores(start_idx:end_idx)     = state.Score;
        case 'done'
            last_idx = (state.Generation + 1) * options.PopulationSize;
            assignin('base', 'ga_history_designs', all_designs(1:last_idx, :));
            assignin('base', 'ga_history_scores',  all_scores(1:last_idx, :));
            clear all_designs all_scores
    end
end

function params = defaultParams()
    params.V_cruise_kt = 350; params.TSFC_hr = 0.55; params.reserve_frac = 0.10; params.stepFrac = 0.01;
    params.useVSP = true; params.polarFolder = "."; params.polarFileName = "x58.polar"; 
    params.failIfOutsidePolar = true; params.badFitness = 1e9; params.badConstraint = 1e6;
    params.lb = [25, 0.2, 25, -5]; params.ub = [45, 0.6, 40, 0];
    params.k_span = 200; params.k_sweep = 100; params.constrainFuelCapacity = true; params.constrainCL = false;  
end

function h = drawWing(ax, data, type)
    span = data(1); sweep = data(2); c_root = data(3); taper = data(4);
    
    b2 = span / 2;
    x_tip_offset = b2 * tan(deg2rad(sweep));
    c_tip = c_root * taper;
    
    X = [0, b2, b2, 0, -b2, -b2];
    Y = [0, -x_tip_offset, -(x_tip_offset + c_tip), -c_root, -(x_tip_offset + c_tip), -x_tip_offset];
    
    if strcmp(type, 'baseline')
        h = patch(ax, X, Y, [0.85 0.85 0.85], 'EdgeColor', 'k', 'LineWidth', 1.5);
    else
        colorOpt = [0 0.447 0.741]; 
        h = patch(ax, X, Y, colorOpt, 'EdgeColor', colorOpt, 'LineWidth', 1.5, ...
                  'LineStyle', '--', 'FaceAlpha', 0.25);
    end
end

function formatAxis(ax)
    axis(ax, 'equal'); 
    xlabel(ax, 'Span (ft)', 'FontWeight', 'bold');
    ylabel(ax, 'Chord (ft)', 'FontWeight', 'bold');
    xlim(ax, [-85 85]);
    ylim(ax, [-45 5]);
end

function polar = loadVSPPolar(filename)
    % Preserved for future VSPAERO integration
    fid = fopen(filename, 'r'); if fid == -1, return; end
    headerLine = '';
    while true
        line = fgetl(fid); if ~ischar(line), fclose(fid); return; end
        if contains(line, 'AoA') && contains(line, 'CLtot') && contains(line, 'CDtot'), headerLine = strtrim(line); break; end
    end
    colNames = strsplit(regexprep(headerLine, '\s+', ' '), ' ');
    idxAoA = find(strcmp(colNames, 'AoA'), 1); idxCL = find(strcmp(colNames, 'CLtot'), 1); idxCD = find(strcmp(colNames, 'CDtot'), 1);
    alphaDeg = []; CL = []; CD = [];
    while true
        line = fgetl(fid); if ~ischar(line), break; end
        if isempty(strtrim(line)), continue; end
        nums = str2double(strsplit(strtrim(regexprep(line, '\s+', ' ')), ' '));
        if any(isnan(nums)) || numel(nums) < max([idxAoA, idxCL, idxCD]), continue; end
        alphaDeg(end+1,1) = nums(idxAoA); CL(end+1,1) = nums(idxCL); CD(end+1,1) = nums(idxCD); 
    end
    fclose(fid);
    [CL_sorted, idx] = sort(CL); polar.alphaDeg = alphaDeg(idx); polar.CL = CL_sorted; polar.CD = CD(idx);
end

function CD = getCDfromCL(CL_query, polar)
    % Preserved for future VSPAERO integration
    CLq = min(max(CL_query, min(polar.CL)), max(polar.CL));
    CD = interp1(polar.CL, polar.CD, CLq, 'pchip'); 
end