% a script to compute the spatial stability of the tanh mixing layer
% using Chebychev collocation method
% and augmented formulation for the polynomial eigenvalue problem
% the code is writen in primitive variables u,v,p

% the eigenvalue problem is (-om*E+A0+k*A1+k^2*A2)*X=0
%----->  k is given and om is the eigenvalue for the temporal problem (usual case)
%----->  om is given and k is the eigenvalue for the spatial problem (polynomial eigenvalue problem)

format compact; clear all

%%%% parameters
N=120 % number of collocation nodes
L=10 % domain from -L to L
Re=200 % Reynolds number based on velocity difference and vorticity thickness
om=4/3 % temporal pulsation (only for spatial stability)
k=1 %  spatial wavenumber (only for temporal stability)
spatemp=1 % type of analysis: 1 for spatial stability, 2 for temporal stability
boucon=1 % boundary conditions: 1 for neumann, 2 for dirichlet

%%%% differentiation matrices 
[y,DM] = chebdif(N,2); D1=DM(:,:,1);D2=DM(:,:,2); % chebychef
scal=L; y=y*scal; D1=D1/scal; D2=D2/scal^2; % scale domain to [-L,L]

%%%% base flow
u=0.5*tanh(2*y)+1.5; % base flow 
up=1./cosh(2*y).^2;% first derivative of base flow
figure(1);plot(u,y,'b',up,y,'r'); legend('u','up');xlabel('profile'); ylabel('y'); grid on; ylim([-L,L]);drawnow

%%%% coordinate selectors
uu=1:N; vv=uu+N; pp=vv+N; 
bou=[1,N]; % boundaries indices

%%%% dynamic matrices: 
Z=zeros(N,N);I=eye(N);
E=blkdiag(-i*I,-i*I,Z); % temporal operator
A0=[D2/Re, -diag(up), Z; Z, D2/Re, -D1; Z, D1, Z]; 
A1=[-i*diag(u), Z, -i*I; Z, -i*diag(u), Z; i*I, Z, Z]; 
A2=blkdiag(-I/Re,-I/Re,Z); 

switch spatemp
 case 1 % spatial stability analysis: k is the eigenvalue
  %%%% augmented formulation to transform polynomial eigenvalue problem into usual one
  %%%% could alternatively be solved using the "polyeig" function from matlab
  AA0=-om*blkdiag(E,Z,Z)+blkdiag(A0,I,I);
  AA1=[A1, [-I/Re, Z; Z, -I/Re; Z, Z]; [-I, Z, Z, Z, Z; Z, -I, Z, Z, Z]];
  bbou=[bou;bou+N; bou+3*N; bou+4*N];II=eye(5*N);DD=blkdiag(D1,D1,D1,D1,D1);
  xlimvec=[0, 2]; ylimvec=[-0.2,0.5];
  
 case 2 % temporal stability analysis: om is the eigenvalue
  AA0=A0+k*A1+k^2*A2;
  AA1=-E; 
  bbou=[bou;bou+N];II=eye(3*N);DD=blkdiag(D1,D1,D1);
  xlimvec=[0, 3]; ylimvec=[-2,1];
end

%%%% impose boundary conditions
AA1(bbou,:)=0;
switch boucon
 case 1; AA0(bbou,:)=DD(bbou,:); % Neumann boundary conditions
 case 2; AA0(bbou,:)=II(bbou,:); % dirichlet boundary conditions 
end

%%%% solve for eigenvalues
disp('computing eigenvalues')
tic; [U,S]=eig(AA0,-AA1); toc

%%%% remove small and large eigenmodes
S=diag(S); rem=abs(S)>100|abs(S)<1e-8; S(rem)=[]; U(:,rem)=[];
U=U(1:3*N,:); % remove ku and kv part of eigenvectors

%%%% plotting the eigenvalues/eigenvectors
figure(1);
for ind=1:length(S);
  %%%% plot one eigenmode
  h=plot(real(S(ind)),imag(S(ind)),'*'); hold on

  %%%%  plotting command for eigenmodes and callback function 
  tt=['figure(2);aa=' num2str(ind) ';plot(real(U(uu,aa)),y,''b'',imag(U(uu,aa)),y,''b--'',real(U(vv,aa)),y,''r'',imag(U(vv,aa)),y,''r--'',real(U(pp,aa)),y,''k'',imag(U(pp,aa)),y,''k--'');legend(''ureal'',''uimag'',''vreal'',''vimag'',''preal'',''pimag'');xlabel(''amplitude'');ylabel(''y'');'];
  set(h,'buttondownfcn',tt);
end
title('Click on eigenvalues to see the eigenmodes')
xlim(xlimvec); ylim(ylimvec); grid on; hold off