% Plot the figures of the streamlines for the conformal mappings as discussed at :
% https://apmr.matelys.com/BasicsAcoustics/LengthCorrections/
% 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

figure(1)
hold on
title("z(x+iy) : w=f(z)=(u+iv)")
xlabel("u");
ylabel("v");

Delta = 0.2; % step for the grids
y_vector = [0:0.2:pi,pi];

for b = 1:1:length(y_vector)
    y = y_vector(b);
    m = 1;
    for a=-5:Delta:5
        x = a; 
        z = x+1j*y;
        
        % Our first conformal mapping
        % (upper half of the 2D flanged duct).
        w = (2/pi) * ( sqrt(exp(z)+1) + 1/2 * log( ( sqrt(exp(z)+1)-1) ./ (sqrt(exp(z)+1)+1) ) );     
        u(m) = real(w);
        v(m) = imag(w);
        m = m+1;        
    end   
    figure(1)
    plot(u,v,'-','Color',[0,204,255]/255)
end

% We plot the length correction from the J. W. Strutt's second formula
% (non-constant velocity accros the duct aperture).
plot(0.82,0,'.','Color',[255,0,205]/255,'MarkerSize',20); % trick to plot a solid disk

axis([-1 5 -3 3])
axis square

% Here are the script for the 2nd plot.
clear all

figure(2)
hold on
title("z(x+iy) : w=f(z)=(u+iv)")
xlabel("u");
ylabel("v");

Delta = 0.2;
y_vector = [-pi,-pi:Delta:pi,pi];

for b = 1:1:length(y_vector)
    y = y_vector(b);
    m = 1;
    for a=-5:Delta:5
        x = a; 
        z = x+1j*y;
        
        % Second mapping (the duct opening in the open plane)
        w = 1/pi * (1+exp(z)+z);      
        u(m) = real(w);
        v(m) = imag(w);
        m = m+1;     
    end
    figure(2)
    plot(u,v,'-','Color',[0,204,255]/255); 
end

plot(0.61,0,'.','Color',[255,0,205]/255,'MarkerSize',20); % trick to plot a solid disk
axis([-1 5 -3 3]) 
axis square

% Here start the plot for the animated conformal mapping
clear all

Delta = 0.2; % our grid step
y_vector = [0:0.2:pi,pi]; % a vector for y positions
nb_frame_max = 2*11;

% t is a parametric value (from 0 to 1 and 1 to 0)
% used for the animation (see mapping below)
for nb_frame = 1:1:nb_frame_max
    if (nb_frame <= nb_frame_max / 2)
        t = (nb_frame-1)/(nb_frame_max/2-1);
    else
        t = (nb_frame_max-nb_frame)/(nb_frame_max/2-1);
    end
    
    figure(3)
    hold on
    grid off
    axis off
    set(gcf,'color','w');
    
    for b = 1:1:length(y_vector) % First step : draw streamlines
        y = y_vector(b);
        
        m = 1;
        for a=-5:Delta:5
            x = a;
            z = x+1j*y;
            
            % We'll use the mapping (1-t)*z + t*f(z) to animate from z (regular
            % mesh grid) to f(z) : our function
            
            % The mapping below maps the upper half plane to 
            % the upper half plane : i.e. only half of the aperture
            w = (1-t)*z + t*( 2/pi * ( sqrt(exp(z)+1) + 1/2 * log( ( sqrt(exp(z)+1)-1) ./ (sqrt(exp(z)+1)+1) ) ) );
            u(m) = real(w);
            v(m) = imag(w);
            m = m+1; 
        end
        figure(3)
        plot(u,v,'Color',[0,204,255]/255)  % Draw upper half of the aperture
        plot(u,-v,'Color',[0,204,255]/255) % Draw bottom half
    end
    
    y = [0:Delta:pi,pi];  % Next step : drawing the equipotentials
    for a=-5:Delta:5
        x = ones(1,length(y))*a;
        z = x+1j*y;
        
        % Again, we use the mapping (1-t)*z + t*f(z) to animate from z (regular
        % mesh grid) to f(z) : our function
        w = (1-t)*z + t*( 2/pi * ( sqrt(exp(z)+1) + 1/2 * log( ( sqrt(exp(z)+1)-1) ./ (sqrt(exp(z)+1)+1) ) ) );
        u = real(w);
        v = imag(w);
        figure(3)
        plot(u,v,'Color',[255,102,0]/255)  % Draw upper half of the aperture
        plot(u,-v,'Color',[255,102,0]/255) % Draw bottom half
    end
    axis([-1 15 -8 8])
    
    % Save animation steps to a gif file
    frame = getframe(gcf);
    img =  frame2im(frame);
    [img,cmap] = rgb2ind(img,256);
    if nb_frame == 1
        imwrite(img,cmap,'FlangedTube_animated.gif','gif','LoopCount',Inf,'DelayTime',1);
    elseif nb_frame == nb_frame_max/2 
        imwrite(img,cmap,'FlangedTube_animated.gif','gif','WriteMode','append','DelayTime',1);
    else
        imwrite(img,cmap,'FlangedTube_animated.gif','gif','WriteMode','append','DelayTime',0.05);
    end
    clf
end
