clear;

% the following simmulates the following reactions
%  X1 --> X1 - X2
%  X2 --> X2 - X1


% in other words, dX2/dt = - s1* X1 and dX1/dt = - s2 * X2

%set the rate constants
speed= [1; 
        0.5];


%specify the reactants  (left hand sides of the reactions)
Reactants = [ 1, 0;              %rows are reactions, columns are reactants
              0, 1];

Products = [ 1,-1;         %(right hand sides of the reactions)
            -1, 1];            %rows are reactions, columns are products
        
%Specify the change matrix (+outcomes, - reactants), or right hand side - left hand side of the reaction
Change = Products-Reactants;        %rows are reactions, columns are outcomes/reactants
          

%initialize the state space
x = [ 50, 70];   %first is how many are of the firct chemical,...



%technical variables, initialization
Number_Of_Reactions = size(speed,1);  %get the number of rows in the speed matrix
%Number_Of_Reactants = size(Reactants,2);     %get the number of collums in the reactant matrix

step = 1;  
X(step,:) = x; %will track the reactants over time (steps)
T(step) = 0;  %will track the time, start at 0

Max_Number_Of_Steps = 10000; % how many steps we want to go, could be replaced by actual time if needed


ONES=ones(Number_Of_Reactions,1); %this is a vector of a lot of ones, ONES*x will create a matrix with a lot of rows, each row consisting of a vector x

%algorithm
while (min(x)>0 & step< Max_Number_Of_Steps) %while the reactions make sense and we did not simulate too much
    %note that the above reactions are exceptional and mix(x) should in
    %most cases be replaced by max(x)
    step = step+1;
    h = speed.*prod((ONES*x).^Reactants,2);  %multiply point by point, i.e. h(2) = speed(2)*(x(1)^(occurence in 2nd reaction)*x(2)^(occurence in 2nd reaction)... 
       
    hc = cumsum(h);  %hc(1) = h1, "urgency" of reaction 1, hc(2) = hc(1)+hc(2), "urgency" of any of reactions 1 and 2
    H = hc(Number_Of_Reactions);  % "urgency of any reaction to occur"

    %based on urgency, calculate when the reaction will occur
    rn = rand(1);
    T(step) = -log(rn)/H +T(step-1);

    %based on urgency, calculate what reaction it will be
    rn = rand(1);
    rk = min(find(rn<=hc/H)); % this determines which reaction occurs, reaction 2 occurs if rn>h(1) and rn<h(1)+h(2)
    
    %perform the reaction
    x = x + Change(rk,:);
    X(step,:) = x; %tracks the reactants

end

figure(1)
plot(T,X) %plot how reactents depend on time

figure(2)
plot(X(:,1),X(:,2))  %plot how reactants "depend" on each other

figure(3)
plot(T, X(:,1)./(X(:,1)+X(:,2)))  %plot proportion of first reactant