# MATLAB: Creating a Code to determine Lift Curve Slope

aerodynamicslifting surfaces

For an Aerodynamics couse, we need to create a code that computes lifting surfaces (wing) aerodynamics. I cannot get the graph to display a nice CL vs AoA curve, but only a single value. I don't know how I can accomplish this. As I am a novice to MAtlab, I would appreciate some pointers or hints. Thank you.
clc;
% Default airfoil or not…
Default_AF = 0;
% A/F Characteristics:
Cl_alpha = 5.79;
% Sectional lift curve slope [deg]:
alpha_0 = -2.0;
Cm_alpha = 0;
% 2D Pitching moment:
Cm_0 = -0.025;
% Coefficient in series expansion of circulation distribution [deg]:
alpha_nl = 8;
% Wing planform Geometry [inches]:
% [root-LE,root-TE,tip-LE,tip-TE]
x1 = 0;
y1 = 0;
x2 = 200;
y2 = 0;
x3 = 0;
y3 = 1000;
x4 = 200;
y4 = 1000;
% Wing twist angle [deg]:
theta_t = 0;
% Number of spanwise coefficients:
N = 100;
% Flight Conditions:
% Free-stream Velocity [ft/sec]
V = 100;
% Density [slugs/ft^3]:
rho = 0.002344;
% Viscous Force [lb.s/ft^2]
mu = 0.0000003737;
% Angle of Attack [deg]
alpha = 5;
% Wing span:
b = sqrt((x2-x1)^2 + (y2-y1)^2);
% Wing chord:
c = sqrt((x3-x1)^2 + (y3-y1)^2);
% Wing area:
S = b*c;
% Wing semi-span:
s = b/2;
% Aspect Ratio:
AR = b^2./S;
% Applying segment length of the wings:
alpha_segment = zeros(1,N);
theta_segment = zeros(1,N);
% Taper Ratio and its range within the wing (which is 0 to 1):
for lambda_r = 0:0.25:1
% Root Chord:
c_root = (2*S) / ((1+lambda_r)*b);
% Tip Chord:
c_tip = ((2*S)/((1+lambda_r)*b)) * (1-((2*(1-lambda_r))/b)*(b/2));
% Mean Aerodynamic Chord:
MAC = (2/3)*(c_root+c_tip-(c_root*c_tip)/(c_root+c_tip));
% Taper ratio:
lambda = c_tip / c_root;
% Providing segment limitation within range:
alpha_segment(1) = alpha;
theta_segment(1) = pi/2;
% Completing the Lifting Line Theory to determine the
% matrix value of the segments:
for N = 2:N
c(N) = c_root+(c_tip-c_root)/N.*N;
alpha_segment(N) = alpha+(theta_t)/N.*N;
theta_segment = pi/(2*N):pi/(2*N):pi/2;
end
% Lifting line theory construction:
xb = 4*b./(pi*c);
% Completing the [B] matrix:
for i = 1:N
for j = 1:N
B(i,j) = (sin(j.*theta_segment(i))).*(xb(i)+j/(sin(theta_segment(i))));
end
slope(i) = (alpha_segment(i)-alpha_0)*(pi/180);
end
A = B\transpose(slope);
end
for fai = 1:N
delta_LE(fai) = fai*(A(fai)/A(1))^2; %#ok<*SAGROW>
end
% Section Lift Coefficient of Airfoil
cl = 0.5*cos(pi/b);
% Wing Lift Coefficient:
CL = pi*AR*A(1);
% Span Efficiency:
delta = sum(delta_LE);
CD_0 = 1/(1+delta);
% Induced drag coefficient:
CD_i = CL.^2 / (pi*CD_0*AR);
% Speed of sound (assuming 20 degree dry air) [ft/sec]:
C = 1125.33;
% Mach Number:
M = V / C;
% Dynamic Pressure:
q = (0.5*rho*V^2);
% Pitching Moment Coefficient:
CM = M / q*S*c;
% Subsonic Lift:
alpha_inf = 2*pi*cos(delta_LE);
% 3D Lift Curve Slope of airfoil:
CL_alpha = 2*pi*A;
% Geometric Angle of Attack of Wing:
alpha = (alpha_inf)/(1+(alpha_inf)/(pi*AR));
%———————————————————-
% Wing lift coefficient (CL) vs Angle of Attack (alpha):
figure(1);
plot(alpha,CL,'-o')
grid on
title('Wing lift coefficient (CL) vs Angle of Attack (alpha)')
ylabel('Wing Lift Coefficient, CL')
xlabel('Angle of Attack, alpha [deg]')
% Pitching moment coefficient(CM) vs Angle of Attack (alpha):
figure(2);
plot(alpha,CM,'-o')
grid on
title('Pitching moment coefficient(CM) vs Angle of Attack (alpha)')
ylabel('Pitching Moment Coefficient, CM')
xlabel('Angle of Attack, alpha [deg]')
% Wing lift coefficient (CL) vs Induced drag coefficient (CD_i):
figure(3);
plot(CD_i,CL,'-o')
grid on
title('Wing lift coefficient (CL) vs Induced drag coefficient (CD_i)')
ylabel('Wing Lift Coefficient, CL')
xlabel('Induced Drag Coefficient, CD_i')

lambda_rv = 0:0.25:1;for k = 1:numel(lambda_rv)    lambda_r = lambda_rv(k);    % Root Chord:    c_root = (2*S) / ((1+lambda_r)*b);    % Tip Chord:    c_tip = ((2*S)/((1+lambda_r)*b)) * (1-((2*(1-lambda_r))/b)*(b/2));    % Mean Aerodynamic Chord:    MAC = (2/3)*(c_root+c_tip-(c_root*c_tip)/(c_root+c_tip));    % Taper ratio:    lambda = c_tip / c_root;    % Providing segment limitation within range:    alpha_segment(1) = alpha;    theta_segment(1) = pi/2;    % Completing the Lifting Line Theory to determine the    % matrix value of the segments:    for N = 2:N                c(N) = c_root+(c_tip-c_root)/N.*N;        alpha_segment(N) = alpha+(theta_t)/N.*N;        theta_segment = pi/(2*N):pi/(2*N):pi/2;            end    % Lifting line theory construction:    xb = 4*b./(pi*c);    % Completing the [B] matrix:    for i = 1:N                for j = 1:N                        B(i,j) = (sin(j.*theta_segment(i))).*(xb(i)+j/(sin(theta_segment(i))));                    end                slope(i) = (alpha_segment(i)-alpha_0)*(pi/180);            end    A(:,k) = B\transpose(slope);    end
fai = 1:N;delta_LE = fai(:).*(A/A(1)).^2; %#ok<*SAGROW>
alpha = (alpha_inf)./(1+(alpha_inf)/(pi*AR));