%Copyright (C) 2004  Mihai Pruna, Alberto Davila

%This program is free software; you can redistribute it and/or
%modify it under the terms of the GNU General Public License
%as published by the Free Software Foundation; either version 2
%of the License, or (at your option) any later version.

%This program is distributed in the hope that it will be useful,
%but WITHOUT ANY WARRANTY; without even the implied warranty of
%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%GNU General Public License for more details.
%http://www.gnu.org/licenses/gpl.html

%You should have received a copy of the GNU General Public License
%along with this program; if not, write to the Free Software
%Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.


%3d panel method for calculating lift of a wing with a cambered airfoil
taper=1;
AR=6;
%semispan of the wing
span=3;
sweep=0;
%# of spanwise panels
n=10;
%# of points of the airfoil
m=6;
% a camberline can be defined
x=zeros(m,1);
z=zeros(m,1);
x(1)=0;
z(1)=0;
x(2)=0.2;
z(2)=0.0;
x(3)=0.4;
z(3)=0.0;
x(4)=0.6
z(4)=0.0;
x(5)=0.8;
z(5)=0.0;
x(6)=1;
z(6)=0;
%Angle of Attack
AOA=5;
%freestream velocity magnitude
Q=1
%freestream velocity vector
Qinf=[Q.*cos(AOA.*3.1416/180),0,Q.*sin(AOA*3.1416/180)];

A=zeros((m-1)*n,(m-1)*n);
Xpanelgrid=zeros(20,20);
Ypanelgrid=zeros(20,20);
Zpanelgrid=zeros(20,20);

%generate grid corner coordinates
y=zeros(20,1);

for j=1:n
    for i=1:m
    Ypanelgrid(i,j)=j./n.*span;
    Zpanelgrid(i,j)=z(i)-z(i).*(1-taper)*Ypanelgrid(i,j)./span;        
    Xpanelgrid(i,j)=Ypanelgrid(i,j).*tan(sweep.*3.1416/180)+x(i)-x(i).*(1-taper)*Ypanelgrid(i,j)./span;
end
end

%generate vortex ring corner coordinates
VRPCX=zeros(20,20);
VRPCY=zeros(20,20);
VRPCZ=zeros(20,20);

for i=1:m-1
    VRPCY(i,1)=0;
    VRPCX(i,1)=x(i)+(x(i+1)-x(i))/4;
    VRPCZ(i,1)=z(i)+(z(i+1)-z(i))/4;
end
VRPCY(m,1)=0;
VRPCX(m,1)=x(m)+(x(m)-x(m-1))/4;
VRPCZ(m,1)=z(m)+(z(m)-z(m-1))/4;

for j=2:n+1
    for i=1:m-1
        VRPCY(i,j)=Ypanelgrid(i,j-1);
        VRPCX(i,j)=Xpanelgrid(i,j-1)+(Xpanelgrid(i+1,j-1)-Xpanelgrid(i,j-1))/4;
        VRPCZ(i,j)=Zpanelgrid(i,j-1)+(Zpanelgrid(i+1,j-1)-Zpanelgrid(i,j-1))/4;
    end
    VRPCY(m,j)=Ypanelgrid(m,j-1);
        VRPCX(m,j)=Xpanelgrid(m,j-1)+(Xpanelgrid(m,j-1)-Xpanelgrid(m-1,j-1))/4;
        VRPCZ(m,j)=Zpanelgrid(m,j-1)+(Zpanelgrid(m,j-1)-Zpanelgrid(m-1,j-1))/4;
    end
    %adding wake points
    for j=1:n
        VRPCX(m+1,j)=20*span;
        VRPCY(m+1,j)=VRPCY(m,j);
        VRPCZ(m+1,j)=0;
    end
    
        %generate collocation points
        
        COLOCX=zeros(20,20);
        COLOCY=zeros(20,20);
        COLOCZ=zeros(20,20);
        
        for j=1:n
            for i=1:m-1
                COLOCX(i,j)=(VRPCX(i,j)+VRPCX(i+1,j)+VRPCX(i,j+1)+VRPCX(i+1,j+1))/4;
                COLOCY(i,j)=(VRPCY(i,j)+VRPCY(i+1,j)+VRPCY(i,j+1)+VRPCY(i+1,j+1))/4;
                COLOCZ(i,j)=(VRPCZ(i,j)+VRPCZ(i+1,j)+VRPCZ(i,j+1)+VRPCZ(i+1,j+1))/4;
        
            end
        end
 %initializing induced velocities    
