% c is the nondimensional phase speed (nondimensionalized by Uo)
% z is the nondim height (nondimensionalized by H)

c = inline('1/2 + sgn*(1./muH).*sqrt((muH/2-coth(muH/2)).*(muH/2-tanh(muH/2)))','muH','sgn');

psihat=inline('cosh(muH*z)-real(c(muH,sgn))*sinh(muH*z)/(muH*abs(c(muH,sgn))^2)+ i*imag(c(muH,sgn))*sinh(muH*z)/(muH*abs(c(muH,sgn))^2)','z','c','muH','sgn');

muH=1.61;
z=[0:.05:1]';

clf
% subplot(4,1,1,'align')
subplot(3,1,1,'align')
ph=psihat(z,c,muH,1); % growing wave
plot(abs(ph),z,'k.-')
hold on
plot(angle(ph),z,'r--')
ph=psihat(z,c,muH,-1); % decaying wave
plot(angle(ph),z,'b-')
title(['\mu H = ',num2str(muH) ', Re(c/Uo) = ',num2str(real(c(muH,1))) ', Im(c/Uo) = ',num2str(imag(c(muH,1)))])
legend('|\psi(z)|','\alpha(z) of unstable wave','\alpha(z) of decaying wave')
ylabel('z/H')
grid

x=[0:.05:7]'; % nondim x (nondim by Ld)
kLd=muH;
sgn=1; 
[X,Z]=ndgrid(x,z);
t=0; % time in units of Ld/Uo
omega=kLd*c(muH,sgn); % frequency 
PHAT=psihat(Z,c,muH,sgn); % growing wave
PSI=real(PHAT.*exp(i*(kLd*X-omega*t)));
% subplot(4,1,2,'align')
subplot(3,1,2,'align')
contourf(x,z,PSI',10),colorbar
title('Perturbation \psi (pressure)')
xlabel('x/\lambda_d')
ylabel('z/H')

% Density perturbation (proportional to -\partial\psi'/\partial z)
RHO=-real(exp(i.*(kLd.*X - muH.*t.*(1./2 + sqrt((muH./2 - coth(muH./2)).*(muH./2 - ...
tanh(muH./2)))./muH))).*((i.*cosh(muH.*Z).*imag(sqrt((muH./2 - ...
coth(muH./2)).*(muH./2 - tanh(muH./2)))./muH))./abs(1./2 + sqrt((muH./2 - ...
coth(muH./2)).*(muH./2 - tanh(muH./2)))./muH).^2 - (cosh(muH.*Z).*(1./2 + ...
real(sqrt((muH./2 - coth(muH./2)).*(muH./2 - tanh(muH./2)))./muH)))./abs(1./2 ...
+ sqrt((muH./2 - coth(muH./2)).*(muH./2 - tanh(muH./2)))./muH).^2 + ...
muH.*sinh(muH.*Z)));

% subplot(4,1,3,'align')
subplot(3,1,3,'align')
contourf(x,z,RHO',10),colorbar
title('Perturbation \rho (density)')
xlabel('x/\lambda_d')
ylabel('z/H')

% subplot(4,1,4,'align')
% plot(x,PSI(:,1))
% hold on
% plot(x,-RHO(:,1),'r--')
% grid
% title('Surface perturbation pressure (\psi) and temperature (-\rho)')
% xlabel('x/\lambda_d')
% legend('\psi','-\rho')

bigfig('p')