%Laboratory of Aerodynamic Design of Air Vehicles
%Developer: Karpouzas Nikolaos

clc
clear
close all

warning('off','all')

delete 'HorizontalAxisData.xlsx'
delete 'VerticalAxisData.xlsx'

%% Specify the range of the axes of the plot

domX=input('Enter the range of the x-axis of the plot [xmin,xmax]: ');
domY=input('Enter the range of the y-axis of the plot [ymin,ymax]: ');

x=min(domX):(max(domX)-min(domX))/2/100:max(domX);
y=min(domY):(max(domY)-min(domY))/2/100:max(domY);

numX=length(x);
numY=length(y);

[X,Y]=meshgrid(x,y);

u=zeros([numX numY]);
v=zeros([numX numY]);

count=1;


%% Specification of the vortices parameters 

while(1)
    
    fprintf("Specify Vortex No.%d [Gamma,x,y] or enter 0 to stop: \n",count);
    user=input(' ');
    
    if(user(1)==0)
        
        break;
        
    end
    
    gamma(count)=user(1);
    XV(count)=user(2);
    YV(count)=user(3);
   
    count=count+1;
    
end



%% Velocity calculation

for i=1:numX
    
    for j=1:numY
        
        for k=1:length(gamma)
        
         [vx,vy]=vortex(gamma(k),X(i,j),Y(i,j),XV(k),YV(k));
         
         u(i,j)=u(i,j)+vx;
         v(i,j)=v(i,j)+vy;
         
        end
        
    end
    
end


%% Calculation of the x-axis and y-axis data

HorizontalAxisData=[X(Y==0) Y(Y==0) u(Y==0) v(Y==0)]; %A table of [X  Y u v] at y=0

VerticalAxisData=[X(X==0) Y(X==0) u(X==0) v(X==0)];   %A table of [X  Y u v] at x=0

B={'x','y','u','v'};


xlswrite('HorizontalAxisData.xlsx',B,'HorizontalAxisData','A1')
xlswrite('HorizontalAxisData.xlsx',HorizontalAxisData,'HorizontalAxisData','A2')

xlswrite('VerticalAxisData.xlsx',{'x','y','u','v'},'VerticalAxisData','A1')
xlswrite('VerticalAxisData.xlsx',VerticalAxisData,'VerticalAxisData','A2')


%% Streamlines plot
figure
s=streamslice(X,Y,u,v);
set(s,'Linewidth',2);
axis equal;
grid on;
xlabel('X - Axis');
ylabel('Y - Axis');
set(gca,'Fontsize',16);


%% Function for calculation of the vortex velocity

function [u,v] = vortex(gamma,XP,YP,XV,YV)

dx=XP-XV;
dy=YP-YV;

r=sqrt(dx^2+dy^2);

u=gamma/(2*pi)*dy/r^2;
v=-gamma/(2*pi)*dx/r^2;

end
