%
%  NRM_LV_timetrack.m	Run Gillespie's algorithm on a Lotka Volterra model.
%      
%  usage:	GILL_LV_timetrack(tend)
%
%               tend = max time for the program to run = 20
%
%  example:	GILL_LV_timetrack(20)
%

function [X] = GILL_LV_timetrack(tend)

%maximum number of steps to run before quitting.
maxi = tend*10^5;

%rate constants.
k = [10, 0.01, 0.01, 0.01, 10];

%reaction vectors (columns);
zeta= [1  0  -1  -1  0;
       0  1   0   0 -1];      


%initialize.  Assume start with 300 animals each.
initial = [300,300]';

%T is time.
T = zeros([maxi,1]);

%state of the system.
x = zeros([2,maxi]);

%initialize the state of the system
x(:,1) = initial;

%holder for intensities.
lambda = zeros([5,1]);


%main loop
for i =1:maxi-1
    
    %set intensity functions.
    lambda(1) = k(1)*x(1,i);
    lambda(2) = k(2)*x(1,i)*x(2,i);
    lambda(3) = k(3)*x(1,i)*x(2,i);
    lambda(4) = k(4)*x(1,i);
    lambda(5) = k(5)*x(2,i);
    
    %find the sum of the intensity functions
    lambda0 = lambda(1) + lambda(2) + lambda(3) + lambda(4) + lambda(5);
    
    %generate a RV for finding when the next reaction takes place.
    r1 = rand;
    
    %find the time of the occurrence.
    t = log(1/r1)/lambda0;
    
    %if the time takes us past the end time, we should quit.
    if T(i) + t > tend
        T(i+1) = tend;
        x(:,i+1) = x(:,i);
        break
    end
    
    %find which reaction occured.
    
    r2 = rand;
    
    if r2 < lambda(1)/lambda0 
        loc = 1;
    elseif r2 < (lambda(1) + lambda(2))/lambda0
        loc = 2;
    elseif r2 < (lambda(1) + lambda(2) + lambda(3))/lambda0
        loc = 3;
    elseif r2 < (lambda(1) + lambda(2) + lambda(3) + lambda(4))/lambda0
        loc = 4;
    else
        loc = 5;
    end
    
    %update the state of the system and the time.
    x(:,i+1) = x(:,i) + zeta(:,loc);
    T(i+1) = T(i) + t;
    
end

count = i;
X = x(count+1);
%plot the timecourse of the state of the system.

figure
plot(T(1:count),x(1,(1:count)),T(1:count),x(2,(1:count)))
     legend('Rabbits','Foxes')

     %     title('Full version')

%Now for the exact solution
% 
% Num = 1000000;
% Delta = tend/Num;
% 
% ex = zeros([2,Num+1]);
% ex(:,1) = initial;
% T = zeros([Num+1,1]);
% for i = 1:Num
%     lambda(1) = k(1)*ex(1,i);
%     lambda(2) = k(2)*ex(1,i)*ex(2,i);
%     lambda(3) = k(3)*ex(1,i)*ex(2,i);
%     lambda(4) = k(4)*ex(1,i);
%     lambda(5) = k(5)*ex(2,i);
%     
%     ex(:,i+1) = ex(:,i) + zeta*lambda*Delta;
%     T(i+1) = T(i) + Delta;
% end
% 
% figure
% plot(T,ex(1,:),T,ex(2,:))
% legend('Rabbits','Foxes')
end %The program