
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Calculate the eigenfrequencies of the LTE
%  using the formulation of Kasahara, MWR, 1976.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input: 
%
% Number of modes in each category
% Lamb parameter.
%  (Sometimes: Zonal wavenumber)
clf %  Clear the plot window.

wg = [];
rh = [];
eg = [];

smax=10
nmodes = 6  
epsr = 8.7484
for s=1:smax
%   Set up the A-matrix
    nn = 3*nmodes;
    lm=nn/3-1;
    
    A=zeros(nn,nn) ;
    
    for i=1:lm
      m=3*i-1;
      t1=s+2*(i-1);
      t2=t1+1;
      t3=t1+2;
      t4=t1+3;
      t5=t1+4;
      A(m,m-1)=t1*(t2-s)/(t2*(2*t1+1));
      A(m,m)=-s/(t2*t3);
      A(m,m+2)=t4*(t3+s)/(t3*(2*t3+1));
      A(m+1,m-1)=-t1*t2/epsr;
      A(m+2,m  )=t2*(t3-s)/(t3*(2*t2+1));
      A(m+2,m+2)=-s/(t3*t4);
      A(m+2,m+3)=t5*(t4+s)/(t4*(2*t4+1));
      A(m+2,m+4)=-1;
    end
    A(1,1)=-1/(s+1);
    A(1,2)=(s+2)*(2*s+1)/((s+1)*(2*s+3));
    A(1,3)=-1;
    t1=s+2*(nn/3-1);
    t2=t1+1;
    t3=t1+2;
    A(nn-1,nn-2)=t1*(t2-s)/(t2*(2*t1+1));
    A(nn,nn-2)=-t1*t2/epsr;
    A(nn-1,nn-1)=-s/(t2*t3);
    
%   Set up the B-matrix
    
    B=zeros(nn,nn) ;
    
    for i=1:lm
      m=3*i-1;
      t1=s+2*(i-1);
      t2=t1+1;
      t3=t1+2;
      t4=t1+3;
      t5=t1+4;
      B(m,m)=-s/(t2*t3);
      B(m,m-1)=t1*(t2-s)/(t2*(2*t1+1));
      B(m,m+1)=-1;
      B(m,m+2)=t4*(t3+s)/(t3*(2*t3+1));
      B(m+1,m)=-t2*t3/epsr;
      B(m+2,m  )=t2*(t3-s)/(t3*(2*t2+1));
      B(m+2,m+2)=-s/(t3*t4);
      B(m+2,m+3)=t5*(t4+s)/(t4*(2*t4+1));
    end
    B(1,1)=-1/(s+1);
    B(1,2)=(s+2)*(2*s+1)/((s+1)*(2*s+3));
    t1=s+2*(nn/3-1);
    t2=t1+1;
    t3=t1+2;
    B(nn,nn-1)=-t2*t3/epsr;
    B(nn-1,nn  )=-1;
    B(nn-1,nn-1)=-s/(t2*t3);
    B(nn-1,nn-2)=t1*(t2-s)/(t2*(2*t1+1));
    
    evalA=eig(A);
    evalB=eig(B);
    
%    subplot(2,2,1), spy(A), title('A')
%    subplot(2,2,2), spy(B), title('B')
    
% Put symmetric and antisymmetric modes together. 
    evalAB = [ evalA ; evalB ];
    evals = sort(evalAB)';
    nas=2*nmodes;
    WG = evals(1:nas);
    RH = evals(nas+1:2*nas);
    EG = evals(2*nas+1:3*nas);
% split into the three types and sort out
    WG = -sort(-WG);
    wg = [wg ; WG ];
    rh = [rh ; RH ];
    eg = [eg ; EG ];
end
%    Plot the eigenvalues against zonal wavenumber.
s=1:smax;
for n=1:2*nmodes
  ww=wg(s,n)';
  rr=rh(s,n)';
  ee=eg(s,n)';
  hold on, subplot(2,2,1), plot(s',ee'), title('eg')
  hold on, subplot(2,2,2), plot(s',ww'), title('wg')
  hold on, subplot(2,2,3), plot(s',rr'), title('rh')
end
% Plot the five day wave.
n=2
rr=rh(s,n)';
hold on, subplot(2,2,4), 
plot(s',rr','y-'), hold on, plot(s',rr','yo'),
axis([0,6,-0.18,-0.09])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Now overplot the approximate values given by the
%   simplified formulae for RH and pure gravity waves.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% nmodes = 6; epsr = 8.7484; these are set above.
mmax=smax;
nmax=2*nmodes;
for m=1:mmax
for n=1:nmax
  n0=n;
  n=n+m-1;
  nurh(m,n0) = -m/(n*(n+1));
  nuwg(m,n0) = -sqrt(n*(n+1)/epsr);
  nueg(m,n0) = +sqrt(n*(n+1)/epsr);
end
end
% First symmetric RH wave.
nufdw = nurh(:,2);
for n0=1:nmax
  for m=1:mmax
    x(m)=m;
    swg(m)=nuwg(m,n0);
    srh(m)=nurh(m,n0);
    seg(m)=nueg(m,n0);
  end
  hold on, subplot(2,2,1)  %, axis square
           plot(x,seg,'--g'), title('\bf Eastward Gravity')
           ylabel('Frequency \sigma') 
  hold on, subplot(2,2,2)  %, axis square
           plot(x,swg,'--g'), title('\bf Westward Gravity')
           ylabel('Frequency \sigma')
  hold on, subplot(2,2,3)  %, axis square
           plot(x,srh,'--g'), title('\bf Rotational Modes')
           xlabel('Zonal wavenumber m');
           ylabel('Frequency \sigma')
end
% Five day wave.
hold on
xax = [ 1 2 3 4 5 6];
PL4 =  subplot(2,2,4)  %, axis square
plot(x,nufdw,'y--'), plot(x,nufdw,'yx'), 
set(PL4,'XTick',xax);
text(1.2,-0.1,'5DW'); text(0.3,-0.158,'3DW');
title('\bf First Symmetric Mode')
xlabel('Zonal wavenumber m');
ylabel('Frequency \sigma')
%%%   Save the image
print -deps 'LTEkas.eps'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
