% 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')