% Compute phase velocities and 1st mode eigenfunctions from Jon Nash's
% Columbia Plume data.
clear

imode = 1;

% Load the data
dat=load('Nash_data.txt');

% Extract the data
V=dat(:,1)*1;
rho=dat(:,2)*1000;
z=dat(:,3);

% Define constants
nu=1.e-6;
kap=nu/7;
g=9.81;
rho_0=mean(rho); 
% (You can choose the characteristic density value in different ways; it doesn't
%  make much difference to the results.)

% Convert density to buoyancy.
B=-g*(rho-rho_0)/rho_0;

% Plot the buoyancy and velocity profiles.
figure
subplot(1,2,1)
plot(B,z,'linewidth',2)
title('B [m/s^2]')
ylabel('z [m]')
subplot(1,2,2)
plot(V,z,'linewidth',2)
title('V [m/s]')

% Specify wavenumber. (For the hydrostatic limit, just choose a small, but
% nonzero, value).
k=.00001;

% Call SSF to compute the eigenvalue spectrum.

[sigs,w,b]=SSF(z,V,B,k,0,nu,kap,[0 0],[0 0],0);

% Convert growth rates to phase speeds.
cphs=-imag(sigs)/k;

% Sort modes by phase speed, fastest first.
[cphs,ind]=sort(cphs,'descend');
w=w(:,ind);b=b(:,ind);

% normalize so that max w = 1.
for i=1:length(ind)
    const=w(abs(w(:,i))==max(abs(w(:,i))),i);
    w(:,i)=w(:,i)/const;
    b(:,i)=b(:,i)/const;
end
u=(sqrt(-1)/k)*ddz(z)*w;

% Plot phase speeds
figure
plot(cphs,'b*')
xlabel('mode number')
ylabel('phase speed [m/s]')

% Plot eigenfunctions for fastest phase speed.
figure
subplot(1,3,1)
plot(real(w(:,imode)),z)
hold on
plot(imag(w(:,imode)),z,'r')
ylabel('z [m]')
xlabel('w eigfn')

subplot(1,3,2)
plot(real(u(:,imode)),z)
hold on
plot(imag(u(:,imode)),z,'r')
xlabel('u eigfn')
title('b=real, r=imag')

subplot(1,3,3)
plot(real(b(:,imode)),z)
hold on
plot(imag(b(:,imode)),z,'r')
xlabel('b eigfn')