Decay

% Cantidad inicial de partículas:
n = 100;
% Tasa de decaimiento:
lambda = 1;
% Tiempos aleatorios de decaimiento:
T1 = rande(n, 1) / lambda;
T2 = rande(n, 1) / lambda;
T3 = rande(n, 1) / lambda;

% Tiempos auxiliares para graficar:
t = linspace(0, max(max([T1 T2 T3]))*1.1, 1000);
% Proceso:
N1 = sum(T1 >= t);
N2 = sum(T2 >= t);
N3 = sum(T3 >= t);

% Media:
p = exp(-lambda*t);
u = n*p;
% Desvío:
d = sqrt( u .* (1-p));

% Probabilidades:
P = [];
for i = 0:n
    P = [P; (p.^i).*(1-p).^(n-i)*bincoeff(n, i)];
end
% P = flipud(P);


figure(1);
colormap(summer);
brighten(0.5);

imagesc(t, 0:n, P, [0 1]);
axis([0 t(end) 0 n], 'xy');
hold on;
plot(   t, [N1; N2; N3],
        t, u, '-g');
legend( 'Proceso 1',
        'Proceso 2',
        'Proceso 3',
        'Media');

colorbar()
xlabel('Tiempo [s]');
ylabel('Cantidad de particulas vivas');
title('Decaimiento');
print('decaimiento.png', '-dpng', '-S1024, 600');
pause