%{
# Free-surface: displaying the eigenmodes using particles

In this code, we take the code that computes the eigenmodes for free surface waves [free_surface_gravity.m]() and we visualize the flow motion of the first eigenmode by expanding the solution to physical space and advecting tracer particles, like we did in a simple case in [particles.m]().
%}

clear all; clf;

n=100;      % number of gridpoints
alpha=1;    % wavenumber in x
L=2;        % Fluid height in y
rho=1;      % fluid density
mu=0.0001;    % fuid viscosity
g=1;        % gravity

% differentiation matrices
scale=-2/L;
[y,DM] = chebdif(n,2); 
D=DM(:,:,1)*scale;    
DD=DM(:,:,2)*scale^2;    
y=(y-1)/scale; 
I=eye(n); Z=zeros(n,n);

% renaming the matrices
dy=D; dyy=DD;
dx=i*alpha*I; dxx=-alpha^2*I;
Delta=dxx+dyy;

% System matrices
A=[mu*Delta, Z, -dx, Z(:,1); ...
   Z, mu*Delta, -dy, Z(:,1); ...
   dx, dy, Z, Z(:,1); ...
   Z(1,:),I(n,:),Z(1,:),0];

E=blkdiag(rho*I,rho*I,Z,1);

% boundary conditions
loc=[1,n,n+1,2*n];
C=[I(1,:),Z(1,:),Z(1,:),0; ... 
   Z(1,:),I(1,:),Z(1,:),0; ...
   Z(1,:),Z(1,:),-I(n,:),rho*g; ...
   dy(n,:),dx(n,:),Z(1,:),0]; 

E(loc,:)=0;  
A(loc,:)=C;

% compute eigenmodes
[U,S]=eig(A,E);
s=diag(S);  [t,o]=sort(-real(s)); s=s(o); U=U(:,o);
rem=abs(s)>1000; s(rem)=[]; U(:,rem)=[];
%{
# The animation

We show the position of the free surface stored in *q(eta)*, and the particles are initialy set as a mesh, and stretched in *y* to fit with the initial position of the free surface for the top prticles.
%}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% particle animation of the eigenmodes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% parameters
nper=1      % number of periods of oscillation
nt=30       % number of time steps per period
nx=40       % number of points in x
modesel=2;  % which more do animate

% select the eigenmode
u=1:n; v=u+n; p=v+n; eta=3*n+1;
q=U(:,modesel); 
lambda=s(modesel);

% The time and x extent
tvec=linspace(0,nper*2*pi/abs(lambda),nper*nt); 
dt=tvec(2)-tvec(1);
Lx=2*pi/alpha; 
x=linspace(0,Lx,nx);

% scale mode amplitude 
q=0.05*q/abs(q(eta)); 

% initialize tracer particles
[px,py]=meshgrid(linspace(0,Lx,60),linspace(0,L,30)); 
py=py.*(1+2*real(exp(lambda*tvec(1))*q(eta)*exp(i*alpha*px))/L);
px=px(:);py=py(:);  

% time loop
for ind=1:nper*nt
    
    
%{
# Expand the solution to physical space

We have assumed a wavelike behaviour of the flow solution in the $x$ direction, thus we have writen for all the flow variables, for instance for the $u$ velocity
$$
u(x,y,t)=\hat{u(y)}\exp(\alpha x+s t)+CC
$$
where $CC$ is the complex conjugate of the other term. This way we have transformed a variable $u$ depending on three coordinates, to a variable $\hat{u}$ depending on one single coordinate, thus transforming everything into a 1D problem. What we do now is to transform back the 1D problem into a full original problem by using the exponential. 

A complex quantity plus its complex conjugate gives twice the real part of it, this is what we code here.

In the code, we do the expansion to hysical space for all the variables $(u,v,p,\eta)$ at once in a single command.

%}
    % expand mode to physical space 
    qq=2*real(exp(lambda*tvec(ind))*q*exp(i*alpha*x));
        
    % plot pressure
    surf(x,y,qq(p,:)-10,'facealpha',0.3); view(2); shading interp; hold on
    
    % plot free surface
    plot(x,L+qq(eta,:),'k-',x,0*x+L,'k--'); hold on
    
    %%%% plot the particles
    plot(mod(px,Lx),py,'k.');

%{
# Interpolation of the particle positions

The particles are not on a cartesian grid, they can be wherever the field has brought them, the only thing that matter here for the interpolation is where they are on the $y$ grid, and then we can do the expansion in $x$ as a second step, this is less computations.
%}
    % compute velocity at position of particles
    pu=interp1(y,q(u),py);
    pv=interp1(y,q(v),py);

%{
# The Taylor expansion

When we wrote the boundary conditionat the free surface for setting up the linear model, we have used Taylor expansions to express the velocity and pressure at the position of the moving interface as a function of the values at the fixed position $L$. here we do again the same thing for the particles that are above $L$.

%}

    % For particles above L, use Taylor expansion for velocity
    pu(py>L)=q(u(n))+q(eta)*dy(n,:)*q(u);
    pv(py>L)=q(v(n))+q(eta)*dy(n,:)*q(v);
    
    % expand to physical space
    puu=2*real(exp(lambda*tvec(ind))*pu.*exp(i*alpha*px));
    pvv=2*real(exp(lambda*tvec(ind))*pv.*exp(i*alpha*px));
 
    % advect particles
    px=px+puu*dt;
    py=py+pvv*dt;
    
    xlabel('x');    
    ylabel('y'); 
    axis equal; axis([0,Lx,0,1.3*L]); 
    grid off
    hold off
    drawnow
 %{   
# To create the GIF animation
In this part, we add a gif animation to show the free surface
%}
    frame = getframe(1);
    im = frame2im(frame);
    [imind,cm] = rgb2ind(im,40);
    outfile = 'freegravitysurfaceparticles.gif';
        
            if ind==1
                imwrite(imind,cm,outfile,'gif','DelayTime',0,'loopcount',inf);
            else
                imwrite(imind,cm,outfile,'gif','DelayTime',0,'writemode','append');
        end

    end

    set(gcf,'paperpositionmode','auto');
    print('-dpng','-r100','free_surface_gravity_particles.png');