% Simulate IVR for a system with two vibrational modes related by
% omega1 = 2omega2, and where reaction occurs if omega2 >= nquant2_min.

% General simulation parameters
% Transition rate in s^{-1}. Can be adjusted to fit experimental data.
w = 1e15;  
% Minimum number of quanta required in mode 2 for reaction:
nquant2_min = 10;    
nsims = 10000;  % Number of simulations at each energy
% Create a vector to hold the values of t from each simulation
% at any given value of nquant2.
t2K = zeros(nsims,1);
% Maximum number of omega2 quanta
nquant2_max = 50;
% nq2count is a counter for the number of nquant2 values 
nq2count = 0;

% Loop over nquant2 values
for nquant2=nquant2_min:nquant2_max
    nq2count = nq2count + 1;
    % Generate nsims simulations.
    for sim=1:nsims
        % Generate a random initial state with nquant2 quanta.
        v1max = floor(nquant2/2);
        v1 = randi([0,v1max]);
        v2 = nquant2 - 2*v1;
        
        % Main simulation loop
        t = 0;
        wtot = w*v1max;
        while (v2 < nquant2_min)
            % Generate the random Delta t.
            r = rand(1);
            Delta_t = -log(r)/wtot;
            % Generate the new values of v1.
            % The while statement makes sure that there is an actual jump,
            % i.e. that the new value of v1 is different from the current 
            % value.
            new_v1 = v1;
            while (new_v1 == v1)
                new_v1 = randi([0,v1max]);
            end
            % Update state and time.
            v1 = new_v1;
            v2 = nquant2 - 2*v1;
            t = t + Delta_t;
        end
        
        % The value of t is how long it took for nquant_min quanta
        % to show up in mode 2. We store this value for later averaging.
        t2K(sim) = t;
    end
    
    % Calculate the rate constant and its error,
    % and store the values for plotting.
    nq2(nq2count) = nquant2;
    k2K(nq2count) = 1.0/mean(t2K);
    stderr_t2K = std(t2K)/sqrt(nsims);
    errk2K(nq2count) = stderr_t2K*k2K(nq2count)^2;
end

figure('DefaultAxesFontSize',18)
errorbar(nq2,k2K,errk2K,'o')
xlabel('Number of quanta')
ylabel('$k_{2K}/$s$^{-1}$','interpreter','latex')