I have a problem updating the "reversalsize" list. It should update every time with "reversal". I need to make a graph to check the power law behavior. It's a simple cellular automata simulating an avalanche behavior, along a channel. I always get the same value in the list, instead it should be an increasing number
reversal=0
reversalsize=[]
def avalanche(M,KH,pmt):
p_scale_factor = 5
global reversal
M_av = M.copy()
M_av_app = M.copy()
for i in range(nx-1):
if i<nx/3:
for j in range(ny-1):
ix = i+1
iy = j+1
p_drop = np.tanh((M[ix,iy]-pmt)/p_scale_factor)
rand_drop = np.random.random()
if rand_drop > 1-p_drop:
if j==0:
if KH[ix,iy] > pmt:
av_mass = (M[ix,iy]-pmt)*ppe
M_av[ix+1,iy] = M_av_app[ix+1,iy]+3*av_mass/4
M_av[ix+1,iy+1] = M_av_app[ix+1,iy+1]+av_mass/4
M_av[ix,iy] = M_av_app[ix,iy]-av_mass
elif 0<j<ny-3:
if KH[ix,iy] > pmt:
av_mass = (M[ix,iy]-pmt)*ppe
M_av[ix+1,iy-1] = M_av_app[ix+1,iy-1]+av_mass/4
M_av[ix+1,iy] = M_av_app[ix+1,iy]+av_mass/2
M_av[ix+1,iy+1] = M_av_app[ix+1,iy+1]+av_mass/4
M_av[ix,iy] = M_av_app[ix,iy]-av_mass
elif j==ny-3:
if KH[ix,iy] > pmt:
av_mass = (M[ix,iy]-pmt)*ppe
M_av[ix+1,iy] = M_av_app[ix+1,iy]+3*av_mass/4
M_av[ix+1,iy-1] = M_av_app[ix+1,iy-1]+av_mass/4
M_av[ix,iy] = M_av_app[ix,iy]-av_mass
M_av_app = M_av
reversal +=1
#print(reversal)
else:
ilim1 = i-nx/3
ilim2 = nx-(i-nx/3)
for j in range(ny-1):
ix = i+1
iy = j+1
p_drop = np.tanh((M[ix,iy]-pmt)/p_scale_factor)
rand_drop = np.random.random()
if rand_drop > 1-p_drop:
if j<ilim1:
pass
if j==ilim1:
if KH[ix,iy] > pmt:
av_mass = (M[ix,iy]-pmt)*ppe
M_av[ix+1,iy+1] = M_av_app[ix+1,iy+1]+av_mass
M_av[ix,iy] = M_av_app[ix,iy]-av_mass
elif ilim1<j<ilim2:
if KH[ix,iy] > pmt:
av_mass = (M[ix,iy]-pmt)*ppe
M_av[ix+1,iy-1] = M_av_app[ix+1,iy-1]+av_mass/4
M_av[ix+1,iy] = M_av_app[ix+1,iy]+av_mass/2
M_av[ix+1,iy+1] = M_av_app[ix+1,iy+1]+av_mass/4
M_av[ix,iy] = M_av_app[ix,iy]-av_mass
elif j==ilim2:
if KH[ix,iy] > pmt:
av_mass = (M[ix,iy]-pmt)*ppe
M_av[ix+1,iy-1] = M_av_app[ix+1,iy-1]+av_mass
M_av[ix,iy] = M_av_app[ix,iy]-av_mass
M_av_app = M_av
reversal +=1
#print(reversal)
return M_av
for i in range(n_iter):
reversal=0
avalanche(M,KH,pmt)
reversalsize.append(reversal)
print(reversalsize)
x, y = np.histogram(reversalsize, 1000)
plt.figure(figsize=(10,7))
plt.clf()
plt.loglog(y[0:-1],x, 'r.')
plt.title("Avalanche Size Distribution", fontsize=14)
plt.xlabel(r"$\log$" + "(Avalanche Size)", fontsize=12)
plt.ylabel(r"$\log$" + "(Frequency)", fontsize=12)
plt.show()
Every time you iterate, you're re-setting reversal to 0. This is why it's not increasing:
for i in range(n_iter):
reversal=0 # Delete this line