% 10.34 - Fall 2006
% HW Set #11 - Kinetic Monte Carlo (Gillespie's Algorithm)
% Rob Ashcraft
% 12/10/2006

% This will generate the data necessary to do part 3.  It calculates the
% rate of failure for D = 250 nm for various k-values to determine the
% dependence of rate of failure on k.  

% This is very similar to "p1_mian.m".  see that code for extensive
% commenting.  


function p1_main_partc

clear; clc;
warning off all;
Nav = 6.022e23;

time_0 = cputime;

line_style_vec = {'o','-','.','--','d',':','^','-.','x','+','s'};
color_vec = {'k','b','r','g','k','b','r','g','k','b','r','g'};

D_vec = [30; 50; 100; 250; 500; 1000; 2500];
k_vec = [1; 10; 100];

% initial conditions
C_ROO_0 = 0;        % molar
C_ROOH_0 = 1e-6;    % molar
C_RH_0 = 10;        % molar, constant value
C_O2_0 = 1e-4;      % molar, constant value

Temp = 313;            % Kelvin, constant value

C_ROOH_cutoff = 0.010;      % concentration at which to end simulation

disp('Initial number of molecules of: ROO, ROOH, RH, and O2 in droplet of size D');
for i=1:length(D_vec)
    [N_ROO_0,N_ROOH_0,N_RH_0,N_O2_0] = calc_N(C_ROO_0,C_ROOH_0,C_RH_0,C_O2_0,D_vec(i));
    disp(['Droplet size (nm): ',num2str(D_vec(i)),' -- ROO = ',num2str(N_ROO_0),...
        ',  ROOH = ',num2str(N_ROOH_0),',  RH = ',num2str(N_RH_0),',  O2 = ',num2str(N_O2_0)]);
    legend_D{i,1} = ['D = ',num2str(D_vec(i)),' nm'];
end
disp('  ');


% Now work on the Kinetic Monte Carlo part of the problem:
t_per_D = 1800;  % wall time you want to spend per droplet size

% Total sim time at a function of number of steps(N):  Time = m*N
% these values will different on a different computer
m =  [1.85e-5;
      8.88e-5;
      0.0026;
      0.113;
      0.898;
      7.188;
      113.1];

for i = 1:length(D_vec)

    MC_steps_temp(i) = t_per_D/m(i);
    
    places(i) = length(num2str(round(MC_steps_temp(i))));  % number of place left of decimal
    if(places(i) > 2)
        MC_steps(i) = ceil(MC_steps_temp(i)/(10^(places(i)-2))) * 10^(places(i)-2);
    elseif(places(i) == 2)
        MC_steps(i) = ceil(MC_steps_temp(i)/10) * 10;
    else
        MC_steps(i) = ceil(MC_steps_temp(i));
    end
end

k_vec_partC = logspace(0,3,13);

