% Matlab code to simulate a fluorescence/phosphorescence system
% Marc R. Roussel

% Start with a clean slate
clear
close all

% Transition rates in s^{-1}
wvr = 1e11;
wf = 1e8;
wisc = 1e8;
wph = 1e-2;

% Initial number of molecules
Nmol = 10000;

% Population vector ordering: [S*,S,T*,T]
N = zeros(4,1);
N(1) = Nmol;
t = 0;

% Counters for the numbers of fluorescence and phosphorescence events.
nfl = 0;
nph = 0;

% Main loop: continue until we run out of molecules.
while (Nmol > 0)
    % Calculate the propensities. There is no super-clever way to do this.
    % Ordering of the reactions: [S*->S, S fluor, ISC, T^*->T, phosphor]
    a(1) = wvr*N(1);
    a(2) = wf*N(2);
    a(3) = wisc*N(2);
    a(4) = wvr*N(3);
    a(5) = wph*N(4);
    
    % KMC algorithm:
    r = rand(2,1);
    a_sums = cumsum(a);
    atot = a_sums(5);
    Delta_t = -log(r(1))/atot;
    k = find(a_sums>r(2)*atot,1,'first');
    
    % Time at next step:
    t = t + Delta_t;
    
    % Take the appropriate action based on the reaction chosen:
    if (k == 1)     % S*->S
        N(1) = N(1) - 1;
        N(2) = N(2) + 1;
    elseif (k==2)   % fluorescence
        N(2) = N(2) - 1;
        Nmol = Nmol - 1;
        nfl = nfl + 1;
        tfl(nfl) = t;
    elseif (k==3)   % ISC
        N(2) = N(2) - 1;
        N(3) = N(3) + 1;
    elseif (k==4)   % T^*->T
        N(3) = N(3) - 1;
        N(4) = N(4) + 1;
    elseif (k==5)   % phosphorescence
        N(4) = N(4) - 1;
        Nmol = Nmol - 1;
        nph = nph + 1;
        tph(nph) = t;
    else    % Defensive programming: we should never get here.
        error('Bad reaction')
    end
end

% Generate histograms of tfl and tph
figure
hist(tfl,20)
xlabel('t/s')
ylabel('Events')
title('Distribution of fluorescence times')

figure
hist(tph,20)
xlabel('t/s')
ylabel('Events')
title('Distribution of phosphorescence times')

% Calculate the fluorescence and phosphorescence decay times:
tfl_decay = mean(tfl)
tph_decay = mean(tph)