velinduced=zeros(3);
velinducedotherhalf=zeros(3);
totalvelinduced=zeros(3,1);
k=0;
l=0;

%computing the influence coefficient matrix
for a=1:m-1
    for b=1:n
        k=k+1;
l=0;
        for i=1:m-1
        for j=1:n
l=l+1;
velinduced=vortxl(COLOCX(a,b),COLOCY(a,b),COLOCZ(a,b),VRPCX(i,j),VRPCY(i,j),VRPCZ(i,j), VRPCX(i,j+1),VRPCY(i,j+1),VRPCZ(i,j+1),1)+vortxl(COLOCX(a,b),COLOCY(a,b),COLOCZ(a,b),VRPCX(i,j+1),VRPCY(i,j+1),VRPCZ(i,j+1), VRPCX(i+1,j+1),VRPCY(i+1,j+1),VRPCZ(i+1,j+1),1)+vortxl(COLOCX(a,b),COLOCY(a,b),COLOCZ(a,b),VRPCX(i+1,j+1),VRPCY(i+1,j+1),VRPCZ(i+1,j+1), VRPCX(i+1,j),VRPCY(i+1,j),VRPCZ(i+1,j),1)+vortxl(COLOCX(a,b),COLOCY(a,b),COLOCZ(a,b),VRPCX(i+1,j),VRPCY(i+1,j),VRPCZ(i+1,j), VRPCX(i,j),VRPCY(i,j),VRPCZ(i,j),1);
velinducedotherhalf=vortxl(COLOCX(a,b),-COLOCY(a,b),COLOCZ(a,b),VRPCX(i,j),VRPCY(i,j),VRPCZ(i,j), VRPCX(i,j+1),VRPCY(i,j+1),VRPCZ(i,j+1),1)+vortxl(COLOCX(a,b),-COLOCY(a,b),COLOCZ(a,b),VRPCX(i,j+1),VRPCY(i,j+1),VRPCZ(i,j+1), VRPCX(i+1,j+1),VRPCY(i+1,j+1),VRPCZ(i+1,j+1),1)+vortxl(COLOCX(a,b),-COLOCY(a,b),COLOCZ(a,b),VRPCX(i+1,j+1),VRPCY(i+1,j+1),VRPCZ(i+1,j+1), VRPCX(i+1,j),VRPCY(i+1,j),VRPCZ(i+1,j),1)+vortxl(COLOCX(a,b),-COLOCY(a,b),COLOCZ(a,b),VRPCX(i+1,j),VRPCY(i+1,j),VRPCZ(i+1,j), VRPCX(i,j),VRPCY(i,j),VRPCZ(i,j),1);
totalvelinduced(1)=velinduced(1)+velinducedotherhalf(1);
totalvelinduced(2)=velinduced(2)-velinducedotherhalf(2);
totalvelinduced(3)=velinduced(3)+velinducedotherhalf(3);
if (i==m-1)
wakevelinduced=vortxl(COLOCX(a,b),COLOCY(a,b),COLOCZ(a,b),VRPCX(i+1,j),VRPCY(i+1,j),VRPCZ(i+1,j), VRPCX(i+1,j+1),VRPCY(i+1,j+1),VRPCZ(i+1,j+1),1)+vortxl(COLOCX(a,b),COLOCY(a,b),COLOCZ(a,b),VRPCX(i+1,j+1),VRPCY(i,j+1),VRPCZ(i+1,j+1), VRPCX(i+2,j+1),VRPCY(i+2,j+1),VRPCZ(i+2,j+1),1)+vortxl(COLOCX(a,b),COLOCY(a,b),COLOCZ(a,b),VRPCX(i+2,j+1),VRPCY(i+2,j+1),VRPCZ(i+2,j+1), VRPCX(i+2,j),VRPCY(i+2,j),VRPCZ(i+2,j),1)+vortxl(COLOCX(a,b),COLOCY(a,b),COLOCZ(a,b),VRPCX(i+2,j),VRPCY(i+2,j),VRPCZ(i+2,j), VRPCX(i+1,j),VRPCY(i+1,j),VRPCZ(i+1,j),1);
wakevelinducedotherhalf=vortxl(COLOCX(a,b),-COLOCY(a,b),COLOCZ(a,b),VRPCX(i+1,j),VRPCY(i+1,j),VRPCZ(i+1,j), VRPCX(i+1,j+1),VRPCY(i+1,j+1),VRPCZ(i+1,j+1),1)+vortxl(COLOCX(a,b),-COLOCY(a,b),COLOCZ(a,b),VRPCX(i+1,j+1),VRPCY(i,j+1),VRPCZ(i+1,j+1), VRPCX(i+2,j+1),VRPCY(i+2,j+1),VRPCZ(i+2,j+1),1)+vortxl(COLOCX(a,b),-COLOCY(a,b),COLOCZ(a,b),VRPCX(i+2,j+1),VRPCY(i+2,j+1),VRPCZ(i+2,j+1), VRPCX(i+2,j),VRPCY(i+2,j),VRPCZ(i+2,j),1)+vortxl(COLOCX(a,b),-COLOCY(a,b),COLOCZ(a,b),VRPCX(i+2,j),VRPCY(i+2,j),VRPCZ(i+2,j), VRPCX(i+1,j),VRPCY(i+1,j),VRPCZ(i+1,j),1);
totalvelinduced(1)=totalvelinduced(1)+wakevelinduced(1)+wakevelinducedotherhalf(1);
totalvelinduced(2)=totalvelinduced(2)+wakevelinduced(2)-wakevelinducedotherhalf(2);
totalvelinduced(3)=totalvelinduced(3)+wakevelinduced(3)+wakevelinducedotherhalf(3);
end
    

