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

Question

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

Answer 1

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.

lines to repair:

re = (xe[k]**2+ye[k]**2)**0.5
mu = k * r**3 * 4 * np.pi**2  # coulomb parameter