import numpy as np
import matplotlib.pyplot as plt

ZZ= 7  #nitrógeno
E0 = 20*10**9 #eV
EC = 600*10**6/ZZ #eV

def rand():return(np.random.uniform(0,1))

E=[E0]
X=[0]
N=[1]
T=[]

tmax = int(np.log(E0/EC)/np.log(2))

for t in range(0,tmax):
    X.append(t+1)

    #s=N[t]
    for i in range(len(E)):
        x=rand()
        e=0.5*x*E[i]
        e1=0.5*(1-x)*E[i]
        if e>EC:
            T.append(e)
        if e1>EC:
            T.append(e1)
        prof=len(T)
        N.append(prof)
        E=T
        T=[]
        
plt.plot(X,np.log(N))
plt.xlabel('Longitudes de radiación')
plt.ylabel('Número de partículas')
plt.title('Cascada electromagnética')
plt.grid(True)