plotting $F=K \frac{q_1 q_2}{r^3}$ in python

I tryin to plot the position of electron in the hidrogen atom by Coulomb's law, $$F=K \frac{q_1 q_2}{r^3}$$

This is mi code

    import numpy as np
import matplotlib.pyplot as plt

def coulomb(re, v0):

#Set parameters:
N = 500
dt = 2.2222222 / N  # Time Step:

#Create an array, for all variables, of size N with all entries equal to zero:
xe = np.zeros((N,))
ye = np.zeros((N,))
vxe = np.zeros((N,))
vye = np.zeros((N,))

# Initial Conditions:
xe[0] = re                   # (x0 = r, y0 = 0)
vye[0] = v0                    #units in m/s^2

#Implement Verlet Algorithm:
for k in range(0, N-1):
re = (xe[k]*2+ye[k]*2)*0.5
vxe[k+1] = vxe[k] - ((mu * xe[k]) / (re**3)) * dt
xe [k+1] = xe[k] + vxe[k+1]*dt
vye[k+1] = vye[k] - ((mu * ye[k]) / (re**3)) * dt
ye [k+1] = ye[k] + vye[k+1]*dt

#Plot:
xi = plt.plot(xe, ye, 'go', markersize = 1)
plt.plot(0,0,'yo')                  # yellow marker
plt.plot(xe[0],0,'bo')  # dark blue marker
plt.axis('equal')
plt.xlabel ('x')
plt.ylabel ('y')

return xi, xe, ye

# average distance electron-nucleus in meter
r = 5.1e-11
k = 8.9e9
mu = r*3 * 4 * k *np.pi*2  # coulomb parameter

coulomb(r, np.sqrt(mu / r));


This works well, I've try to changing the value of Permittivity in l, but the graph becomes just 2 dots (the values of xe and vye remain constant) instead of a circular orbit and I don't know why this is? this happens

• Where do you compute mu? It should be mu = k*el^2/me where el is the elementary charge and me the electron mass? What variable is or is influenced by "Permittivity in l"? – Lutz Lehmann Jan 30 at 11:09
• oh sorry man, i just realize that i put the wrong code, i upload the correct version – Rei D Gar Jan 30 at 18:21

The Kepler laws, or simple insertion of $$z(t)=re^{iωt}$$, tell that for a circular orbit of radius $$r$$ in the central field $$\ddot z=-\mu\frac{z}{|z|^3}$$ the angular speed $$\omega$$ is given by $$r^3ω^2=μ.$$ That is how you get the initial velocity $$(0,\sqrt{\mu/r})$$ at the point $$(r,0)$$. The period of the orbit is $$T=\frac{2\pi}{ω}=2π\sqrt{r^3/μ}.$$ Inserting $$\mu=k·4\pi^2·r^3$$ as coded gives $$T=k^{-1/2}=1.06·10^{-5}$$. Your sampling rate is just not high enough. Or you need to reduce the integration time to about $$2.2·10^{-5}$$ to integrate over two periods.
But most importantly, check your code and formula everywhere that you use the ** operator where you want to compute a power, at the moment there are many places where you have just one multiplication sign.
re = (xe[k]**2+ye[k]**2)**0.5