from random import random
from numpy import arange
from pylab import plot,xlabel,ylabel,show

#Constantes
NTl = 1000       #numero inicial de thalios
NPb = 0          #numero inicial de plomos
tau = 3.053*60   #vida media del thalio en seg.
th = 1.0         #tamano de paso en el tiempo seg.
p = 1-2**(-th/tau) #probabilidad de que decaiga un Tl
tmax = 1000       #tiempo total seg
    
#los puntos a graficar
tpuntos = arange(0.0,tmax,th)
Tlpuntos = []
Pbpuntos = []

#ciclo principal
for t in tpuntos:
    Tlpuntos.append(NTl)
    Pbpuntos.append(NPb)

    #calculando el numero de nucleos que decayeron
    decay=0
    for i in range(NTl):
        if random()< p:
          decay+=1
    NTl -= decay
    NPb += decay

#La grafica
plot(tpuntos,Tlpuntos)
plot(tpuntos,Pbpuntos)
xlabel("Tiempo")
ylabel("Numero de nucleos")
show()