SS_Airfoil_Obsiana_1001900496.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                       Supersonic Airfoil Design                         %
%                               FALL 2025                                 %
%                    MAE-3303-001 COMPRESSIBLE FLOW                       %
%                       Code by Aaron A. Obsiana                          %
%                          UTA ID: 1001900496                             %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear global;
N = 50;    % Number of iterations to generate.
nsu = 35;   % Number of SEGMENTS on the UPPER surface
nsl = 40;   % Number of SEGMENTS on the LOWER surface

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                           DESIGN CONSTRAINTS                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Minf =                            2.0;    % Freestream Mach Number
cmle_less_than_equal_to =         0.1;    % Constraint for abs(Cm) <= #
cl_greater_equal_to =             0.25;   % Constraint for Cl >= #
max_thickness =                   0.1;    % Maximum Thickness for #%*chord
minimum_AoA =                     1.0;    % Minimum AoA Constraint
maximum_AoA =                     10.0;   % Maximum AoA Constraint
chord =                           1.0;    % Airfoil Chord Length

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fprintf(['\n\n\n- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ' ...
    '- - - - -\n']);
fprintf(['Supersonic Airfoil Design Program for MAE 3303-001 Compressible ' ...
    'Flow\n']);
fprintf(['- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ' ...
    '- - - -\n']);
fprintf(['Code by Aaron Obsiana, The University of Texas at Arlington, FALL ' ...
    '2025\n']);
