The most basic approach is given here:
clc
clear
N = 9;
S = 18.1;
AR = 7;
lambda = 0.8;
alpha_twist = -1.5;
i_w = 10;
a_2d = 6.3;
a_0 = -3;
a_0_fd = -6;
b = sqrt(AR*S);
bf_b=0.6;
MAC = S/b;
Croot = (1.5*(1+lambda)*MAC)/(1+lambda+lambda^2);
theta = pi/(2*N):pi/(2*N):pi/2;
alpha=i_w+alpha_twist:-alpha_twist/(N-1):i_w;
alpha_0 = zeros(1,N);
for i=1:N
if (i/N)>(1-bf_b)
alpha_0(i)=a_0_fd;
else
alpha_0(i)=a_0;
end
end
z = (b/2)*cos(theta);
c = Croot * (1 - (1-lambda)*cos(theta));
mu = c * a_2d / (4 * b);
LHS = mu .* (alpha-alpha_0)/57.3;
B = zeros(N,N);
for i=1 : N
for j=1 : N
B(i,j) = sin((2*j-1) * theta(i)) * (1 + (mu(i) *(2*j-1)) / sin(theta(i)));
end
end
A=B\transpose(LHS);
sum1 = zeros(1, N);
sum2 = zeros(1, N);
for i = 1 : N
for j = 1 : N
sum1(i) = sum1(i) + (2*j-1) * A(j)*sin((2*j-1)*theta(i));
sum2(i) = sum2(i) + A(j)*sin((2*j-1)*theta(i));
end
end
CL_TO = pi * AR * A(1);
