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

clc
clear
close all


%% Input Parameters (Specified by User)

% Free-Stream Velocity
Vinf = 1;

% Dipole Strength
k = 1;

% Vortex Strength
gamma = 0;


%% Cylinder Radius & Coordinates

R = sqrt(k/(2*pi*Vinf));

theta = linspace(0,2*pi,100);
x_cyl = R*cos(theta);
y_cyl = R*sin(theta);


%% Computational Grid

% Grid Points
numX = 501;
numY = 501;

% Domain
x = linspace(-10,10,numX);
y = linspace(-10,10,numY);

% Create Structured Grid
[X,Y] = meshgrid(x,y);

% Initialize Grid Velocities
vx = zeros([numX numY]);
vy = zeros([numX numY]);


%% Grid Velocities Calculation

for i=1:numX
    
    for j=1:numY
        
       r = sqrt(X(i,j)^2+Y(i,j)^2);
     
       vx(i,j) = (k/(2*pi))*(Y(i,j)^2-X(i,j)^2)/r^4 + ...
           (gamma*Y(i,j))/(2*pi*r^2) + Vinf;
       
       vy(i,j) = -(k/(2*pi))*(2*X(i,j)*Y(i,j))/r^4 - ...
           (gamma*X(i,j))/(2*pi*r^2);
       
    end

end


%% Velocity and Pressure Coeffcient Calculations

vXY=sqrt(vx.^2+vy.^2);

for i=1:numX

    for j=1:numY
        
        r = sqrt(X(i,j)^2+Y(i,j)^2);
        
        if(r<R)
            
            vXY(i,j) = 0;
            
        end
        
    end
end

CpXY=1-vXY.^2/Vinf^2;


%% Velocity Magnitude and Pressure Coeffcient Contour Plots

figure                                                                                                            
set(gcf,'Color','White');                                               
set(gca,'FontSize',12);                                                 
contourf(X,Y,CpXY,100,'EdgeColor','none');
hold on; grid on;
c=colorbar;
colormap(jet(256));
c.Label.String = 'Pressure Coefficient (Cp)';
title("Cylinder - \Gamma = "+num2str(gamma));
fill(x_cyl,y_cyl,'k','Linewidth',2);
axis equal;
lim=R+3;
xlim([-lim lim]);
ylim([-lim lim]);
xlabel('X - Axis');
ylabel('Y - Axis');
set(gca,'Fontsize',20);


figure                                                              
cla; hold on; grid on;                                                  
set(gcf,'Color','White');                                               
set(gca,'FontSize',12);                                                 
contourf(X,Y,vXY,100,'EdgeColor','none');
c=colorbar;
c.Label.String="Velocity Magnitude (m/s)";
colormap(jet(256))
fill(x_cyl,y_cyl,'k','Linewidth',2);
axis equal;
lim=R+3;
xlim([-lim lim]);
ylim([-lim lim]);
xlabel('X - Axis');
ylabel('Y - Axis');
title("Cylinder - \Gamma = "+num2str(gamma));
set(gca,'Fontsize',20);