cp=zeros(3,1);
normal=zeros(3,1);
cp=cross([VRPCX(a,b+1)-VRPCX(a+1,b),VRPCY(a,b+1)-VRPCY(a+1,b),VRPCZ(a,b+1)-VRPCZ(a+1,b)],[VRPCX(a+1,b+1)-VRPCX(a,b),VRPCY(a+1,b+1)-VRPCY(a,b),VRPCZ(a+1,b+1)-VRPCZ(a,b)]);
cpmag=sqrt(cp(1).*cp(1)+cp(2).*cp(2)+cp(3).*cp(3));
normal=cp/cpmag;
A(k,l)=dot(totalvelinduced,normal);
end
end
end
end
k=0;
RHS=zeros((m-1)*n,1);
GAMMA=zeros((m-1)*n,1);
%computing the right hand side
for a=1:m-1
    for b=1:n
    k=k+1;
    cp=cross([VRPCX(a,b+1)-VRPCX(a+1,b),VRPCY(a,b+1)-VRPCY(a+1,b),VRPCZ(a,b+1)-VRPCZ(a+1,b)],[VRPCX(a+1,b+1)-VRPCX(a,b),VRPCY(a+1,b+1)-VRPCY(a,b),VRPCZ(a+1,b+1)-VRPCZ(a,b)]);
cpmag=sqrt(cp(1).*cp(1)+cp(2).*cp(2)+cp(3).*cp(3));
normal=cp/cpmag;
RHS(k)=-dot(Qinf,normal);
end
end
GAMMA=inv(A)*RHS;
%now we calculate lift loading distribution for each wing section
GAMMA
for j=1:n
L(j)=GAMMA(j);
    for i=2:m-1
        L(j)=L(j)+GAMMA(j+(i-1)*n)-GAMMA(j+(i-1)*n-n);
    end
end
L

 %file vortxl.m
function vortex=vortxl(x,y,z,x1,y1,z1,x2,y2,z2,gamma)
rcut=1*10^(-10);
r1r2x=(y-y1)*(z-z2)-(z-z1)*(y-y2);
r1r2y=-((x-x1)*(z-z2)-(z-z1)*(x-x2));
r1r2z=(x-x1)*(y-y2)-(y-y1)*(x-x2);
square=r1r2x.*r1r2x+r1r2y.*r1r2y+r1r2z.*r1r2z;
r1=sqrt((x-x1)*(x-x1)+(y-y1)*(y-y1)+(z-z1)*(z-z1));
r2=sqrt((x-x2)*(x-x2)+(y-y2)*(y-y2)+(z-z2)*(z-z2));
if (r1<rcut) || (r2<rcut) || (square<rcut)
    vortex(1)=0;
    vortex(2)=0;
    vortex(3)=0;
else
    ror1=(x2-x1)*(x-x1)+(y2-y1)*(y-y1)+(z2-z1)*(z-z1);
    ror2=(x2-x1)*(x-x2)+(y2-y1)*(y-y2)+(z2-z1)*(z-z2);
    coeff=gamma./(4.0*3.1416*square)*(ror1./r1-ror2./r2);
    vortex(1)=r1r2x*coeff;
    vortex(2)=r1r2y*coeff;
    vortex(3)=r1r2z*coeff;
end