
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Calculate the eigenfrequencies of the LTE
%  using the formulation of Kasahara, MWR, 1976.
%
%  This version plots dispersion diagrams.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Input: 
%
% Number of modes in each category
% Lamb parameter.
% 

clear
clf %  Clear the plot window.

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

s = 1    %   Zonal wavenumber 1.

nmodes = 50 % Reduce to 5 for quick run;
 
neps = 40;
epsmax=10000; epsmin=0.01;
repsmax=sqrt(epsmax);repsmin=sqrt(epsmin);
lrepsmax=log10(repsmax);lrepsmin=log10(repsmin);
logrooteps = linspace(lrepsmin,lrepsmax,neps);
lreps = logrooteps;
rooteps = 10.^(lreps);
eps = rooteps.^2;

for ne=1:neps
  epsr = eps(ne);

%   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);
    
% 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.
ne=1:neps;
xx=lreps;
lx=rooteps;
nplotmodes = min(10,nmodes);

for n=1:2*nplotmodes
  ww=wg(ne,n)';
  rr=rh(ne,n)';
  ee=eg(ne,n)';

  hold on, subplot(2,2,1), loglog(lx,ww'), 
           title('\bf Westward Gravity'), axis([0.1,100,-10,-0.01])
           xlabel('\surd\epsilon'),
           ylabel('Frequency \sigma')
  hold on, subplot(2,2,2), loglog(lx,ee'), 
           title('\bf Eastward Gravity'), axis([0.1,100,0.01,+10])
           xlabel('\surd\epsilon'),
           ylabel('Frequency \sigma'), text(5,0.05,'KW')
  hold on, subplot(2,2,3), loglog(lx,rr'), 
           title('\bf Rotational'), axis([0.1,100,-1.0,-0.001])
           xlabel('\surd\epsilon'),
           ylabel('Frequency \sigma'), text(5,-0.5,'RG') 

  hold on, subplot(2,2,4), 
    arr = asinh(rr); aee = asinh(ee); aww = asinh(ww);
%   semilogx(lx,arr','y-',lx,aee','y--',lx,aww','y--'), 
    semilogx(lx,arr','b-',lx,aee','b--',lx,aww','b--'), 
    title('\bf All Modes'), axis([0.1,100,-2.0,+2.0])
%   xlabel('root(eps)'), ylabel('arcsinh(frequency)')  
    xlabel('\surd\epsilon'), ylabel('sinh^{-1}\sigma')  
    text(.3,-0.7,'RG'), text(.3,1.1,'KW')

end

%%%   Save the image
%%%   print -deps 'eigen1.eps'
      pause; figure  %%%   Pause and then open new plot window.

for n=1:2*nplotmodes
  ww=wg(ne,n)';
  rr=rh(ne,n)';
  ee=eg(ne,n)';

  hold on, subplot(2,2,1), loglog(lx,ww'), 
           title('\bf Westward Gravity'), axis([0.1,100,-10,-0.01])
           xlabel('\surd\epsilon'),
           ylabel('Frequency \sigma')
  hold on, subplot(2,2,2), loglog(lx,ee'), 
           title('\bf Eastward Gravity'), axis([0.1,100,0.01,+10])
           xlabel('\surd\epsilon'),
           ylabel('Frequency \sigma'), text(5,0.05,'KW')
  hold on, subplot(2,2,3), loglog(lx,rr'), 
           title('\bf Rotational'), axis([0.1,100,-1.0,-0.001])
           xlabel('\surd\epsilon'),
           ylabel('Frequency \sigma'), text(5,-0.5,'RG') 

  hold on; subplot(2,2,4) 
    ah = atanh(0.5); sig0 = 0.25;
    arr = tanh(ah*rr/sig0); aee = tanh(ah*ee/sig0); aww = tanh(ah*ww/sig0);
%   semilogx(lx,arr','y-',lx,aee','y--',lx,aww','y--'), 
    semilogx(lx,arr','b-',lx,aee','b--',lx,aww','b--'), 
    title('\bf All Modes'), axis([0.1,100,-1.0,+1.0])
%   xlabel('root(eps)'), ylabel('tanh(frequency)')
    xlabel('\surd\epsilon'), ylabel('tanh(\sigma/\sigma_0)')
    text(1.0,-0.65,'RG'), text(1.0,0.6,'KW')

end

%%%   Save the image
      print -deps 'LTEeig.eps'

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
