Creating a Code to determine Lift Curve Slope

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:
% Lift coefficient [1/rad]:
Cl_alpha = 5.79;
Cl_alphaInRadians = deg2rad(Cl_alpha);

% Sectional lift curve slope [deg]:
alpha_0 = -2.0;

% Pitching moment [1/rad]:
Cm_alpha = 0;
Cm_alphaInRadians = deg2rad(Cm_alpha);

% 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 = Btranspose(slope);

end

% Leading edge sweep angle
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’)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:
% Lift coefficient [1/rad]:
Cl_alpha = 5.79;
Cl_alphaInRadians = deg2rad(Cl_alpha);

% Sectional lift curve slope [deg]:
alpha_0 = -2.0;

% Pitching moment [1/rad]:
Cm_alpha = 0;
Cm_alphaInRadians = deg2rad(Cm_alpha);

% 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 = Btranspose(slope);

end

% Leading edge sweep angle
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 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’) 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:
% Lift coefficient [1/rad]:
Cl_alpha = 5.79;
Cl_alphaInRadians = deg2rad(Cl_alpha);

% Sectional lift curve slope [deg]:
alpha_0 = -2.0;

% Pitching moment [1/rad]:
Cm_alpha = 0;
Cm_alphaInRadians = deg2rad(Cm_alpha);

% 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 = Btranspose(slope);

end

% Leading edge sweep angle
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));

%———————————————————-

Cart

Cart