% Script Eady_ageostrophic.m
% Stability analysis of agestrophic Eady model 
%   using 2nd-order finite difference matrix method
% calls functions: baro, ddz, ddz2
% Used for problem A.20 of
%   Instability in Geophysical Flows
%   W.D. Smyth & J.R. Carpenter
%   Cambridge University Press, 2019
%
% Written using Matlab R2017a 
%   by W. Smyth, Oregon State University 2018
%   smythw@oregonstate.edu

%% preliminaries
clear
close all
fs=18;
ms=14;
lw=2;

%% set parameters, 
P=.01;
Rib1=100;
Rib2=.75;

%Rib1=.945;Rib2=Rib1;
Nk=45;
Nl=45;

% define vertical grid
N=60;
z_st=linspace(-1/2,1/2,N+1)';z_st(1)=[];z_st(end)=[];
dz=mean(diff(z_st))

fac=P; % plots will be scaled using P

%% loop over k 

k_sts=linspace(.1,3,Nk)*P;
for i=1:Nk
    disp(num2str(i/Nk));
    k_st=k_sts(i);
    [sig_k(i)]=baro(z_st,-1,Rib1,P*sqrt(Rib1),k_st,0);
end
% identify fastest-growing mode
[sig_k_mx,i_k_mx]=max(real(sig_k));
% refine fastest-growing mode using quadratic fit
ss=real(sig_k(i_k_mx+[-1 0 1]))
kk=k_sts(i_k_mx+[-1 0 1])
p=polyfit(kk-kk(2),ss,2);
kmx=kk(2)-.5*p(2)/p(1)
sig_kmx=baro(z_st,-1,Rib1,P*sqrt(Rib1),kmx,0)

%% loop over l 
l_min=1;
l_max=pi/dz;
l_sts=P*10.^linspace(log10(l_min),log10(l_max),Nl);
for i=1:Nl;
    disp(num2str(i/Nl));
    l_st=l_sts(i);
    [sig_l(i)]=baro(z_st,-1,Rib2,P*sqrt(Rib2),0,l_st);
end
% identify fastest-growing mode 
[sig_l_mx,i_l_mx]=max(real(sig_l));
% refine fastest-growing mode using quadratic fit
ss=real(sig_l(i_l_mx+[-1 0 1]))
ll=log10(l_sts(i_l_mx+[-1 0 1]))
p=polyfit(ll-ll(2),ss,2);
llmx=ll(2)-.5*p(2)/p(1);
lmx=10^llmx;
sig_lmx=baro(z_st,-1,Rib2,P*sqrt(Rib2),0,lmx)


%% plot results

figure
% l=0 (baroclinic modes)
subplot(2,1,1)
% plot numerical eigenvalues for l*=0
leg1(1)=plot(k_sts/fac,real(sig_k)/fac,'linewidth',lw) 
hold on
title(sprintf('(a) P=%.2f, N=%.0f, Ri_b=%.0f, FGM: \\sigma_r*/P=%.3f at k*/P=%.2f',P,N,Rib1,sig_kmx/P,kmx/P),'fontsize',fs,'Fontweight','normal')
ylabel('\sigma_r*/P','fontsize',fs)
set(gca,'fontsize',fs-2)
xlabel('k*/P','fontsize',fs)
set(gca,'fontsize',fs-2)

% add theoretical QG value
kmax_QG=1.61*P;
sigmax_QG=0.31*P;
plot(kmax_QG/fac,sigmax_QG/fac,'o','markersize',16,'linewidth',lw)

% k=0 (symmetric modes)
subplot(2,1,2)
% plot numerical eigenvalues 
leg1(2)=semilogx(l_sts/fac,real(sig_l)/fac,'linewidth',lw) % numerical eigenvalues
hold on
ylabel('\sigma_r*/P','fontsize',fs)

set(gca,'fontsize',fs-2)
xlabel('l*/P','fontsize',fs)
set(gca,'fontsize',fs-2)

smax_l=max(real(sig_l))/fac;
title(sprintf('(b) Ri_b=%.2f,     FGM: \\sigma_r*/P=%.3f at l*/P=%.2f',Rib2,sig_lmx/P,lmx/P),'fontsize',fs,'Fontweight','normal')

% add theoretical value
sUzP_theo=sqrt(1-Rib2);
if Rib2>1; sUzP_theo=0; end
semilogx(xlim,sUzP_theo*[1 1],'r','linewidth',1)
xlim([1 1000])
return