%% ELEC 3004 
% Quick & Dirty Code for Question 8
% For a review of Bateman's Equation. See:
% -- http://www.nucleonica.com/wiki/index.php?title=Help:Decay_Engine#Radioactive_Decay_Chains
% -- http://chemistry.sfu.ca/assets/uploads/file/Course%20Materials%2012-1/NUSC%20342/L9.pdf
% -- http://en.wikipedia.org/wiki/Radioactive_decay

%% Have Matlab Solve the ODE
% The solution can be found from the Bateman's Equation
% But this is a "lazy" way to type in all the equations
syms N1(t) N2(t) N3(t) N4(t) l1 l2 l3 l4 N0

q1=dsolve(diff(N1)==-l1*N1, N1(0) == N0);
q2=dsolve(diff(N2)==-l2*N2+l1*q1, N2(0) == 0);
q3=dsolve(diff(N3)==-l3*N3+l2*q2, N3(0) == 0);
q4=dsolve(diff(N4)==l3*q3, N4(0)==0);


%% Handle Normalzied Case (part c)
% Input lamba values
k1=log(2)/36;
k2=log(2)/2.1;
k3=log(2)/(4.8);
k4=0;
N0val = 1;  % value for N1(0) 

% Substitutions
eq1=subs(q1, [l1,l2,l3,l4, N0], [k1,k2,k3,k4, N0val])
eq2=subs(q2, [l1,l2,l3,l4, N0], [k1,k2,k3,k4, N0val])
eq3=subs(q3, [l1,l2,l3,l4, N0], [k1,k2,k3,k4, N0val])
eq4=subs(q4, [l1,l2,l3,l4, N0], [k1,k2,k3,k4,N0val])

% Output in Latex
% (speaking of which -- 
% also check out "publish('elec3004q8graph.m','latex')"
texeq1=latex(simplify(eq1));
texeq2=latex(simplify(eq2));
texeq3=latex(simplify(eq3));
texeq4=latex(simplify(eq4));

% Set graph limits
tmin=0;
tmax=240;
ymin=0;
ymax=1;

% Reset figure
figure(123)
close(123)
figure(123)

% Plot easily via ezplot
hold on;
h1=ezplot(char(eq1), [tmin,tmax,ymin,ymax] );
set(h1, 'Color', 'r', 'LineWidth', 3);

h2=ezplot(char(eq2), [tmin,tmax,ymin,ymax]);
set(h2, 'Color', 'g', 'LineWidth', 2);

h3=ezplot(char(eq3), [tmin,tmax,ymin,ymax]);
set(h3, 'Color', 'b', 'LineWidth', 2);

h4=ezplot(char(eq4), [tmin,tmax,ymin,ymax]);
set(h4, 'Color', 'm', 'LineWidth', 1);

hold off;

legend('^{211}Pb', '^{211}Bi', '^{207}Tl', '^{207}Pb');

title('')
xlabel('Time (min)')
ylabel('Normalized Concentrations')

%% Handle case with 1 kg of atoms 

% N1(0) for 1 kg of Pb-211:
N0val = (1000/211)*(6.022E23);  % (1000 g/211 g/mol)*(Avogadro's Constant)

% Solve for New Concetration
eq4_1kg=subs(q4, [l1,l2,l3,l4, N0], [k1,k2,k3,k4,N0val]);
solve(eq4_1kg == ((900/207)*(6.022E23)), t)

% New lambda terms in seconds
k1sec=log(2)/(36*60);
k2sec=log(2)/(2.1*60);
k3sec=log(2)/((4.8)*60);
k4sec=0;

% New Substitutions (Assuming all nuclei are radioactive -- the curve will
% have the same shape even with an approriate fraction of them)
radiation=l1*q1+l2*q2+l3*q3;
eq_rad=subs(radiation, [l1,l2,l3,l4, N0], [k1sec,k2sec,k3sec,k4sec, N0val]);

syms tm   % time in minutes
eqradmin=subs(eq_rad, t, tm*60)
eradcurie=eqradmin*(1/3.7E10);

figure(234);
ezplot(eradcurie, [0 150])
title('');
xlabel('Time (min)')
ylabel('Radiation (Curie)')

