% Control3_4.m
% Solution to problem 3 of control 3
% JMS, May.2024

% Set parameters
clear
tau = 1.5e3;  % mean life time
N = 1e4;      % number of particles

% Find random times for nuclei desintegration
F = rand(1,N);       % probability distribution values
t = -tau*log(1-F);   % desintegration times

% Find required magnitudes
t = sort(t);              % times t are now in increasing order
tP05 = t(N/2);            % (N/2)'th desintegration time (N must be even)
PtGTtau = sum(t>tau)/N;                          % fraction of times t>tau
PtLTtau3GT5dt3 = ( sum(t>tau/3 & t<5*tau/3) )/N; % fraction with tau/3<t<5*tau/3
fprintf('     t for (P=0.5) = %4.0f s\n',tP05)
fprintf('          P(t>tau) = %5.3f\n',PtGTtau)
fprintf('P(tau/3<t<5*tau/3) = %5.3f\n',PtLTtau3GT5dt3)

% Desintegration time intervals
dt = diff(t);             % time intervals between succesive desintegrations
H = (N:-1:2).*dt/tau;     % normalized intervals
Hmean = meanValue(H);
DH = standDev(H);

% Plot histogram of H
figure(2)
histogram(H,50), hold on
plot([Hmean,Hmean],[0,1e3],'r','LineWidth',3)
plot([Hmean-DH,Hmean-DH],[0,1e3],'g','LineWidth',2)
plot([Hmean+DH,Hmean+DH],[0,1e3],'g','LineWidth',2)
hold off
xlabel('H')
ylabel('Frequency')
grid on
set(gca,'FontSize',15)
