% code to simulate tiger bush via Method of lines, using data 
function bush_sim
 global a m v N dx 
 %dimensional parameters from Table 12.1
 A = 300;
 E = 4;
 V = 200;
 D = 1;
 
 %set the case
 veggies = 2; %grass
 
 veggies = 1 % trees
 
 switch veggies
     case 1
 % for trees
 B = 0.002;
 M = 0.18;
 R = 1.5;
 
     case 2
 % for grass
 B = 0.003;
 M = 1.8;
 R = 100;
 
 end
 %A=250;
 E_crit = R*(A*B/(2*M) )^2 
 
 %non-dimensional parameters
 eta = sqrt(E/R);
 rho = E/(eta*R*B);
 a = A/(rho*E)  %water source rate
 m = M/E % plant death (or harvest) rate
 v = V/sqrt(D) % advective velocity
 
 
% for these parameters estimate the speed and wavelength from linear theory
at = a;
mt = m;
%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:.002:omlst];
for j = 1:length(oml);
    om = oml(j);  
fn =  mt-om.^2;
fw = nm.^2;
gn = -2*mt;
gw=-1-nm.^2+i*v*om;
lam(:,j)=eig([fn,fw;gn,gw]);  % eigen values of the Jacobian matrix
end
figure(1)
% plot the dispersion curve
plot(oml,real(lam(1,:)),oml,zeros(length(oml),1),'--',oml,real(lam(2,:)),'r','linewidth',2)
 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,'*','linewidth',2);
 
 bt =  imag(lam(1,Jmx));
 ndwavelength = 2*pi/ommx  
 ndspeed = bt/ommx
 hold off
 

 L = ndwavelength; % length of domain; pick this to be  the wavelength of the fastest growing mode
 N = 201; %number of cells 
 dx = L/(N-1);
 x = [0:N-1]*dx;
 
 %nullclines
 w = [0.001:.01:a];
n1 = m./w;
w2 = [0.001:0.01:a];
n2 = sqrt(a./w2-1);

wp = (a+sqrt(a^2-4*m^2))/2;
wm = (a-sqrt(a^2-4*m^2))/2;

np = m/wp;
nm = m/wm;

% setup for pde solve
 % specify initial data
 w0 = a*ones(N,1);
 u0 = m./w0 + sin(2*pi*x'/L) ;
 [umx,jmx] = max(u0);
 X(1) = x(jmx);
 s0 = [w0;u0];
 switch veggies
     case 2
tstep =   0.5;
t_end = 40;
     case 1
         tstep =   5;
t_end = 400;
 end

         
%specify the output points
tspan = [0:tstep:t_end];
 
[T,S] = ode15s(@deRHS,tspan, s0);  

 %plot the solution as a movie
figure(2)
 
for j = 1:length(T)
    plot(x,S(j,N+1:2*N),x,S(j, 1:N),'b' )
    
    % here we assume there is a unique location of the maximum
    %i.e., a single wavelength
    
    [mxu,mxj]=max(S(j,N+1:2*N));  % find the position of maximum
   % X tracks the location of the max
     
     X(j) = x(mxj);
     
    pause(.1)
end
 
% calculate the wavespeed
shift = 0;
Xt(1) = X(1);
for j = 2:length(T)
  if(X(j)<X(j-1))
      shift = shift+1;
  end
  Xt(j) = X(j) + L*shift;
end

% show the trajectory as a function of time
figure(4)
plot(T,Xt,'linewidth',2)
xlabel('t','fontsize',20)
xmax = Xt(end);
xmin = Xt(4);
 vel = (xmax-xmin)/(T(end)-T(4));
 
 % plot the final pattern
 
 % convert back to dimensional units
speed = vel*sqrt(E*D) %units of meters per year 
xscale = sqrt(D/E);
wavelength = ndwavelength*xscale % units of meters

 figure(3)
 formatSpecF = '%6.2f\n';
plot(x*xscale,S(length(T),N+1:2*N),'b',x*xscale,S(length(T), 1:N),'r','linewidth',2 )
 title(strcat('Velocity = ',sprintf(formatSpecF,vel*sqrt(D*E)),'meters/year'),'fontsize',18)
  xlabel('x (meters)','fontsize',20)
legend('plant density','water density','fontsize',20)
axis([0 wavelength 0 fix(max(S(length(T),N+1:2*N))+1)])
 set(gca,'Ytick',0:2.5:15)

%the right hand side for ode simulation:
function s_prime=deRHS(t,y)
global a m v N dx
 w = y(1:N); %water content
 u = y(N+1:2*N); %number density of plants
 
 fw = a -w-w.*u.^2;
 
 fu = w.*u.^2-m*u;
%using periodic boundary conditions 
Fw = fw + v*([w(2:N);w(1)]-w)/dx;% forward difference is upwinfing assuming v>0
 
 Fu = fu+ (-2*u+[u(end);u(1:end-1)]+[u(2:end);u(1)])/dx^2;
 
s_prime = [Fw;Fu];

