% Plot the figures of the streamlines, the equipotentials & the pressure field
% for the potential flow around a disk as discussed at :
% https://apmr.matelys.com/BasicsMechanics/Fluid/FlowAroundADisk.html
% luc.jaouen@matelys.com in Acoustical Porous Media Recipes (ISSN 2606-4138)
% Copyleft CC BY 3.0, 2021 
%
% Enter the 21st century, cite webpage that desserves it.


clear all
close all

%% Constants
U = 20;  % Velocity far from disk/cylinder
R = 4;   % cylinder radius

%% Building the "mesh" used to compute and display the results
x_min = -20;  % min x mesh grid results
x_max =  20;  % max x 
y_min = -20;  % min y
y_max =  20;  % max y

n = 800;  % Number of intervals in grid

% Building the mesh on which results will be plotted
[x,y] = meshgrid([x_min:(x_max-x_min)/n:x_max],[y_min:(y_max-y_min)/n:y_max]');

% We set x and y to 0 for points inside the cylinder
% thus r = sqrt(x.^2+y.^2) will be 0 and 1/r : NaN 
% and NaN values are not plotted (they remain white).
for i = 1:length(x)
    for k = 1:length(x) % mesh is square
        f = sqrt(x(i,k).^2 + y(i,k).^2);
        if f < R
            x(i,k) = 0;
            y(i,k) = 0;
        end
    end
end

%% Computation of the polar variables
r = sqrt(x.^2+y.^2);                            
sin_th = y./r;
cos_th = x./r;

%% Expression of the velocity potential function
phi = U*r.*(1+R^2./r.^2).*cos_th;

%% Expression of the streamline function 
psi = U*r.*(1-R^2./r.^2).*sin_th;

%% Computation & plot of the "normalized pressure"
cos_2th = cos_th.^2 - sin_th.^2;
Twop_pinf2_div_rhoU2 = 2*R^2./r.^2.*cos_2th - R^4./r.^4;

figure(1)
colormap jet
surf(x/R,y/R,Twop_pinf2_div_rhoU2,'EdgeColor','None');
view(2);  
grid off
hold on
colorbar
xlabel('x/R')
ylabel('y/R')
title('Normalized pressure')

%% Trick to show a disk at (0,0) with radius R
 % We'll use this in the next plots.
m = 100; 
s = ones(1,m+1)*R;
t = [0:2*pi/m:2*pi];

%% Plot of velocity potentials
 % We'll use the contour function.
 % which is based on the Marching squares algorithm.
 % An alternative is to use Conformal mapping (see end of script).
figure(2)
contour(x,y,phi,50); % 50 is the number of velocity potential curves
hold on 
polar(t,s,'k-');
axis equal
title('Velocity potentials');
grid off

%% Plot of streamlines
 % Again, we'll use the contour function.
 % An alternative is to use Conformal mapping (see end of script).
figure(3)
contour(x,y,psi,50);
hold on 
polar(t,s,'k-');
axis equal
title('Streamlines');
grid off

 
%% Conformal mapping - for streamlines

% The "mesh" is changed for plotting streamlines
% 30 streamlines with 30 points along them
n = 30;
m = 30;

figure(4)
hold on
title("z = f(-1)(w=u+iv)")
xlabel("x");
ylabel("y");

for b=1:1:m
    v = y_min + (b-1)*(y_max-y_min)/(m-1);
    
    x = zeros(1,n);
    y = zeros(1,n);
        
    for a = 1:1:n
        u = x_min + (a-1)*(x_max-x_min)/(n-1);
        
        w = u+1j*v;
        
        % At this point, we have created a rectangular map
        % w = u + iv
        
        signe = -1; 
        if (u > 0) signe = 1; end 
  
        z = 0.5*(w + signe* sqrt(w.^2 -4*R^2)); 
        % Above is f^(-1).
        % Here, we only care about the position of the points but not
        % the exact values of the streamlines function so U is taken to 1.
        
        x(a) = real(z);
        y(a) = imag(z);
    end
    plot(x/R,y/R,'k-')
      
    % Mapping the interior of the disk 
    % (just for fun as no physical meaning)
       for a = 1:1:n
        u = x_min + (a-1)*(x_max-x_min)/(n-1);
        
        w = u+1j*v;
        
        signe = 1;
        if (u > 0) signe = -1; end
        z = 0.5*(w + signe* sqrt(w.^2 -4*R^2));
        
        x(a) = real(z);
        y(a) = imag(z);
        
    end
    plot(x/R,y/R,'k-')
    
end

polar(t,s/R,'r-');

axis equal;

