% stochastic reaction simulation, PHY825/BMS925

Na = 1000;   % molecules
Nb = 1000;
Nab = 0;

V = 1;      % uL

C1 = 1e-4;   % uL/(s*molecules*molecules)
Cn1 = 0.02;  % 1/(s*molecules)

clear molecdata
global molecdata

molecdata(1) = Na;
molecdata(2) = Nb;
molecdata(3) = Nab;

clear molectimedata
clear timedata

rxn{1}.subs = [1 1;
                1 2];
rxn{1}.prods = [1 3];
rxn{1}.rate = C1./V;
rxn{1}.time = inf;

rxn{2}.subs = [1 3];
rxn{2}.prods = [1 1;
                1 2];
rxn{2}.rate = Cn1;
rxn{2}.time = inf;

timeend = 1000;

tnext = 0;

i = 1;

timedata(i) = tnext;
molectimedata(i,:) = molecdata;

rxn = calc_rxn_times(rxn,tnext);
[tnext, rxn_num] = find_next_reaction(rxn);
    
while tnext < timeend
    
    tnext
    
    i = i+1;
    
    execute_reaction(rxn,rxn_num);
    
    molectimedata(i,:) = molecdata;
    timedata(i) = tnext;
    
    rxn = calc_rxn_times(rxn,tnext);
    
    [tnext, rxn_num] = find_next_reaction(rxn);
    
end

plot(timedata,molectimedata(:,1));