% script SSF_plunging
% Stability analysis of viscous stratified shear flow
% Calls functions: SSF, ddz, ddz2
% Used for problem A.18 of
%   Instability in Geophysical Flows
%   W.D. Smyth & J.R. Carpenter
%   Cambridge University Press, 2019
%
% Written using Matlab R2017a 
%   W. Smyth, Oregon State University, 2018
%   smythw@oregonstate.edu

clear

% set parameters
Lz_st=8;
Ri_bs=[0 0.1 0.1 0.1];
Pr=1;
Res=[1e6 1e6 100 20];

% define z values
del=0.2;
z_st=[-Lz_st/2+del:del:Lz_st/2-del]'; % exclude boundaries
nz=length(z_st)
labs=['(a) ';'(b) ';'(c) ';'(d) '];

figure

% loop over Ri, Re
for j=1:4
    Ri_b=Ri_bs(j);
    Re=Res(j);
    
    % define profiles
    U_st=-tanh(z_st);
    
    B_st=Ri_b*tanh(z_st);
    Bz_st=Ri_b*sech(z_st).^2;
    
    % loop over k
    ks=[0.0:.01:1];nks=length(ks)
    l_st=0;
    for i=1:nks
        k_st=ks(i)
        [s_st,w_st,beta_st]=SSF(z_st,U_st,Bz_st,k_st,l_st,1/Re,1/(Re*Pr));
        sig_st(i)=real(s_st(1));
    end
    
    % plot results
    subplot(2,2,j)
    hold on
    plot(ks,sig_st,'k','linewidth',1.5)
    xlabel('k*')
    ylabel('\sigma*')
    title('LSA of stratified shear flow')
    axis tight
    yl=ylim;
    xlim([0 1])
    ylim([0 yl(2)])
    
    srmax=max(sig_st);
    kmax=ks(sig_st==srmax);
    ttle=sprintf('Ri_b=%.2f, Re=%.0e, Pr=%.0f, \\sigma_r*=%.4f, k*=%.2f',Ri_b,Re,Pr,srmax,kmax)
    title([labs(j,:) ttle],'fontweight','normal')
end
