% % dispersion curves for tiger bush
clear all
 set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.2, ...
   'defaultlinelinewidth', 2.0, ...
   'defaultpatchlinewidth', 0.7); 

 

v = 20;
mt=2;
alist = [ 7,8];  % values of a to try
for k = 1:length(alist)
    at = alist(k);
%specify the steady state solutions    
wp = (at+sqrt(at^2-4*mt^2))/2;  
wm = (at-sqrt(at^2-4*mt^2))/2;  

np = mt/wp;
nm = mt/wm;

% calculate eigenvalues of the linearized system, as a function of omega
 
omlst = 3;  % largest value of omega to include
oml=[0:.02:omlst];
for kk = 1:length(oml);
    om = oml(kk);  
fn =  mt-om.^2;
fw = nm.^2;
gn = -2*mt;
gw=-1-nm.^2+i*v*om;
lam(:,kk)=eig([fn,fw;gn,gw]);  % eigenvalues of the Jacobian matrix
end
figure(1)
% plot the dispersion curve
plot(oml,real(lam(1,:)), oml,zeros(length(oml),1),'--','linewidth',2)
axis([0 2.5 -1.5 1.5])
xlabel('\omega','fontsize',20)
ylabel('Re(\lambda)','fontsize',20)
hold on

% now find the most unstable mode:
[mxlam,Jmx] = max(real(lam(1,:)))
al = mxlam
ommx = oml(Jmx)
 plot(ommx,al,'k*','linewidth',2);
 
bt = imag(lam(1,Jmx))
 wavelength = 2*pi/ommx  
 speed = -bt/ommx
 
end
 
text(1.1,0.3,'a=7','fontsize',20)
text(0.43,-1,'a=8','fontsize',20)
hold off
 
%the stability diagram
  am = [0:.01:10];
ma = [2:.01:11];
atr = ma.^2./(sqrt(ma-1));  %critical spiral point stability curve

% now compute the spatially inhomogeneous critical curve
gmlst = [0.01:0.025:5.5];
mlst = [0.01:0.025:5.];
[gm,m]=meshgrid(gmlst,mlst);
 
a=2*m+gm;
al=a./m;
n0=(al+sqrt(al.^2-4))/2;
n02=n0.^2;
b = 1+n02;
 
%the resultant has 24 terms; this was calculated using maple
R=b.^2.*m.^2.*v.^6 + 2.*b.^4.*m.*v.^4 + 4.*b.^3.*m.^2.*v.^4 ...
    - 8.*b.^3.*m.*n02.*v.^4 - 4.*b.^2.*m.^3.*v.^4 + 8.*b.^2.*m.^2.*n02.*v.^4 ...
    + 4.*b.*m.^3.*n02.*v.^4 + b.^6.*v.^2 + 8.*b.^5.*m.*v.^2 - 8.*b.^4.*m.^2.*v.^2 ...
    - 40.*b.^4.*m.*n02.*v.^2 + 76.*b.^3.*m.^2.*n02.*v.^2 - 40.*b.^2.*m.^3.*n02.*v.^2 ...
    - 32.*b.^2.*m.^2.*n02.^2.*v.^2 + 32.*b.*m.^3.*n02.^2.*v.^2 + 4.*m.^4.*n02.^2.*v.^2 ...
    + 4.*b.^7 - 4.*b.^6.*m - 24.*b.^5.*m.*n02 ...
    + 24.*b.^4.*m.^2.*n02 + 48.*b.^3.*m.^2.*n02.^2 - 48.*b.^2.*m.^3.*n02.^2 ...
    - 32.*b.*m.^3.*n02.^3 + 32.*m.^4.*n02.^3;
  
% plot the zero set of R
figure(2)
[C,H] = contour(gm,m,R,[0 0],'k','linewidth',2);
   % shift the zero set A values and then replot
hold on
 A=2*C(2,:)+C(1,:);
 
J = C(2,1);
   % the curve has a legitimate component and a spurious component
 % plot only the legitimate component
   
  figure(3)
Mx = max(A(2:J));
 plot(A(2:J),C(2,2:J),[0,2*Mx],[0,Mx],'linewidth',2)
  
 hold on
 % just to check, plot a ``*" at each point considered above, to see where
 % on the stability diagram they occur.
% for k = 1:length(alist)
%     plot(alist(k),mt,'k*')
% end
% % 
  
% add the spiral point critical curve and labels
plot(am,am/2,atr,ma,'--', 'linewidth',2)
axis([0 10 0 5])
text(3,2.8,'I','fontsize',20)
%text(2,3,'no vegetation','fontsize',18)
text(8,1.,'II','fontsize',20)
%text(4,1.5,'stable vegetation','fontsize',18)
text( 9,4.25,'III','fontsize',20)
%text(4,1.5,'unstable homogeneous vegetation','fontsize',18)
text(8,2.8,'IV','fontsize',20)
xlabel('a','fontsize',20)
ylabel('m','fontsize',20)
 

% 
% now find the speed and wavelength of the most unstable mode
 
 
for j = 1:length(gmlst)
     gam = gmlst(j);
    for k = 1:length(mlst)
     
        mt = mlst(k);
        at = 2*mt+gam;
% 
%find the steady state solutions    
wm = (at-sqrt(at^2-4*mt^2))/2;  
nm = mt/wm;

% calculate eigenvalues of the linearized system, as a function of omega
  
for kk = 1:length(oml);
    om = oml(kk);  
fn =  mt-om.^2;
fw = nm.^2;
gn = -2*mt;
gw=-1-nm.^2+i*v*om;
lam(:,kk)=eig([fn,fw;gn,gw]);  % eigenvalues of the Jacobian matrix
end
 
% % for this dispersion curve find the largest real(lambda)
% 
[mxlam1,Jmx1] = max(real(lam(1,:)));
 
[mxlam2,Jmx2] = max(real(lam(2,:)));
if(mxlam1>=mxlam2)
  
    lammx(j,k) = mxlam1; 
    btmx(j,k) = imag(lam(1,Jmx1));
     
    omndx(j,k) = Jmx1;
    ommx(j,k) = oml(Jmx1);
    wvl(j,k)=2*pi*(oml(Jmx1)>0)/((oml(Jmx1)<=0)+oml(Jmx1));
else
    
    lammx(j,k) = mxlam2;
    btmx(j,k) = imag(lam(2,Jmx2));
     omndx(j,k) = Jmx2;
    ommx(j,k) = oml(Jmx2);
    wvl(j,k)=2*pi*(oml(Jmx2)>0)/((oml(Jmx2)<=0)+oml(Jmx2));
end


    end
end

 %contour plot of wavelength
 figure(2)
 contour(gm,m,max(0,ommx'.*(lammx'>0)))
 xlabel('m','fontsize',20)
 ylabel('g','fontsize',20)
 hold off
 figure(3)
 a = 2*m+gm;
 contour(a,m,max(0,ommx'.*(lammx'>0)),'linewidth',2)
 colorbar
 hold off
  
%  