fprintf('P-code, clcdcm.p, provided by Dr. Zhen Xue Han\n');
fprintf(['- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -' ...
    ' - - -\n']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Below generates tu,tl,xtu,xtl values that will be used to find the 
% optimal values through a 'linearly generated' set, based on the required
% design constraints using the 'linspace' function:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% UPPER THICKNESS AND LENGTH-TO-THICKNESS VALUES GENERATED USING 'LINSPACE'
% WHILE ALSO SATISFYING THE DESIGN VARIABLE PARAMETERS:
tu_gen = linspace(0.05, 0.1, N); 
tl_gen = linspace(0.05, 0.1, N);
xtu_gen = linspace((1/3)*chord, (2/3)*chord, N);
xtl_gen = linspace((1/3)*chord, (2/3)*chord, N);
xu = linspace(0, chord, N);
xl = xu;

% USING 'MESHGRID' FUNCTION TO FIND ALL VALID COMBINATIONS TO SATISFY THE
% MAXIMUM THICKNESS CONSTRAINT FOR THE NESTED FOR-LOOP:
[tu_combinations, xtu_combinations] = meshgrid(tu_gen, xtu_gen);
[tl_combinations, xtl_combinations] = meshgrid(tl_gen, xtl_gen);


% Prealocating all the valid results to be stored (for optimization)
results = zeros([]);
valid_results = 1;
valid = false;
design_number = 0;
valid_designs = zeros([]);
best_designs = 50;
CL_design = zeros([]);
CD_design = zeros([]);

fprintf('50 TOP VALID DESIGNS LOADING AND GENERATING...\n')
% NESTED FOR-LOOP
for i = 1:N
    for j = 1:N
                tu = tu_combinations(i,j); % LOOPS TU COMBINATIONS
                tl = tl_combinations(i,j); % LOOPS TL COMBINATIONS
                xtu = xtu_combinations(j);
                xtl = xtl_combinations(i);

        if (tu + tl) < max_thickness*chord % MAX THICKNESS CONSTRAINT 
            continue
        end
            
            % Y-UPPER AND Y-LOWER VALUES CALCULATED AS AN ARRAY USING THE
            % PCHIP POLYNOMIAL EQUATION IN THE PROBLEM STATEMENT:
            yu_gen = ((tu./xtu).*(xu./xtu)).*((1-2.*xtu).*xu.^2 ...
                + ((3.*xtu.^2)-1).*xu + ...
                (xtu.*(2-3.*xtu))./((1-xtu).^2));

            yl_gen = (((-tl)./xtl).*(xl./xtl)).*((1-2.*xtl).*xl.^2 ...
                + ((3.*xtl.^2)-1).*xl ...
                + (xtl.*(2-3.*xtl))./((1-xtl).^2));

            [yu, yl] = meshgrid(yu_gen, yl_gen);

            % FOR LOOP TO FIND AOA THATS VALID FOR THE DESIGN.
            for alphad = linspace(minimum_AoA,maximum_AoA,N)

                % clcdcm.p p-code below utilized.
                [cl,cd,cm,xmu,cpu,xml,cpl]=clcdcms(Minf,alphad,nsu,xu,yu,nsl,xl,yl);

                CL_design(i) = cl; % USED FOR PLOTTING CL VS AOA
                CD_design(i) = cd; % USED FOR PLOTTING CL VS AOA

                % IF STATEMENTS TO FILTER OUT CL AND CM VALUES THAT DO NOT
                % SATISFY THE DESIGN CONSTRAINTS
            if isempty(cl) || isempty(cd) || isempty(cm)
                continue
            end

            if abs(cm) > cmle_less_than_equal_to
                continue
            end

            if cl >= cl_greater_equal_to
                valid = true;
                
                LD = cl/cd;
                design_number = design_number + 1;

                % STORES ALL VALID RESULTS FOR OPTIMIZATION AND PRINTOUT
                results(valid_results).design_number = design_number;
                results(valid_results).cl = cl;
                results(valid_results).cd = cd;
                results(valid_results).cm = cm;
                results(valid_results).AoA = alphad;
                results(valid_results).tu = tu;
                results(valid_results).tl = tl;
                results(valid_results).xtu = xtu;
                results(valid_results).xtl = xtl;
                results(valid_results).LD = LD;
                valid_results = valid_results + 1;
                

                design_cl = [results.cl];
                design_cd = [results.cd];
                design_cm = [results.cm];
                design_AoA = [results.AoA];
                design_tu = [results.tu];
                design_tl = [results.tl];
                design_xtu = [results.xtu];
                design_xtl = [results.xtl];
                design_LD = [results.LD];

                valid_designs(design_number, :) = [design_number, cl, ...
                    cd, cm, LD, alphad, tu, tu, xtu, xtl];
            end
            end
    end
end




% Filters out the numerous valid designs to find the 50 most optimal airfoils
% and also finds the best design out of ALL valid designs.
if exist('valid_designs', 'var') && ~isempty(valid_designs)
    [~, design_order] = sort(valid_designs(:,5), 'descend');
    sorted_designs = valid_designs(design_order, :);

    [~, least_ordered] = sort(valid_designs(:,5), 'ascend');
    worst_designs = valid_designs(least_ordered, :);

    top50 = min(50, size(sorted_designs, 1));
    lowest_bunch = min(1, size(sorted_designs, 1));
    best_designs = sorted_designs(1:top50,:);
    worst_design = worst_designs(1:lowest_bunch,:);


    for n = 1:top50
        fprintf(['|| FEASIBLE DESIGN #%5d || Cl = %.5f | Cd = %.5f | Cm = %.5f |'...
            ' LD = %.5f | AoA = %.2f | tu = %.5f | tl = %.5f | xtu = %.5f | xtl = %.5f |\n'], ...
            best_designs(n,1), best_designs(n,2), best_designs(n,3), best_designs(n,4), ...
            best_designs(n,5), best_designs(n,6), best_designs(n,7), best_designs(n,8), ...
            best_designs(n,9), best_designs(n,10))
    end

    chosen_design = best_designs(1,:);
    minimum_design = worst_design(1,:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                    PLOTTING PORTION OF THE CODE                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PRINTOUTS THE 'chosen_design' FUNCTION THAT IS THE MOST OPTIMAL OUT OF
% ALL THE COMBINATIONS GENERATED. ALSO CONDUCTS PERFORMANCE ANALYSIS FOR
% FINDING OSW AND PM EXPANSION.
    
    writematrix(best_designs, 'best_designs.xlsx');
    writematrix(chosen_design, 'optimal_design.xlsx');
    writematrix(valid_designs, 'valid_designs.xlsx');
    
    % PRINTS OUT OPTIMIZED DESIGN
    fprintf('\n|| OPTIMIZED DESIGN - VALID DESIGN #%d ||\n', chosen_design(1));
    fprintf('Cl = %.5f, Cd = %.5f, Cm = %.5f\n', chosen_design(2), chosen_design(3), chosen_design(4));
    fprintf('LD = %.5f\n', chosen_design(5));
    fprintf('Angle-of-Attack = %.1f\n', chosen_design(6));
    fprintf('tu = %.5f, tl = %.5f\n', chosen_design(7), chosen_design(8));
    fprintf('xtu = %.5f, xtl = %.5f\n', chosen_design(9), chosen_design(10));
    
    % OPTIMIZED DESIGN VALUES STORED TO BE USED TO PLOT
    tu_optimal = chosen_design(7);
    tl_optimal = chosen_design(8);
    xtu_optimal = chosen_design(9);
    xtl_optimal = chosen_design(10);
    xu_panels = linspace(0, chord, nsu + 1);
    xl_panels = linspace(0, chord, nsl + 1);

    yu_optimal = ((tu_optimal ./ xtu_optimal) .* (xu_panels ./ xtu_optimal)) .* ... 
             ((((1 - (2 .* xtu_optimal)).* (xu_panels.^2)) + ...
             (((3 .* (xtu_optimal.^2)) - 1).* xu_panels) + (xtu_optimal .* ...
             (2 - (3.*xtu_optimal)))) ./ ((1 - xtu_optimal).^2));

    yl_optimal = (((-tl_optimal ./ xtl_optimal) .* (xl_panels ./ xtl_optimal)) .* ... 
             ((((1 - (2 .* xtl_optimal)).* (xl_panels.^2)) + ...
             (((3 .* (xtl_optimal.^2)) - 1).* xl_panels) + (xtl_optimal .* ...
             (2 - (3.*xtl_optimal)))) ./ ((1 - xtl_optimal).^2)));

    % TO CALCULATE INCLINATION ANGLES AND STORE THEM TO USE IN THE
    % REPORT. YNODES INTERPOLATED BY 1D. MATLAB 'PCHIP' FUNCTION IS USED 
    % TO MAKE CALCULATIONS EASIER WITHOUT HAVING TO USE THE PCHIP EQUATION
    % PROVIDED IN THE PROJECT.
    upper_ynodes = interp1(xu, yu_gen, xu_panels, 'pchip');
    lower_ynodes = interp1(xl, yl_gen, xl_panels, 'pchip');
            
    % SOLVES FOR THE PANEL SLOPES AND ITS RESPECTIVE ANGLES.
    dxu = diff(xu_panels);
    dxl = diff(xl_panels);
    dyu = diff(upper_ynodes);
    dyl = diff(lower_ynodes);

    % INVERSE TANGENT TO GET DEFLECTION ANGLES.
    theta_upper = atan2(dyu, dxu);
    theta_lower = atan2(dyl, dxl);

    % CONVERTING AOA FROM DEGREES TO RADIANS FOR USE IN OSW AND PM
    % RELATION EQUATIONS.
    alphar = deg2rad(alphad);

    % COMPUTES PHI (DEFLECTION ANGLE)
    phi_upper = alphar - theta_upper;
    phi_lower = alphar - theta_lower;

    % STORES AN CATEGORIZES IF ITS COMPRESSION OR EXPANSION.
    waves_upper = cell(size(phi_upper));
    waves_lower = cell(size(phi_lower));
    for u = 1:length(phi_upper)
        if phi_upper(u) > 0
            waves_upper{u} = 'Compression';
        elseif phi_ipper(u) < 0
            waves_upper{u} = 'Expansion';
        else
            waves_upper{u} = 'Neutral';
        end
    end

    for l = 1:length(phi_lower)
        if phi_lower(l) > 0
            waves_lower{l} = 'Compression';
        elseif phi_lower(l) < 0
            waves_lower{l} = 'Expansion';
        else
            waves_lower{l} = 'Neutral;';
        end
    end

    % OSW ANGLES AND PM 
    gamma = 1.4;
    minimum_turning = 0.01;
    betaSW_upper = nan(size(phi_upper));
    betaSW_lower = nan(size(phi_lower));
    PrandtlMeyer_delta_upper = nan(size(phi_upper));
    PrandtlMeyer_delta_lower = nan(size(phi_lower));

    % FOR THE UPPER SURFACE
    for g = 1:length(phi_upper)
        turning = rad2deg(abs(phi_upper(g)));
        % FOR OBLIQUE WAVES
        if phi_upper(g) > 0 && turning > minimum_turning
            betaSW_upper(g) = OSW(Minf, phi_upper(g), gamma);
        
        % FOR EXPANSION WAVES
        elseif phi_upper(g) < 0 && turning > minimum_turning

            PrandtlMeyer_delta_upper(g) = abs(PM_Expansion(AoA_Expansion(Minf, ...
                phi_upper(g), gamma)));
        end
    end

    % FOR THE LOWER SURFACE
    for h = 1:length(phi_lower)
        turning = rad2deg(abs(phi_lower(g)));
        % FOR OBLIQUE WAVES
        if phi_lower(h) > 0 && turning > minimum_turning
            betaSW_lower(h) = OSW(Minf, phi_lower(h), gamma);
        
        % FOR EXPANSION WAVES
        elseif phi_lower(h) < 0 && turning > minimum_turning

            PrandtlMeyer_delta_lower(h) = abs(PM_Expansion(AoA_Expansion(Minf, ...
                phi_lower(h), gamma)));
        end
    end

    % VALID RESULTS STORED FOR OPTIMIZED AIRFOIL SHOCK-EXPANSION METHOD
    results(valid_results).theta_upper = theta_upper;
    results(valid_results).theta_lower = theta_lower;
    results(valid_results).phi_upper = phi_upper;
    results(valid_results).phi_lower = phi_lower;
    results(valid_results).waves_upper = waves_upper;
    results(valid_results).waves_lower = waves_lower;
    results(valid_results).betaSW_upper = betaSW_upper;
    results(valid_results).betaSW_lower = betaSW_lower;
    results(valid_results).PrandtlMeyer_delta_upper = PrandtlMeyer_delta_upper;
    results(valid_results).PrandtlMeyer_delta_lower = PrandtlMeyer_delta_lower;

    
    % Plotting the OPTIMIZED Airfoil using nsu and nsl panels using the
    % results obtained.
    figure(Name = 'Most Optimized Supersonic Airfoil with Highest L/D');
    title('Optimized Supersonic Airfoil Design (Highest Lift-to-Drag Ratio)');
    hold on;

    for k = 1:nsu
        plot([xu_panels(k), xu_panels(k+1)], [yu_optimal(k), yu_optimal(k+1)],'blue');
    end

    for k = 1:nsl
        plot([xl_panels(k), xl_panels(k+1)], [yl_optimal(k), yl_optimal(k+1)],'blue');
    end
    
    plot(xtu_optimal, tu_optimal, '^', 'MarkerFaceColor', 'green');
    text(xtu_optimal, tu_optimal+0.05, 'Max Upper Thickness (t_u)', AffectAutoLimits="on");
    plot(xtl_optimal, -tl_optimal, 'v', 'MarkerFaceColor', 'green');
    text(xtl_optimal, -tl_optimal-0.05, 'Max Lower Thickness (t_l)', AffectAutoLimits="on");
    plot(xu_panels, yu_optimal, 'ro', 'MarkerSize', 3, 'DisplayName', 'Upper Panels');
    plot(xl_panels, yl_optimal, 'ro', 'MarkerSize', 3, 'DisplayName', 'Lower Panels');
    
    xlabel('Chord (x)');
    ylabel('Airfoil Thickness (t)');
    xlim([-0.2*chord, chord*1.2]);
    ylim([-0.5, 0.5]);
    hold off

    
    % Plotting the LEAST OPTIMIZED Airfoil using nsu and nsl panels using the
    % results obtained. (NOT REQUIRED*)
    fprintf('\n|| LEAST OPTIMIZED DESIGN - VALID DESIGN #%d ||\n', minimum_design(1));
    fprintf('Cl = %.5f, Cd = %.5f, Cm = %.5f\n', minimum_design(2), minimum_design(3), minimum_design(4));
    fprintf('LD = %.5f\n', minimum_design(5));
    fprintf('Angle-of-Attack = %.1f\n', minimum_design(6));
    fprintf('tu = %.5f, tl = %.5f\n', minimum_design(7), minimum_design(8));
    fprintf('xtu = %.5f, xtl = %.5f\n', minimum_design(9), minimum_design(10));


    tu_min = minimum_design(7);
    tl_min = minimum_design(8);
    xtu_min = minimum_design(9);
    xtl_min = minimum_design(10);
    xu_min = linspace(0, chord, nsu + 1);
    xl_min = linspace(0, chord, nsl + 1);

    yu_min = ((tu_min ./ xtu_min) .* (xu_panels ./ xtu_min)) .* ... 
         ((((1 - (2 .* xtu_min)).* (xu_panels.^2)) + ...
         (((3 .* (xtu_min.^2)) - 1).* xu_panels) + (xtu_min .* ...
         (2 - (3.*xtu_min)))) ./ ((1 - xtu_min).^2));

    yl_min = (((-tl_min ./ xtl_min) .* (xl_panels ./ xtl_min)) .* ... 
         ((((1 - (2 .* xtl_min)).* (xl_panels.^2)) + ...
         (((3 .* (xtl_min.^2)) - 1).* xl_panels) + (xtl_min .* ...
         (2 - (3.*xtl_min)))) ./ ((1 - xtl_min).^2)));

    figure(Name = 'Supersonic Airfoil with the Lowest L/D');
    title('Supersonic Airfoil with the Lowest Lift-to-Drag Ratio');
    hold on;

    for k = 1:nsu
        plot([xu_panels(k), xu_panels(k+1)], [yu_min(k), yu_min(k+1)],'blue');
    end

    for k = 1:nsl
        plot([xl_panels(k), xl_panels(k+1)], [yl_min(k), yl_min(k+1)],'blue');
    end

    plot(xtu_min, tu_min, '^', 'MarkerFaceColor', 'green');
    text(xtu_min, tu_min+0.05, 'Max Upper Thickness (t_u)', AffectAutoLimits="on");
    plot(xtl_min, -tl_min, 'v', 'MarkerFaceColor', 'green');
    text(xtl_min, -tl_min-0.05, 'Max Lower Thickness (t_l)', AffectAutoLimits="on");
    plot(xu_panels, yu_min, 'ro', 'MarkerSize', 3, 'DisplayName', 'Upper Panels');
    plot(xl_panels, yl_min, 'ro', 'MarkerSize', 3, 'DisplayName', 'Lower Panels');

    xlabel('Chord (x)');
    ylabel('Airfoil Thickness (t)');
    xlim([-0.2*chord, chord*1.2]);
    ylim([-0.5, 0.5]);
    hold off
    else
        fprintf('Results are Invalid.\n');
end

% FUNCTIONS USED FOR OSW AND PRANDTL-MEYER EXPANSION 
function beta_s = OSW(M, theta, gamma)
    theta = abs(theta);
    if theta <= 0
        beta_s = NaN;
        return
    end

    % EQN. 9.23
    f = @(b) tan(theta) - 2.*cot(b).*(M.^2.*sin(b).^2-1)./(M.^2.*(gamma + cos(2.*b)) + 2);

    b_lower = asin(1/M) + 1e-6;
    b_higher = pi/2 - 1e-6;

    if f(b_lower)*f(b_higher) > 0
        beta_s = NaN;
        return
    end

    beta_s = fzero(f, [b_lower b_higher]);
end


function PMF = PM_Expansion(M, gamma)
    if M <= 1
        PMF = 0;
        return
    end
    
    PMF = sqrt((gamma+1)/(gamma-1)).*atan(sqrt(((gamma-1)/gamma+1).*(M^2-1))) - atan(sqrt(M^2-1));

end