% Program to simulate a simple model of a unimolecular reaction with 5
% equally spaced energy levels, and reaction from the last level.

% Make sure that we don't start with any predefined values.
%clear

% Transition rates in s^{-1}
wdown = 1e9
wup = 0.8*wdown
wr = 1e13

% Number of molecules and initial populations:
t = 0;
Ntot = 10000
N = zeros(10,1);
N(1) = Ntot;

% A counter for the number of simulation steps
nsteps = 0;

% Main simulation loop: continue until we run out of molecules.
while (Ntot > 0)
    % Calculate propensities.
    a(1:4) = wup*N(1:4);    % You can only go up from levels 1 to 4.
    a(5:8) = wdown*N(2:5);  % You can only go down from levels 2 to 5.
    a(9) = wr*N(5);         % Reaction from level 5.
    
    % KMC selection of time and next jump.
    a_sums = cumsum(a);
    atot = a_sums(end);
    r = rand(2,1);
    Delta_t = -log(r(1))/atot;
    k = find(a_sums > r(2)*atot,1,'first');
    
    % What reaction did we pick? Change the populations accordingly.
    if (k<5)
        % This is an upward jump.
        start_state = k;
        N(start_state) = N(start_state) - 1;
        N(start_state+1) = N(start_state+1) + 1;
    elseif (k<9)
        % Downward jump.
        start_state = k - 3;
        N(start_state) = N(start_state) - 1;
        N(start_state-1) = N(start_state-1) + 1;
    else    % k = 9: reaction
        N(5) = N(5) - 1;
        Ntot = Ntot - 1;
    end
    
    % Increment t, then store the current total number of molecules.
    t = t + Delta_t;
    nsteps = nsteps + 1;
    tval(nsteps) = t;
    Ntotval(nsteps) = Ntot;
end

% Simulation done. Plot ln(Ntot) vs t.
% The last value of Ntotval will always be zero so we can't take its log.
% Get rid of this point.
nsteps = nsteps - 1;
tval = tval(1:nsteps);
Ntotval = Ntotval(1:nsteps);
plot(tval,log(Ntotval))
xlabel('t/s')
ylabel('ln(N_{total})','interpreter','tex')
hold    % This will allow us to add a line to the plot.

% Fit the data to get a rate constant.
% line_coeffs(1) will be the slope and line_coeffs(2) the intercept.
line_coeffs = polyfit(tval,log(Ntotval),1);
% The first time point is at t=0. The last time point is stored in t.
t_endpoints = [0,t];
% Calculate the corresponding N values from the equation of the line:
N_endpoints = line_coeffs(1)*t_endpoints + line_coeffs(2);
plot(t_endpoints,N_endpoints)
hold off
% The rate constant is -(slope).
k = -line_coeffs(1)