Chladni Square Plate Normal Modes Simulation

chladni.py

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))