for j=1:length(k_vec_partC);    % for k = 10
for i = 4:4%length(D_vec)
    time_in = cputime;
    
    disp(['Now running MC simulation for D = ',num2str(D_vec(i)),' nm']);
    started_dead(j) = 0;
    started_but_died(j) = 0;
    went_too_long(j) = 0;

    tc = 1;
    [N_ROO,N_ROOH,N_RH_0,N_O2_0] = calc_N(C_ROO_0,C_ROOH_0,C_RH_0,C_O2_0,D_vec(i));

    % defined vector sizes to make the code run faster
    if(N_ROOH > 1)
        N_ROOH_temp{j} = zeros(MC_steps(i),1);
        time_fail_temp{j} = zeros(MC_steps(i),1);
        rate_fail_temp{j} = zeros(MC_steps(i),1);
    else        
        N_ROOH_temp{j} = zeros(MC_steps(i)*N_ROOH*1.25,1)*NaN;
        time_fail_temp{j} = zeros(MC_steps(i)*N_ROOH*1.25,1)*NaN;
        rate_fail_temp{j} = zeros(MC_steps(i)*N_ROOH*1.25,1)*NaN;
    end
    
    Vol = pi/6*D_vec(i)^3;  % volume in nm^3
    
    sizect = 0;
    for mc = 1:MC_steps(i)

        [N_ROO,N_ROOH,N_RH_0,N_O2_0] = calc_N(C_ROO_0,C_ROOH_0,C_RH_0,C_O2_0,D_vec(i));
        
        % since the N_0's calculated from the concentrations will be non-integer,
        % we need to round them to the nearest integer based on a probability.
        % This is accomplished with the following lines:
        N_ROO = round(N_ROO + (rand-.5));
        N_ROOH = round(N_ROOH + (rand-.5));
        N_RH_0 = round(N_RH_0 + (rand-.5));
        N_O2_0 = round(N_O2_0 + (rand-.5));
        
        if(N_ROOH == 0)
            %rates = [0;0;0;0];   tau = inf;
            started_dead(j) = started_dead(j)+1;
            mctime = NaN;
            %N_ROOH_final{i}(mc,1) = NaN;
            %time_fail_final{i}(mc,1) = mctime;
            
        else
            
            rates(1) = (1e15*exp(-15000/Temp)*N_ROOH/Vol)*Vol;                           % molecules/nm^3-s      ROOH + 2_RH + 2_O2 ==> 2_ROO + ROH + H2O
            rates(2) = (1e9*(1e24/6.022e23)*exp(-6000/Temp)*N_ROO/Vol*N_RH_0/Vol)*Vol;     % molecules/nm^3-s      ROO + RH + O2 ==> ROOH + ROO
            rates(3) = (1e6*(1e24/6.022e23)*(N_ROO/Vol)*((N_ROO-1)/Vol))*Vol;                       % molecules/nm^3-s      2_ROO ==> ROH + aldehyde + O2
            rates(4) = (0.10/D_vec(i)*k_vec_partC(j)*N_ROO/Vol)*Vol;                      % D in nm       % molecules/nm^3-s      ROO + X(aq) ==> ROOH + X_dot(aq)

            tau = 1/(rates(1) + rates(2) + rates(3) + rates(4));
            
            mctime = 0;
            
            C_ROOH = C_ROOH_0;
            dead = 0;
            while(C_ROOH <= C_ROOH_cutoff  &  mctime <= 1e8  &  dead ~= 1  &  tau < inf)
                
                
                if(N_ROO + N_ROOH == 0)
                    dead = 1;
                    N_ROOH = N_ROOH;
                    N_ROO = N_ROO;
                elseif(N_ROO + N_ROOH ~= 0)            

                    rates(1) = (1e15*exp(-15000/Temp)*N_ROOH/Vol)*Vol;                           % molecules/nm^3-s      ROOH + 2_RH + 2_O2 ==> 2_ROO + ROH + H2O
                    rates(2) = (1e9*(1e24/6.022e23)*exp(-6000/Temp)*N_ROO/Vol*N_RH_0/Vol)*Vol;     % molecules/nm^3-s      ROO + RH + O2 ==> ROOH + ROO
                    rates(3) = (1e6*(1e24/6.022e23)*(N_ROO/Vol)*((N_ROO-1)/Vol))*Vol;                       % molecules/nm^3-s      2_ROO ==> ROH + aldehyde + O2
                    rates(4) = (0.10/D_vec(i)*k_vec_partC(j)*N_ROO/Vol)*Vol;                      % D in nm       % molecules/nm^3-s      ROO + X(aq) ==> ROOH + X_dot(aq)

                    tau = 1/(rates(1) + rates(2) + rates(3) + rates(4));
                    
                    mctime = mctime - tau*log(rand);
                    
                    rand2 = rand;
                    %determine the fate next reaction that happens
                    if(rand2 <= tau*rates(1))         % ROOH + 2_RH + 2_O2 ==> 2_ROO + ROH + H2O
                        N_ROOH = N_ROOH - 1;
                        N_ROO = N_ROO + 2;
                    elseif(rand2 <= tau*(rates(1) + rates(2)))    % ROO + RH + O2 ==> ROOH + ROO
                        N_ROO = N_ROO - 0;
                        N_ROOH = N_ROOH + 1;
                    elseif(rand2 <= tau*(rates(1) + rates(2) + rates(3)))    % 2_ROO ==> ROH + aldehyde + O2
                        N_ROO = N_ROO - 2;
                        N_ROOH = N_ROOH;
                    elseif(rand2 <= tau*(rates(1) + rates(2) + rates(3) + rates(4)))    % ROO + X(aq) ==> ROOH + X_dot(aq)
                        N_ROO = N_ROO - 1;
                        N_ROOH = N_ROOH + 1;
                    end

                    %[rates,tau] = rate_calc_N(N_ROO,N_ROOH,N_RH_0,D_vec(i),k_vec(j),Temp);
                end

                C_ROOH = N_ROOH/((Nav*pi/6*D_vec(i)^3)/1e24);

            end  % while loop

            if(N_ROOH == 0 & mctime > 0)% & mctime <= 1e7)  % started by died
                %time_fail_final{i}(sizect,1) = NaN;
                started_but_died(j) = started_but_died(j) + 1;
            elseif(C_ROOH < 0.01 & mctime > 1e7)        % went too long
                %time_fail_final{i}(sizect,1) = 1e8;
                went_too_long(j) = went_too_long(j) + 1;
            else
                sizect = sizect + 1;
                time_fail_temp{j}(sizect,1) = mctime;
                rate_fail_temp{j}(sizect,1) = 1/mctime;
                N_ROOH_temp{j}(sizect,1) = N_ROOH;
            end

            %N_ROOH_final{i}(mc,1) = N_ROOH;
        
        end   % end else N_ROOH ~= 0 statement
        
        if(places(i) == 1)
            disp(['Step Number ',num2str(mc),' of ',num2str(MC_steps(i)),';  Current wall-time = ',num2str(cputime-time_in),' seconds']);            
        else
            if(rem(mc,MC_steps(i)/10) == 0)
                disp(['Step Number ',num2str(mc),' of ',num2str(MC_steps(i)),';  Current wall-time = ',num2str(cputime-time_in),' seconds']);
            end
        end
        
    end   % end iterations of N_droplets

    max_not_nan = max(find(isnan(N_ROOH_temp{j})==0));
    
    N_ROOH_final{j} = N_ROOH_temp{j}(1:max_not_nan);
    time_fail_final_250{j} = time_fail_temp{j}(1:max_not_nan);
    rate_fail_final_250{j} = rate_fail_temp{j}(1:max_not_nan);
    
    avg_time(j) = exp(mean(log(time_fail_final_250{j})));
    
    
    %Num_died(i) = length(find(N_ROOH_final{i}(:) == 0));
    alive(j) = MC_steps(i) - started_dead(j);
    disp('  ');
    disp(['Number of steps that reached critical C_ROOH: ',num2str(length(N_ROOH_final{j}(:))),' of ',num2str(MC_steps(i)),' steps']);
    disp(['Went too long (but did not die): ',num2str(went_too_long(j)),' of ',num2str(alive(j)),' steps']);
    disp(['Number of droplets that died (but started alive): ',num2str(started_but_died(j)),' of ',num2str(alive(j))]);
    disp(['Fraction of droplets that started dead: ',num2str(started_dead(j)),' of ',num2str(MC_steps(i)),' (',num2str(started_dead(j)/MC_steps(i)*100),'%)']);
    disp('  ');
    disp('  ');
    pause(1);
