Commit 091ba0a1 authored by Luuk Balkenende's avatar Luuk Balkenende

Fixed small error in bead calculation

parent 3610c510
......@@ -3,33 +3,39 @@ import sys
def add_bead(pos, n, T_dim, j, w_jl, w_l):
# Standard phi, KAN NAAR MAIN.PY
# Standard phi, KAN NAAR MAIN.PY
phi_0 = np.linspace(0, 2 * np.pi - (2 / n) * np.pi, n)
# Preallocate arrays
E = np.zeros(n)
# Calculate energy of different angles. DEZE FOR LOOP KAN WEG.
for i in range(n):
phi = phi_0 + (1-2*np.random.rand(n))*np.pi/n
pos_new = np.zeros(2)
pos_new[0] = pos[j-1, 0] + np.cos(phi[i])
pos_new[1] = pos[j-1, 1] + np.sin(phi[i])
rx = pos_new[0] - pos[0:j, 0]
ry = pos_new[1] - pos[0:j, 1]
r = np.sqrt(rx ** 2 + ry ** 2)
r[r == 0] = np.nan
E[i] = np.sum(4 * ((1 / r) ** 12 - (1 / r) ** 6))
#print('E = ', E)
# Calculate energy of different angles.
phi = phi_0 + (1 - 2 * np.random.rand(n)) * np.pi / n
pos_new = np.zeros([n, 2])
pos_new[:, 0] = pos[j - 1, 0] + np.cos(phi)
pos_new[:, 1] = pos[j - 1, 1] + np.sin(phi)
rx = np.zeros([n, j])
ry = np.zeros([n, j])
for i in range(n):
rx[i, :] = pos_new[i, 0] - pos[0:j, 0]
ry[i, :] = pos_new[i, 1] - pos[0:j, 1]
print('rx = ', rx)
r = np.sqrt(rx ** 2 + ry ** 2)
r[r == 0] = np.nan
E = np.sum(4 * ((1 / r) ** 12 - (1 / r) ** 6), axis=1)
print('E = ', E)
# for i in range(n):
# phi = phi_0 + (1 - 2 * np.random.rand(n)) * np.pi / n
# pos_new = np.zeros(2)
# pos_new[0] = pos[j - 1, 0] + np.cos(phi[i])
# pos_new[1] = pos[j - 1, 1] + np.sin(phi[i])
# rx = pos_new[0] - pos[0:j, 0]
# ry = pos_new[1] - pos[0:j, 1]
# r = np.sqrt(rx ** 2 + ry ** 2)
# r[r == 0] = np.nan
# E[i] = np.sum(4 * ((1 / r) ** 12 - (1 / r) ** 6))
# print('E = ', E)
w = np.exp(-E / T_dim)
#if w == NaN:
print('w = ', w)
#w[np.isnan(w)] = 0
#print('w = ', w)
w[np.isnan(w)] = 0
w_tot = np.sum(w)
if w_tot == 0:
......@@ -38,13 +44,14 @@ def add_bead(pos, n, T_dim, j, w_jl, w_l):
sys.exit()
P = w / w_tot
print('p = ', P)
final_phi = np.random.choice(phi, p=P)
pos[j] = np.array([[pos[j-1, 0] + np.cos(final_phi), pos[j-1, 1] + np.sin(final_phi)]])
pos[j] = np.array([[pos[j - 1, 0] + np.cos(final_phi), pos[j - 1, 1] + np.sin(final_phi)]])
w_idx = np.where(phi == final_phi)
print('w_idx = ', w_idx[0])
w_jl = w_jl*w[w_idx[0]]
w_jl = w_jl * w[w_idx[0]]
w_l = w_l * w_tot
return pos, w_jl, w_l
\ No newline at end of file
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment