%% Foreword
 % Written by: A.M. Wouden, Penn State University
 % Tailored from MathCad code of: J.M. Cimbala, Penn State University
 
 % Theory can be found in:
 %            Kundu, Cohen, Dowling, "Fluid Mechanics", 6ed., p544-53

%% Preliminary Setup - Convienience arguments

clc
clear
close all
format compact
format shortG
options = optimset('Display','off');

cbrt = @(x) x^(1/3);
sq = @(x) x^2;

%% Wave Number Setup
ntotal = 500;
Kmin = 1;
Kmax = 8;
K=Kmin:(Kmax-Kmin)/(ntotal-1):Kmax;

%% Since this appears frequently
RaK = @(K,Ra) cbrt(Ra/(K^4));

%% Define the roots q0, q, and qs as functions of the wavenumber K and
 % Rayleigh number Ra:
q0  = @(K,Ra) K*sqrt(RaK(K,Ra)-1);
   
q   = @(K,Ra) sqrt(K^2*(1+0.5*RaK(K,Ra)*(1+sqrt(3)*1i)));
    
qs  = @(K,Ra) sqrt(K^2*(1+0.5*RaK(K,Ra)*(1-sqrt(3)*1i)));

%% Now define the matrix as given in the class notes
 % (Matrix M is also a function of K and Ra):
M   = @(K,Ra) ...
    [ 
        cos(q0(K,Ra)/2) ...
       cosh( q(K,Ra)/2) ...
       cosh(qs(K,Ra)/2);
      -q0(K,Ra)* sin(q0(K,Ra)/2) ...
        q(K,Ra)*sinh( q(K,Ra)/2) ...
       qs(K,Ra)*sinh(qs(K,Ra)/2);
       sq(q0(K,Ra)^2+K^2)* cos(q0(K,Ra)/2) ...
       sq( q(K,Ra)^2-K^2)*cosh( q(K,Ra)/2) ...
       sq(qs(K,Ra)^2-K^2)*cosh(qs(K,Ra)/2);
    ];
 % Define the determinant of this matrix
dm = @(K,Ra) det(M(K,Ra));

%% Use fsolve(func,guess) to solve for the corresponding Ra which yields
 % non-zero roots
RaGuess = 4000;
RaSol = zeros(ntotal,1);
for iter=1:ntotal
    RaSol(iter) = fsolve(@(Ra) dm(K(iter),Ra), RaGuess, options);
end

%% Write the minimum Ra for stability, along with its wavenuber, K
[RaMin,id] = min(RaSol);

display('Minimum Ra for Stability');
disp('K = ')
disp(K(id))
disp('Ra = ')
disp(real(RaSol(id)))
%% Plot the marginal (neutral) stability curve
hold on
plot(real(RaSol),K,...
     '-k',...
     'LineWidth', 3, 'Color', 'red')
grid on
xlim([0 8000])
ylim([K(1)-1 K(ntotal)+1])
xlabel('Rayleigh number, \it Ra \rm')
ylabel('Wavenumber, \it K \rm')
title({'-- Benard Convection Problem --';...
       'Eigenvalue Solution for Marginal Stability'})
   
str1 = {'STABLE';'\sigma < 0'};
text(2000,8,str1,'horizontalAlignment','center')
str2 = {'UNSTABLE';'\sigma > 0'};
text(5000,4,str2,'horizontalAlignment','center')
str3 = {'marginallly stable';'(\sigma = 0)'};
text(real(RaSol(1))+1000,K(1),str3,'horizontalAlignment','center')

rectangle('Position',[0 K(1)-1 real(RaSol(id)) K(id)])
str4 = strcat('K = ',num2str(K(id)));
text(800,K(id+20),str4,'horizontalAlignment','center')
str5 = strcat('  Ra = ',num2str(real(RaSol(id))));
text(real(RaSol(id)),0.5,str5)