function temple9200_velocity_field
%TEMPLE9200_VELOCITY_FIELD
%   Constructs an example of an unsteady velocity field by
%   differentiating a stream function. The stream function
%   and the velocity field are plotted (the latter as a
%   quiver plot). In addition, markers are placed in a
%   part of the domain. These markers are then moved with
%   the flow, and their positions in time are shown.

% 08/2022 by Benjamin Seibold
%            http://www.math.temple.edu/~seibold/

%-----------------------------------------------------------------------
% Problem parameters
ax = [0 3 0 1]; % domain
nx = [ax(2)-ax(1),ax(4)-ax(3)]*20; % mesh points for plotting
tf = 5; % final time
nt = 200; % number of time steps
%-----------------------------------------------------------------------

% Plotting grid
x = linspace(ax(1),ax(2),nx(1)); % x-values
y = linspace(ax(3),ax(4),nx(2)); % y-values
[X,Y] = meshgrid(x,y); % 2d grid

[px,py] = meshgrid(linspace(.1,.5,80),... % insert markers
                   linspace(.4,.8,80)); % into domain

e = 1e-6; % small number (for differentiation)
dt = tf/nt; % time step size
for t = linspace(0,tf,nt+1) % time loop
    streamfun = @(x,y) y.^2/2-y.^3/3 ... % stream function: channel flow
        +cos(2*pi*x).*sin(2*pi*y).*(x-ax(1)).*(ax(2)-x)... % vortices
        *sin(pi*t)*.05; % time modulation
    dxsf = @(x,y) (streamfun(x+e,y)... % dx of
        -streamfun(x-e,y))/(2*e); % stream function
    dysf = @(x,y) (streamfun(x,y+e)... % dy of
        -streamfun(x,y-e))/(2*e); % stream function
    S = streamfun(X,Y); % evaluate stream function
    U = dysf(X,Y); V = -dxsf(X,Y); % evaluate velocity field

    % Plotting
    clf
    subplot(1,3,1)
    surf(x,y,S)
    axis([ax -.05 1/6+.05])
    xlabel('x'), ylabel('y')
    title(sprintf('Stream function at t=%0.1f',t))
    subplot(1,3,2:3)
    plot(ax([1 2 2 1 1]),ax([3 3 4 4 3]),'k-','LineWidth',1)
    hold on
    quiver(x,y,U,V)
    plot(px,py,'k-',px',py','k.-','LineWidth',.25,'MarkerSize',4)
    hold off
    axis equal, axis(ax+[-1 1 -1 1]*.02)
    xlabel('x'), ylabel('y')
    title(sprintf('Velocity field at t=%0.1f and Lagrangian markers',t))
    pause(1e-3)
    
    % Move markers with flow using forward Euler method
    px = px+dt*dysf(px,py);
    py = py-dt*dxsf(px,py);
end