end
end % loop over j (k values)

plot(k_vec_partC,avg_time,'o');

time_total = cputime - time_0;
disp(['Total wall-time = ',num2str(time_total/3600),' hours']);

save('workspace_whatever_part3');

%========================================
function [rates,tau] = rate_calc_N(N_ROO,N_ROOH,N_RH,D,k,T)
% rates are returned in molecules/nm^3-sec

Vol = pi/6*D^3;

rates(1) = (1e15*exp(-15000/T)*N_ROOH/Vol)*Vol;                           % molecules/nm^3-s      ROOH + 2_RH + 2_O2 ==> 2_ROO + ROH + H2O
rates(2) = (1e9*(1e24/6.022e23)*exp(-6000/T)*N_ROO/Vol*N_RH/Vol)*Vol;     % molecules/nm^3-s      ROO + RH + O2 ==> ROOH + ROO
rates(3) = (1e6*(1e24/6.022e23)*(N_ROO/Vol)*((N_ROO-1)/Vol))*Vol;                       % molecules/nm^3-s      2_ROO ==> ROH + aldehyde + O2
rates(4) = (0.10/D*k*N_ROO/Vol)*Vol;                      % D in nm       % molecules/nm^3-s      ROO + X(aq) ==> ROOH + X_dot(aq)
% 
% if(N_ROO < 2)
%     rates(3) = 0;
% end

tau = 1/sum(rates);

return;

%========================================
function dCdt = continuum_solver(time,C,T,D,k,C_RH,C_O2)

C_ROO = C(1);
C_ROOH = C(2);

rate1 = 1e15*exp(-15000/T)*C_ROOH;      % mole/L-s      ROOH + 2_RH + 2_O2 ==> 2_ROO + ROH + H2O
rate2 = 1e9*exp(-6000/T)*C_ROO*C_RH;    % mole/L-s      ROO + RH + O2 ==> ROOH + ROO
rate3 = 1e6*C_ROO^2;                    % mole/L-s      2_ROO ==> ROH + aldehyde + O2
rate4 = 0.10/D*k*C_ROO;    % D in nm    % mole/L-s      ROO + X(aq) ==> ROOH + X_dot(aq)

dROO_dt = 2*rate1 - 2*rate3 - rate4;
dROOH_dt = -rate1 + rate2 + rate4;

dCdt = [dROO_dt; dROOH_dt];
return;


%=========================================
function [N_ROO,N_ROOH,N_RH,N_O2] = calc_N(C_ROO,C_ROOH,C_RH,C_O2,D)

% diameter should be in nanometers and the concentration in moles/liter
Nav = 6.022e23;
% number of molecules in the droplet: N = C*Nav*pi/6*D^3*(1 dm^3/10^24 nm^3)
N_ROO = (C_ROO*Nav*pi/6*D^3)/1e24;
N_ROOH = (C_ROOH*Nav*pi/6*D^3)/1e24;
N_RH = (C_RH*Nav*pi/6*D^3)/1e24;
N_O2 = (C_O2*Nav*pi/6*D^3)/1e24;

return;