Skip to content

Instantly share code, notes, and snippets.

@profConradi
Last active May 6, 2024 20:58
Show Gist options
  • Save profConradi/202139aa56af21eef1dad1c155a83421 to your computer and use it in GitHub Desktop.
Save profConradi/202139aa56af21eef1dad1c155a83421 to your computer and use it in GitHub Desktop.
Chladni Square Plate Normal Modes Simulation
def energy(z, n, m, L):
return np.cos(n* np.pi *np.real(z) / L) *np.cos(m *np.pi*np.imag(z) / L) - np.cos(m*np.pi*np.real(z)/ L) *np.cos(n* np.pi *np.imag(z) / L)
class ChladniPlate:
def __init__(self, n, m, L=1, n_particles=10000):
self.L = L
self.n_particles = n_particles
self.n = n
self.m = m
self.positions = np.random.uniform(0, self.L, (self.n_particles,))+1.j*np.random.uniform(0, self.L, (self.n_particles,))
self.history = [self.positions]
def metropolis_step(self, z, eps, beta):
energies = energy(z, self.n, self.m, self.L)**2
delta = np.random.uniform(-eps,eps,(self.n_particles,)) + 1.j*np.random.uniform(-eps,eps,(self.n_particles,))
new_energies = energy(z+delta, self.n, self.m, self.L)**2
r = np.exp(-beta*(new_energies-energies))
r[r>1]=1
x = np.random.uniform(0., 1., self.n_particles)
new_positions = z+delta
update = (x<r) & (np.real(new_positions)>0) & (np.real(new_positions)<self.L) & (np.imag(new_positions)>0) & (np.imag(new_positions)<self.L)
new_positions[np.logical_not(update)] = z[np.logical_not(update)]
return new_positions
def simulate(self, n_iter, eps, beta):
for _ in range(1, n_iter):
self.history.append(self.metropolis_step(self.history[-1], eps, beta))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment