function lmethod

% function for compute the lift and drag forces on a trapezoidal wing using
% lattice method, ref. 'LOW - SPEED AERODYNAMICS, From Wing Theory to Panel
% Methods' J. Katz, A. Plotkin

rt=input('root(m)=');
t=input('tip(m)=');
s=input('span(m)=');
nrt=input('root elements number=');
ns=input('span elements number=');
v_inf=input('free stream velocity(m/s)=');
alpha=input('angle of attack(deg)=');
density=input('air density(kg/m^3)=');

% identify the meshed wing (3/4-chord locations) matrix %

w=s/ns;

for i=1:1:ns
    
    red(i)=(rt-t)/ns*(i-1)+(rt-t)/(2*ns);
    h(i)=(rt-red(i))/nrt;
    
    for j=1:1:nrt
        
        x(i,j)=-w/2+w*(i);
        y(i,j)=rt-red(i)-(j- 1)*h(i)-.75*h(i);
                  
   end
end
plot(x,y,'*');

hold on

line([0 0 s s],[0 rt t 0],'LineWidth',4)

title('the wing shape with vortex 3/4-chord locations') 
xlabel('span(meter)')
ylabel('root(meter)')

% calculate induced velocities for vortex powers equal to 1 %

for i1=1:1:ns
    for j1=1:1:nrt
        
        % set boundary conditions %
        
        v(1,j1+nrt*(i1-1))=-v_inf*sin(alpha*pi/180);
        
        for i2=1:1:ns
            for j2=1:1:nrt
                r1x=x(i2,j2)-w/2;
                r1y=y(i2,j2)+h(i2)/2;
                r1z=0;
                r2x=x(i2,j2)+w/2;
                r2y=y(i2,j2)+h(i2)/2;
                r2z=0;
                rx=x(i1,j1);
                ry=y(i1,j1);                
                rz=0;
                r1=[r1x r1y r1z];
                r2=[r2x r2y r2z];
                r=[rx ry rz];
                r2p=-r+r2;
                r1p=-r+r1;
                            
                ind=horseshoe(r1,r2,r);
                induced_velocity(j1+nrt*(i1-1),j2+nrt*(i2-1))=ind(1,3);
            end
        end
    end
end

% obtain vortex powers , lift and drag %

vortex_power=v*inv(induced_velocity);

delta_lift=0;

for i3=1:1:nrt*ns
    delta_lift=density*v_inf*vortex_power(i3)*w+delta_lift;
end

lift=delta_lift
drag=lift*sin(alpha*pi/180)

% subfunction for calculation of induced velocity at location r due to a horseshoe vortex at (r1,r2) %

function k=horseshoe(r1,r2,r)

vp=1;
n=[0 1 0];
l=r2-r1;
r2p=-r+r2;
r1p=-r+r1;

s=cross(r1p,r2p);

if norm(s)==0
   vtip=0;
else
r12=r1p/norm(r1p)-r2p/norm(r2p);
vtip=-(vp/(4*pi))*s/(norm(s))^2*(dot(l,r12));
end

h1=norm(cross(n,r1p));

if h1==0
    vl1=0;
else
s1=-(cross(n,r1p))/norm(r1p);
e1=s1/norm(s1);
vl1=(vp/(4*pi*h1))*(1+dot(n,r1p)/norm(r1p))*e1;
end

h2=norm(cross(-n,r2p));

if h2==0
    vl2=0
else
s2=-(cross(-n,r2p))/norm(r2p);
e2=s2/norm(s2);
vl2=(vp/(4*pi*h2))*(1-dot(-n,r2p)/norm(r2p))*e2;
end

k=vtip+vl1+vl2;