CO67_2

function [ pop_time_data ] = simulate_decay ( init_pop , half_lives , t)
% Author: Sneha Ramshanker  , Date: ??/??/????
% Simulating decay of D elements given the initial population, half lives, and the time grid.
% Input/output vectors in this function are sorted with the first element decays to the second
% element, then to the third element, and so on until the last element.
% So the last element is the stable element.
% Input:
% ? init_pop: D?elements vector of the initial population of each element, sorted
% such that the stable element is at the last.
% ? half_lives: (D?1)?elements vector of half lives of each element.
% ? t: N?elements vector of the time grid where the population is returned.
%
% Output:
% ? pop_time_data: (N x D) matrix that contains the population of each element in every time step.
%
% Constraints:
% ? The half life of an element can be very small, so in this case, the function must
% make an appropriate adjustment/approximation to handle this.
%
% Example use:
% >> init_pop = [1000, 0, 0, 0, 0];
% >> half_lives = [2.5e5, 8e4, 1.62e3, 4/365];
% >> t = linspace(0, 1e6, 1000);
% >> pop = simulate_decay(init_pop, half_lives, t);
% >> plot(t, pop(:,1)); % plot the population of the 1st element vs time

%Defining initial constraints
    dt = t(2)-t(1);                     %time interval for the integration
    n_points = length(t);               %number of time intervals
    pop_time_data = [init_pop];         %initializing output vector
    pop = init_pop;                     %Initializing current population vector
    N = sum(init_pop);                  %Total number of atoms

    half_lives = half_lives(1:3);       %accounting for stable conditions

    l =  log(2)./half_lives;            %conversion to activity

    for i = 1:n_points-1                %iterating through time intervals
        for j = 1: length(l)            %iterating through the different particles

            %prob = l*dt;
            prob = (1-exp(-l*dt));      %setting probability of decay

            for k = 1:pop(j)            % Iterating through each atom in each particle type
                if prob(j)>rand()       % Applying monte-carlo integration
                    if pop(j)>0
                        pop(j) = pop(j)-1;
                        pop(j+1) = pop(j+1)+1;
                    end
                end
            end
        end
        pop(4) = N - (pop(1)+pop(2)+pop(3));
        pop_time_data = [pop_time_data;pop]; %Storing population values
    end

end