% Chebyshev spectral solution to Benard convection linear stability
% EVP

% We define a concatenated solution vector [u w T p] and set up
% EVP in blocks

  nu = 1/sqrt(1707.76/8);
  kappa = nu;
  kc = 3.117/2;

  % First, contour plot psi(x,z) and T(x,z) fields for fastest
  % growing mode at k=kc (which in this case has growthrate sigma=0).
  
  N = 8;
  getq = 1;
  [sigma, q] = benard(N,kc,nu,kappa,getq);
  real(sigma(1))
  [D,z] = cheb(N); % Get Chebyshev points in vertical
  z_intp = -1:0.1:1;
  x = (0:0.05:1); % x in units of wavelength lambda = 2*pi/kc
  psi = [0; q(1:N-1); 0];
  T = [0; q(N:2*N-2); 0];
  psi_intp=polyval(polyfit(z,psi,N),z_intp);
  T_intp = polyval(polyfit(z,T,N),z_intp);
  psi_field = real(exp(2*pi*1i*x')*psi_intp);
  T_field = real(exp(2*pi*1i*x')*T_intp);
  
  % Contour plot of psi
  
  subplot(2,2,1)
  contour(x,z_intp,psi_field',0:0.1:1)
  hold on
  contour(x,z_intp,psi_field',-0.1:-0.1:-1,'--')
  xlabel('x/\lambda')
  ylabel('z')
  title('\psi, CI = 0.1, negative contours dashed')
  hold off
  
  % Contour plot of T
  
  subplot(2,2,2)
  contour(x,z_intp,T_field',0:0.2:1)
  hold on
  contour(x,z_intp,T_field',-0.2:-0.2:-1,'--')
  xlabel('x/\lambda')
  ylabel('z')
  title('T, CI = 0.2, negative contours dashed')
  hold off

  % Plot convergence of numerically computed growth rate
  % vs. N. Assume (20/20 hindsight) that N=16 solution is
  % essentially exact and plot N = 2,...,16 sigma differences
  % from this solution. The N = 12-16 points all give zero values
  % and don't show up on the semilogy plot.
  
  subplot(2,2,3)
  getq = 0;
  mmax = 8;
  growthrate = zeros(1,mmax);
  dN = 2;
  for m = 1:mmax
    [sigma, junk] = benard(m*dN,kc,nu,kappa,getq);
    growthrate(m) = real(sigma(1));
  end
  semilogy((1:mmax)*dN,abs(growthrate-growthrate(mmax)),'x')
  xlabel('N')
  ylabel('|\sigma_f - \sigma_f(N=16)|, k = k_c')
  title('Convergence vs. Chebyshev order N')

  % Lastly, plot growthrate vs. wavenumber k for 0.5 < k < 3,
  % for default, 0.5x, and 2x nu and kappa.
  subplot(2,2,4)
  N = 8;
  krange = 0.5:0.1:3;
  nk = length(krange);
  gr = zeros(3,nk);
  dampfact = [1 0.5 2];
  for j = 1:nk
    for ifact =1:3
      [sigma, junk] ...
        = benard(N,krange(j),nu*dampfact(ifact),kappa*dampfact(ifact),getq);
      gr(ifact,j) = real(sigma(1));
    end % ifact loop
  end % j loop
  plot(krange,gr)
  xlabel('k')
  ylabel('\sigma_f')
  legend('\nu_n = 0.0684','\nu = 0.5\nu_n','\nu = 2\nu_n',3)
  title('Growth rate (\nu = \kappa) vs. wavenumber (N=8)')
  hold on
 % line([kc kc],[0 -0.01]) % Draw a tick mark at k = kc
  plot(kc,0,'bx')
  text(kc,-0.2,'k_c')
  hold off