% Lifting Line code in MATLAB
% Coded by L. sankar, December 1998
%
croot=input('Please input root chord   ');
ctip=input('Please input tip chord   ');
span=input('Please input span   ');
thetaroot=input('Please input Root twist angle in degrees   ');
thetatip=input('Please input tip twist angle in degrees   ');
a0root=input('Please input root lift curve slope in units/ radian   ');
a0tip=input('Please input lift curve slope at the tip, in units/radian   ');
alpha=input('Please input angle of attack, in degrees   ');
alpha0root=input('Please input zero-lift angle at the root   ');
alpha0tip=input('Please input zero lift angle at the tip   ');
thetaroot=thetaroot * atan(1.)/45.;
thetatip=thetatip * atan(1.)/45.;
alpha = alpha * atan(1.)/45.;
alpha0root=alpha0root * atan(1.)/45.;
alpha0tip=alpha0tip * atan(1.)/45.;
n = 21;
theta=zeros(1,n);
y=zeros(1,n);
c=zeros(1,n);
alp=zeros(1,n);
a=zeros(1,n);
rhs=zeros(n,1);
b=zeros(n,n);
a=zeros(n,1);
%
% Define properties at n span stations
%
pi = 4. * atan(1.);
for i=1:n
  theta(i) = (i) * pi/(2. * n);
  y(i) = span * 0.5 * cos(theta(i));
  c(i) = croot+(ctip-croot)*y(i)*2./span;
  alp(i) = alpha+(alpha0tip-alpha0root-thetaroot+thetatip)*y(i)*2./span;
  a(i) = a0root+(a0tip-a0root)*y(i)*2./span;
 end
pi = 4. * atan(1.);
% Set up 2n x 2n system of equations for A1, A3 , ... A2n-1
for j=1:n
rhs(j,1)=alp(j);
for i=1:n
l = 2 * i-1;
b(j,i)=sin(l*theta(j))*2.*span/pi/c(j)+l * sin(l*theta(j))/sin(theta(j));
end
end
%
% Solve for the Fourier Coefficients
%
a=b\rhs
% Compute wing area and aspect ratio
S=(croot+ctip)/2.*span;
AR=span*span/S;
% Compute CL
CL=a(1)*pi*AR
% Compute CD
CD0=1.;
for i=2:n
CD0=CD0+(2.*i-1)*a(i)*a(i)/(a(1)*a(1));
end
CD = pi * AR * a(1) * a(1) * CD0
% Compute spanwise load distribution, normalized by Freestream velocity times Span
gamma=zeros(1,n);
for i=1:n
gamma(i)=0.0;
y(i)=y(i)*2./span;
for j=1:n
gamma(i)=gamma(i)+2.*a(j)*sin((2*j-1)*theta(i));
end
end
%
% Plot the circulation distribution from root to tip with labels
%
plot(y,gamma);
xlabel('y/Semi-Span');
ylabel('GAMMA/VINF/SPAN');
title('Circulation distribution');