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

clc
clear
close all;

%% Input Parameters (Specified by User)

% Free-Stream Velocity
Vinf = 1;

% Source Strength
lambda = 1;

% Distance between sources
d = 1;

%% Rankine Body Dimensions

% First Source Coordinates
x1 = -d/2; 
y1 = 0;

% Second Source Coordinates
x2 = d/2;
y2 = 0;

% Length & Thickness
[L,t] = dims(Vinf,lambda,x1,x2);

% Ellipse Coordinates
theta = linspace(0,2*pi,50);
xc = (x1+x2)/2;
yc = (y1+y2)/2;
xr = L/2*cos(theta)+xc;
yr = t/2*sin(theta)+yc;


%% Computational Grid

% Grid Points
numX = 501;
numY = 501;

% Domain Range
MinRange = -2*max([10,abs(x1),abs(x2)]);
MaxRange = 2*max([10,abs(x1),abs(x2)]);

x = MinRange:(MaxRange-MinRange)/(numX-1):MaxRange;
y = MinRange:(MaxRange-MinRange)/(numY-1):MaxRange;

% 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
        
        % Check if grid points are inside rankine body
        [in,on] = inpolygon(X(i,j),Y(i,j),xr,yr);
      
        if (in == 1 || on == 1)

            vx(i,j) = 0;
            vy(i,j) = 0;
        
        else
            
               dx1 = X(i,j)-x1;
               dy1 = Y(i,j)-y1;

               dx2 = X(i,j)-x2;
               dy2 = Y(i,j)-y2;

               r1 = sqrt(dx1^2+dy1^2);
               r2 = sqrt(dx2^2+dy2^2);

               vx(i,j) = (lambda*dx1)/(2*pi*r1^2) + ...
                   (-lambda*dx2)/(2*pi*r2^2) + Vinf;

               vy(i,j) = (lambda*dy1)/(2*pi*r1^2) + ...
                   (-lambda*dy2)/(2*pi*r2^2);
            
        end
       
       
    end
end

%% Streamlines

numSL  = 50; % # of streamlines
xStart = -2*max([10,abs(x1),abs(x2)]).*ones(numSL,1);
yStart = linspace(-2*max([10,abs(x1),abs(x2)]), ...
    2*max([10,abs(x1),abs(x2)]),numSL);

figure
sl=streamline(X,Y,vx,vy,xStart,yStart);
set(sl,'Color','red','Linewidth',2);
hold on;
fill(xr,yr,'k');
sources = plot([x1 x2],[y1 y2],'ko','MarkerFaceColor','c','Markersize',10);
axis equal;
grid on;
title('Rankine Body');
xlabel('X - Axis');
ylabel('Y - Axis');
legend(sources,'Sources');
set(gca,'Fontsize',16);

% Domain Span
a = 15;
xlim([-a a])
ylim([-(a-5) (a-5)])


%% Velocity Magnitude and Pressure Coefficient Contour Plots

vXY = sqrt(vx.^2+vy.^2);
CpXY = 1-vXY.^2/Vinf^2;

figure                                                                                                         
contourf(X,Y,vXY,100,'EdgeColor','none');
colormap(jet(256));
c = colorbar;
c.Label.String = 'Velocity Magnitude (m/s)';
hold on;
grid on;
axis equal;
fill(xr,yr,'k');
sources = plot([x1 x2],[y1 y2],'ko','MarkerFaceColor','c','Markersize',10);
title('Rankine Body');
xlabel('X - Axis');
ylabel('Y - Axis');
set(gca,'Fontsize',16);
legend(sources,'Sources');

% Domain Span
a = 15;
xlim([-a a])
ylim([-(a-5) (a-5)])


figure                                                                                                         
contourf(X,Y,CpXY,100,'EdgeColor','none');
colormap(jet(256));
c=colorbar;
c.Label.String='Pressure Coefficient Cp';
hold on;
grid on;
axis equal;
fill(xr,yr,'k');
sources=plot([x1 x2],[y1 y2],'ko','MarkerFaceColor','c','Markersize',10);
title('Rankine Body');
xlabel('X - Axis');
ylabel('Y - Axis');
set(gca,'Fontsize',16);
legend(sources,'Sources');

xlim([-a a])
ylim([-(a-5) (a-5)])


% Function that calculates dimensions of Rankine Body 
function [L,t] = dims(Vinf,lambda,x1,x2)

b = abs(x1);

L = 2*sqrt(lambda/(2*pi*Vinf)*(x2-x1)+abs(x1*x2));

syms t
t = vpasolve(t/(2*b)==(lambda/(2*Vinf*b)-lambda/(pi*Vinf*b)*atan(t/(2*b))));
t = double(t);


end