Project 02: What set of upward velocity and position between two protons will yield a stable electron's orbit?

Name: Nam Bui

Introduction:

What set of upward velocity and position between two protons will yield a stable electron's orbit? under these conditions?

1. Two protons are 10 meters apart on x-y plane, located at (-5,0)m and (5,0)m respectively.
2. The initial position of the electron will be varied between -5 to 5 meters to test its orbit's stability.
3. The electron's velocity in the x direction is fixed at 0.
4. The velocity in the y direction will be varied between 5 - 15 m/s to test its orbit's stability. I chose this small range to avoid long code running time and it is just enough to see unbounded and bounded orbits of the electron.
5. The time interval analyzed is 10 seconds only (since I varied 3 sets of initial conditions which might initiate a long running time).
6. Two scenarios that define unstable orbit in this simulation: the electron has enough kinetic energy to escape the system (unbounded) and the electron gets too close to either proton (less than the Bohr's radius: $ 5.3 \times 10^{-3}m $). Otherwise, the orbit is stable.
7. The Bohr's radius is from an excited atom (Reference 3) that has a very large electron-orbiting radius compared to its size. This helps my simulation run smoother since too small radius demands more precision and running time.

The Physics Simulation:

Firstly, the electron is accelerated by the electric fields originating from two protons, which are governed by Coulomb's Law: $$ \mathbf{F} = k \cdot \frac{q_1 q_2}{r^2} \hat{\mathbf{r}} = -k \cdot \frac{q^2}{|\mathbf{r}|^3} $$ The minus sign is due to opposite charge between electron and protons. This total force will be used to define the total acceleration for solve_ivp function. Next, to find out if the electron is unbounded, I will calculate its total mechanical energy by: $$E = PE + KE$$ $$ PE = -\left( \frac{k \cdot q^2}{r_1} \right) - \left( \frac{k \cdot q^2}{r_2} \right) $$

$$ KE = 0.5 \cdot m_e \cdot \left( v_x^2 + v_y^2 \right) $$ Since only the coulomb's force is considered, this is a conservative force field and E is always constant/conserved. Therefore, to define whether the electron is unbounded, I only need to consider its initial total energy. If E >= 0, it is bounded and E < 0, it is unbounded. This is because E = PE + KE and PE is negative, KE is positive. If KE >= |PE|, E>=0. In other words, the electron has enough kinetic energy to escape the system.

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Constants
k = 8.99e9  # N⋅m^2/C^2 
q = 1.6e-19  # C 
m_e = 9.11e-31  # kg (mass of the electron)
bohr_r = 5.3e-3 #m (bohr's radius)

# Nucleus 1 & 2 Position:
x_atom1 = -5  # meter 
y_atom1 = 0  # meter
x_atom2 = 5 #meter
y_atom2 = 0


#Flag stable orbit at the beginning:
unstable_collide = False
unstable_E = False

# Differential equations
def diff_eqns(t, state):
    x, y, vx, vy = state
    r1 = np.sqrt((x - x_atom1)**2 + (y - y_atom1)**2)
    r2 = np.sqrt((x - x_atom2)**2 + (y - y_atom2)**2)
    ax1 = -k * q**2 * (x - x_atom1) / (m_e * r1**3)
    ay1 = -k * q**2 * (y - y_atom1) / (m_e * r1**3)
    ax2 = -k * q**2 * (x - x_atom2) / (m_e * r2**3)
    ay2 = -k * q**2 * (y - y_atom2) / (m_e * r2**3)
    ax = ax1 + ax2
    ay = ay1 + ay2
    return [vx, vy, ax, ay]

#Detection of electron crashes (gets closer to Bohr's radius):
def event1(t, state):
    x, y = state[0], state[1]
    r1 = np.sqrt((x - x_atom1)**2 + (y - y_atom1)**2)
    return r1 - bohr_r  # Event triggers when r1 equals bohr_r
event1.terminal = True
event1.direction = -1   #Only consider distance decreases toward bohr radius

def event2(t, state):
    x, y = state[0], state[1]
    r2 = np.sqrt((x - x_atom2)**2 + (y - y_atom2)**2)
    return r2 - bohr_r  # Event triggers when r2 equals bohr_r
event2.terminal = True
event2.direction = -1   #Only consider distance decreases toward bohr radius

x_list=[]
y_list=[]
z_list=[]
unstable_list=[]

#Loop varying initial velocity:
for v0 in [5,10,15]:
    # Loop testing initial positions:
    for i in np.arange(-4.9, 4.9, 1):
        for j in np.arange(-4.9, 4.9, 1):
            x_list.append(i)
            y_list.append(j)
            z_list.append(v0)
            #Reset the flag after each iteration:
            unstable_collide = False
            unstable_E = False
            state0 = [i, j, 0, v0]  
            #Time settings
            t0, tmax = 0, 10  # seconds
            dt = 0.1
            t_span = (t0, tmax)
            
            # Solve the system
            sol = solve_ivp(diff_eqns, t_span, state0, events=(event1, event2), rtol=1e-6, atol=1e-8, max_step=dt)
            
            #Detection of crashing events:
            if len(sol.t_events[0]) > 0:   #Check if event 1 happens (t_events array is not empty)
                unstable_collide = True            #Flag unstable orbit happened
            elif len(sol.t_events[1]) > 0:   #Check if event 2 happens (t_events array is not empty)
                unstable_collide = True
            else:
                unstable_collide = False
            #Detect when & where electron crashed to proton
            # if len(sol.t_events[0]) > 0:
            #     print(f"Event 1 (Nucleus 1) detected at t={sol.t_events[0][0]:.3e} with position x={sol.y_events[0][0][0]:.3e}, y={sol.y_events[0][0][1]:.3e}")
            # elif len(sol.t_events[1]) > 0:
            #     print(f"Event 2 (Nucleus 2) detected at t={sol.t_events[1][0]:.3e} with position x={sol.y_events[1][0][0]:.3e}, y={sol.y_events[1][0][1]:.3e}")
                

            # Check if the solver succeeded
            if not sol.success:
                continue  # Skip this iteration if the solver failed
            
            #Mechanical Energy:
            x, y, vx, vy = sol.y
            r1 = np.sqrt((x - x_atom1)**2 + (y - y_atom1)**2)
            r2 = np.sqrt((x - x_atom2)**2 + (y - y_atom2)**2)
            PE = -(k * q**2 / r1)  - (k * q**2 / r2)
            KE = 0.5 * m_e * (vx**2 + vy**2)
            E = PE + KE
            
            #For debugging:
            #print(f"x={i}, y={j}, k={k}, PE_min={np.min(PE):.3e}, PE_max={np.max(PE):.3e}, KE_min={np.min(KE):.3e}, KE_max={np.max(KE):.3e}, E_min={np.min(E):.3e}, E_max={np.max(E):.3e}, unstable_E={unstable_E}")
            

            
            #Detection of escaping trajectory:
            if (E >= 0).any():     #E>=0 represents unbounded trajectory
                unstable_E = True
            else:
                unstable_E = False
                
            #Decide if the set of initial conditions yield stable/unstable orbit
            if unstable_E or unstable_collide:
                unstable_list.append(True)
            else:
                unstable_list.append(False)
            # For debugging:
            # print(f"x={i}, y={j}, k={k}, PE_min={np.min(PE):.3e}, PE_max={np.max(PE):.3e}, KE_min={np.min(KE):.3e}, KE_max={np.max(KE):.3e}, E_min={np.min(E):.3e}, E_max={np.max(E):.3e}, unstable_E={unstable_E}")
            # print(unstable_E)
            # print(unstable_collide)


# Assign colors based on the independent 'status' variable
colors = ['blue' if not s else 'red' for s in unstable_list]

# Create a multipanel plot
fig = plt.figure(figsize=(15, 7))

# 3D plot (Panel 1)
ax1 = fig.add_subplot(121, projection='3d')  # 1 row, 2 columns, 1st plot
scatter = ax1.scatter(x_list, y_list, z_list, c=colors, marker='o', s=10)

ax1.set_title("Electron Stability by Initial Conditions", pad=20)
ax1.set_xlabel("Initial x position (m)", labelpad=15)
ax1.set_ylabel("Initial y position (m)", labelpad=15)
ax1.set_zlabel("Initial velocity (m/s)", labelpad=15)

# 2D plot (Panel 2)
ax2 = fig.add_subplot(122)  # 1 row, 2 columns, 2nd plot

# Solve and plot trajectories for two specific initial conditions
initial_conditions = [
    [3.1, 2.1, 0, 10],  # (x, y, vx, vy)
    [3.1, 1.1, 0, 10],  # (x, y, vx, vy)
]

# Time settings
t0 = 0
tmax = 10
dt = 0.1
t_span = (t0, tmax)

for idx, state0 in enumerate(initial_conditions):
    sol = solve_ivp(diff_eqns, t_span, state0, rtol=1e-6, atol=1e-8, max_step=dt)
    x, y = sol.y[0], sol.y[1]
    #Plot the trajectories
    ax2.plot(x, y, label=f"Trajectory {idx + 1}: (x, y, vy) = ({state0[0]}, {state0[1]}, {state0[3]})")

# Proton positions
ax2.scatter([x_atom1, x_atom2], [y_atom1, y_atom2], color='red', s=30, label="Protons")

# Trajectory plot settings:
ax2.set_title("Electron Trajectories for Two Specific Initial Conditions")
ax2.set_xlabel("x position (m)")
ax2.set_ylabel("y position (m)")
ax2.axhline(0, color='gray', linestyle='--', linewidth=0.5)
ax2.axvline(0, color='gray', linestyle='--', linewidth=0.5)
ax2.legend()
ax2.grid(True)
plt.tight_layout()
plt.show()





        

Conclusion

The simulation shows that a stable orbit is supported in multiple regions/sets of initial conditions, for instance, (x,y,vy) = (0, -4.9, 5), (4.9, 4.9, 10), etc. We can also notice that as vy increases, more unstable orbits are detected. This is expected since higher vy means higher total energy and a higher chance of unbounded orbit (E>0).

Indeed, this is restricted to one small region of (x,y,vy) space, in which x and y positions are restricted between -4.9m and 4.1m, and only 3 values of vy are tested (5m/s, 10m/s, 15m/s). The motion of electrons is also restricted to 2-D, which realistically they are in 3-D (or more?). Another restriction of this simulation is the time interval analyzed is 10 seconds, thus, some stable orbits might have not been unstable yet (they might crash to a proton after 10 seconds). But still, the simulation helped answer my research question, in which I can find a set of initial conditions that yield a stable orbit.

To leverage my simulation, I will try to shorten the running time so that I can look at a more dense region of points in the plot. Also, the plot is quite hard to distinguish between points, I will try to find a way to better display the data points.

References

(Each reference in this section should include a citation associated with that reference in the body of your project. You may find it helpful to look at a typical Wikipedia page as an example of how to do this. Make sure you include references for source of datas, literature cited, and any python packages used beyond our standard ones. Additionally, you must also cite the sources for any code that you found on the internet or from peers.)

1. https://chem.libretexts.org/Courses/Chabot_College/Introduction_to_General_Organic_and_Biochemistry/06%3A_Ionic_and_Molecular_Compounds/6.11%3A_Covalent_Bonds#:~:text=The%20bond%20in%20a%20hydrogen,%C3%97%2010%E2%88%9212%20m).
2. https://en.wikipedia.org/wiki/Proton#:~:text=The%20root%20mean%20square%20charge,uncertainty%20of%20%C2%B10.010%20fm.
3. https://electricfusionsystems.com/references/rydberg-physics/

Appendix 1: Code validation

A1.1: Validation of a stable and a unstable orbit:

In this appendix, I will use the code to look closer at two sets of initial conditions in order to validate their result in the previous plot. I picked these two points as shown below:(x,y,vy) = (3.1, 2.1, 10) which is the blue/stable point and (x,y,vy) = (3.1, 1.1, 10) which is red/unstable point.

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Constants
k = 8.99e9  # N⋅m^2/C^2 
q = 1.6e-19  # C 
m_e = 9.11e-31  # kg (mass of the electron)
bohr_r = 5.3e-3 #m (bohr's radius)



# Nucleus 1 & 2 Position:
x_atom1 = -5  # meter 
y_atom1 = 0  # meter
x_atom2 = 5 #meter
y_atom2 = 0

#Flag orbit is stable at beginning:
unstable = False

# Differential equations
def diff_eqns(t, state):
    x, y, vx, vy = state
    r1 = np.sqrt((x - x_atom1)**2 + (y - y_atom1)**2)
    r2 = np.sqrt((x - x_atom2)**2 + (y - y_atom2)**2)
    ax1 = -k * q**2 * (x - x_atom1) / (m_e * r1**3)
    ay1 = -k * q**2 * (y - y_atom1) / (m_e * r1**3)
    ax2 = -k * q**2 * (x - x_atom2) / (m_e * r2**3)
    ay2 = -k * q**2 * (y - y_atom2) / (m_e * r2**3)
    ax = ax1 + ax2
    ay = ay1 + ay2
    return [vx, vy, ax, ay]

def event1(t, state):
    x, y = state[0], state[1]
    r1 = np.sqrt((x - x_atom1)**2 + (y - y_atom1)**2)
    return r1 - bohr_r  # Event triggers when r2 equals bohr_r
event1.terminal = True
event1.direction = -1   #Only consider distance decreases toward bohr radius

def event2(t, state):
    x, y = state[0], state[1]
    r2 = np.sqrt((x - x_atom2)**2 + (y - y_atom2)**2)
    return r2 - bohr_r  # Event triggers when r2 equals bohr_r
event2.terminal = True
event2.direction = -1   #Only consider distance decreases toward bohr radius

# Tested Initial conditions:
for state0 in [[3.1, 1.1, 0, 10],[3.1,2.1,0,10]]: 
    print(f"For (x,y,vx,vy) = {state0}")
    #Time settings
    t0, tmax = 0, 10  # seconds
    dt = 0.1
    t_span = (t0, tmax)
    
    # Solve the system
    sol = solve_ivp(diff_eqns, t_span, state0, events=(event1, event2), rtol=1e-6, atol=1e-8, max_step=dt)
    
    #Detection of crashing events:
    if len(sol.t_events[0]) > 0:   #Check if event 1 happens (t_events array is not empty)
        print(f"Electron approached nucleus 1 at t = {sol.t_events[0][0]:.6f} seconds.")
        print(f"Position: x = {sol.y_events[0][0][0]:.6e}m, y = {sol.y_events[0][0][1]:.6e}m")  #sol.y_events[0] = [[vx,vy,x,y]] of events 1
        unstable = True
    elif len(sol.t_events[1]) > 0:   #Check if event 2 happens (t_events array is not empty)
        print(f"Electron approached nucleus 2 at t = {sol.t_events[1][0]:.6f} seconds.")
        print(f"Position: x = {sol.y_events[1][0][0]:.6e}m, y = {sol.y_events[1][0][1]:.6e}m")
        unstable = True
    else:
        print("Electron did not approach within the Bohr radius of either nucleus.")
        
    #Mechanical Energy:
    x, y, vx, vy = sol.y
    r1 = np.sqrt((x - x_atom1)**2 + (y - y_atom1)**2)
    r2 = np.sqrt((x - x_atom2)**2 + (y - y_atom2)**2)
    PE = -k * q**2 / r1 + -k * q**2 / r2
    KE = 0.5 * m_e * (vx**2 + vy**2)
    E = PE + KE
    
    #Detection of escaping trajectory:
    if (E >= 0).any():
        print("Electron escaped the orbit (E>0)")
        unstable = True
    else:
        print("Electron is bounded (E<= 0)")
    
    # Plot the trajectory
    plt.figure(figsize=(8, 8))
    plt.plot(x_atom1, y_atom1, 'ro', label='Proton 1')  # Nucleus 1 position
    plt.plot(x_atom2, y_atom2, 'ro', label='Proton 2')  # Nucleus 2 position
    plt.plot(state0[0], state0[1], color='blue', label='Electron initial position')  # Electron initial position
    plt.plot(sol.y[0], sol.y[1], 'b-',  marker='o', markersize = 3, markevery=10, label='Electron Trajectory')  # Electron trajectory
    plt.xlabel('x (m)')
    plt.ylabel('y (m)')
    plt.legend()
    plt.grid(True)
    plt.axis('equal')
    plt.title(f"Electron Trajectory Around Two Protons for (x, y, vy) = ({state0[0]}, {state0[1]}, {state0[3]})")
    plt.xlim(0, 10)  
    plt.ylim(-5.5, 6)
    plt.show()

    if state0 == [3.1,2.1,0,10]:
        plt.figure(figsize=(8, 6))
        plt.plot(x_atom1, y_atom1, 'ro', label='Proton 1')  # Nucleus 1 position
        plt.plot(x_atom2, y_atom2, 'ro', label='Proton 2')  # Nucleus 2 position
        plt.plot(state0[0], state0[1], color='blue', label='Electron initial position')  # Electron initial position
        plt.plot(sol.y[0], sol.y[1], 'b-',  marker='o', markersize = 3, markevery=10, label='Electron Trajectory')  # Electron trajectory
        plt.xlabel('x (m)')
        plt.ylabel('y (m)')
        plt.legend()
        plt.grid(True)
        plt.axis('equal')
        plt.title(f"The Zoom into A Region Closed to Proton 2 of the Electron Trajectory for (x, y, vy) = ({state0[0]}, {state0[1]}, {state0[3]})")
        plt.xlim(4, 6)  
        plt.ylim(-1, 1)
        plt.show()
        
    
    # Plot energy
    plt.figure(figsize=(8, 6))
    plt.plot(sol.t, PE, label="Potential Energy (PE)")
    plt.plot(sol.t, KE, label="Kinetic Energy (KE)")
    plt.plot(sol.t, PE + KE, label="Total Energy", linestyle="--")
    plt.xlabel("Time (s)")
    plt.ylabel("Energy (J)")
    plt.legend()
    plt.grid(True)
    plt.title("Energy of the Electron")
    plt.show()

For (x,y,vx,vy) = [3.1, 1.1, 0, 10]
Electron approached nucleus 2 at t = 2.712486 seconds.
Position: x = 5.004690e+00m, y = 2.467710e-03m
Electron is bounded (E<= 0)

For (x,y,vx,vy) = [3.1, 2.1, 0, 10]
Electron did not approach within the Bohr radius of either nucleus.
Electron is bounded (E<= 0)

Indeed, the orbit stability status of these two sets of initial conditions agrees with the finding in the body of the project. Moreover, by this simulation, I could see that the orbit is unstable due to the electron's proximity rather than being unbounded.

A1.2: Comparing the Simulation and By-hand Calculation of total energy:

In this appendix, I want to verify that the total energy of the electron was conserved and correctly calculated by the simulation. This is important because the system is conservative and must yield a constant total energy. The precision of the energy calculation is also essential because it was used to decide whether the electron is bounded or not. Since we have not looked at the unbounded case where E>0, I will pick a set of initial conditions that yield E>0, then pick two coordinates in their trajectory to calculate the total energy by hand. I can then compare if the energy is conserved and if it agrees with the simulation.

Picked set: (x,y,vx,vy) = [1.1, 4.1, 0, 15]

Energy and coordinate from simulation:

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Constants
k = 8.99e9  # N⋅m^2/C^2 
q = 1.6e-19  # C 
m_e = 9.11e-31  # kg (mass of the electron)
bohr_r = 5.3e-3 #m (bohr's radius)



# Nucleus 1 & 2 Position:
x_atom1 = -5  # meter 
y_atom1 = 0  # meter
x_atom2 = 5 #meter
y_atom2 = 0

#Flag orbit is stable at beginning:
unstable = False

# Differential equations
def diff_eqns(t, state):
    x, y, vx, vy = state
    r1 = np.sqrt((x - x_atom1)**2 + (y - y_atom1)**2)
    r2 = np.sqrt((x - x_atom2)**2 + (y - y_atom2)**2)
    ax1 = -k * q**2 * (x - x_atom1) / (m_e * r1**3)
    ay1 = -k * q**2 * (y - y_atom1) / (m_e * r1**3)
    ax2 = -k * q**2 * (x - x_atom2) / (m_e * r2**3)
    ay2 = -k * q**2 * (y - y_atom2) / (m_e * r2**3)
    ax = ax1 + ax2
    ay = ay1 + ay2
    return [vx, vy, ax, ay]

def event1(t, state):
    x, y = state[0], state[1]
    r1 = np.sqrt((x - x_atom1)**2 + (y - y_atom1)**2)
    return r1 - bohr_r  # Event triggers when r2 equals bohr_r
event1.terminal = True
event1.direction = -1   #Only consider distance decreases toward bohr radius

def event2(t, state):
    x, y = state[0], state[1]
    r2 = np.sqrt((x - x_atom2)**2 + (y - y_atom2)**2)
    return r2 - bohr_r  # Event triggers when r2 equals bohr_r
event2.terminal = True
event2.direction = -1   #Only consider distance decreases toward bohr radius

# Tested Initial conditions:
state0 = [1.1, 4.1, 0, 15]
print(f"For (x,y,vx,vy) = {state0}")
#Time settings
t0, tmax = 0, 10  # seconds
dt = 0.1
t_span = (t0, tmax)

# Solve the system
sol = solve_ivp(diff_eqns, t_span, state0, events=(event1, event2), rtol=1e-6, atol=1e-8, max_step=dt)

#Detection of crashing events:
if len(sol.t_events[0]) > 0:   #Check if event 1 happens (t_events array is not empty)
    print(f"Electron approached nucleus 1 at t = {sol.t_events[0][0]:.6f} seconds.")
    print(f"Position: x = {sol.y_events[0][0][0]:.6e}m, y = {sol.y_events[0][0][1]:.6e}m")  #sol.y_events[0] = [[vx,vy,x,y]] of events 1
    unstable = True
elif len(sol.t_events[1]) > 0:   #Check if event 2 happens (t_events array is not empty)
    print(f"Electron approached nucleus 2 at t = {sol.t_events[1][0]:.6f} seconds.")
    print(f"Position: x = {sol.y_events[1][0][0]:.6e}m, y = {sol.y_events[1][0][1]:.6e}m")
    unstable = True
else:
    print("Electron did not approach within the Bohr radius of either nucleus.")
    
#Mechanical Energy:
x, y, vx, vy = sol.y
r1 = np.sqrt((x - x_atom1)**2 + (y - y_atom1)**2)
r2 = np.sqrt((x - x_atom2)**2 + (y - y_atom2)**2)
PE = -k * q**2 / r1 + -k * q**2 / r2
KE = 0.5 * m_e * (vx**2 + vy**2)
E = PE + KE

#Detection of escaping trajectory:
if (E >= 0).any():
    print("Electron escaped the orbit (E>0)")
    unstable = True
else:
    print("Electron is bounded (E<= 0)")

print(f"Energy of the system: E={E[0]:5g}J")
print(f"Picked position and velocity for by-hand calculation: (x1,y1,vx1,vy1)=[1.1m,4.1m,0m/s,15m/s]; (x2,y2,vx2,vy2) ={x[-1]:5g}m,{y[-1]:5g}m,{vx[-1]:5g}m/s,{vy[-1]:5g}m/s")


# Plot the trajectory
plt.figure(figsize=(8, 8))
plt.plot(x_atom1, y_atom1, 'ro', label='Proton 1')  # Nucleus 1 position
plt.plot(x_atom2, y_atom2, 'ro', label='Proton 2')  # Nucleus 2 position
plt.plot(state0[0], state0[1], color='blue', label='Electron initial position')  # Electron initial position
plt.plot(sol.y[0], sol.y[1], 'b-',  marker='o', markersize = 3, markevery=10, label='Electron Trajectory')  # Electron trajectory
plt.scatter(sol.y[0][0], sol.y[1][0], color='orange', label='Picked Point', zorder=5)  # Picked first point
plt.scatter(sol.y[0][-1], sol.y[1][-1], color='orange', zorder=5)  # Picked last point
plt.xlabel('x (m)')
plt.ylabel('y (m)')
plt.legend()
plt.grid(True)
plt.axis('equal')
plt.title("Electron Trajectory Around Two Protons")
plt.show()

# Plot energy
plt.figure(figsize=(8, 6))
plt.plot(sol.t, PE, label="Potential Energy (PE)")
plt.plot(sol.t, KE, label="Kinetic Energy (KE)")
plt.plot(sol.t, PE + KE, label="Total Energy", linestyle="--")
plt.xlabel("Time (s)")
plt.ylabel("Energy (J)")
plt.legend()
plt.grid(True)
plt.title("Energy of the Electron")
plt.show()

For (x,y,vx,vy) = [1.1, 4.1, 0, 15]
Electron did not approach within the Bohr radius of either nucleus.
Electron escaped the orbit (E>0)
Energy of the system: E=3.05033e-29J
Picked position and velocity for by-hand calculation: (x1,y1,vx1,vy1)=[1.1m,4.1m,0m/s,15m/s]; (x2,y2,vx2,vy2) =0.59532m,100.767m,-0.0733501m/s,8.77365m/s

By-hand Calculation:

From the energy calculation of two different positions in trajectory, the total energy of the system is bigger than 0, which agrees to the behaviour of the simulation. While E2 agrees with the simulation which is roughly 3e-29J, there is a small discrepancy for E1 value which makes it bigger by 7e-27J. Nevertheless, this small discrepancy does not significantly affect the simulation since the orbit is still unbounded.

Appendix 2: Reflection questions

Reflection 1: Coding Approaches (A)

(How well did you apply and extend your coding knowledge in this project? Consider steps you took to make the code more efficient, readable and/or concise. Discuss any new-to-you coding techniques, functions or python packages that you learned how to use. Reflect on any unforeseen coding challenges you faced in completing this project.)

In my opinion, I effectively applied coding techniques taught in class, as well as new techniques I learned from online sources and generative AI. I have applied Joss's tips for solve_ivp, in which I use max_step instead of t_eval for better flexibility in the time step's size of the function. Other typical in-class techniques I have applied include event detection function, adjusting precision with rtol and atol, 3D plot. The new technique I learned when doing this project is using some print statements and diving the code into smaller chunks for debugging, which is my biggest challenge in this project. My simulation got into an issue where the colour of all points in the scatter plot was red (unstable). I then learned from generative AI the print statement that detects any iteration with energy smaller than 0, of which there was none. Meanwhile, I used another code to plot smaller chunks of iterations and calculate energy values, then finally found out the values were inconsistent, meaning energy was incorrectly calculated. Nevertheless, the issue came from the fact that I used iterator k to vary the initial speed, but it was also unintentionally presented in the energy calculation (k constant.

Reflection 2: Coding Approaches (B)

(Highlight an aspect of your code that you feel you did particularily well. Discuss an aspect of your code that would benefit the most from further effort.)

The part of the code I feel like I did well is finding a way to flag and plot the colour for each stable/unstable orbit in 3D space of initial conditions. An aspect that might be most improved if I put more effort into is that finding a way to incorporate larger range of initial conditions without making the code run forever. Since I used nested for loops, it made the simulation run exponentially longer when I increased the range of each initial condition (x,y,vy).

Reflection 3: Simulation phyiscs and investigation (A)

(How well did you apply and extend your physical modelling and scientific investigation skills in this project? Consider the phase space you chose to explore and how throroughly you explored it. Consider how you translated physics into code and if appropriate any new physics you learned or developed a more thorough understanding of.)

To be honest, I think this is a very simple model of physics since it is based mostly on Coulom's Law and some vector manipulation for the direction of acceleration. One new physics I learned when looking for an atom with a larger Bohr's radius for my simulation is the Rydberg atom. Though, I think I effectively simulated and modelled the physics in this project.

Reflection 4: Simulation phyiscs and investigation (B)

(Highlight something you feel you did particularily well in terms of the context of your simulation, the physical modelling that you did or the investigation you performed. Discuss an aspect of these dimensions of your project that would benefit the most from further effort.)

The aspect I feel I did well is to define stable and unstable orbit where I used electron's proximity and bounded energy. Though, I could come up with more definition for unstable orbit like the shape of orbit such as circular or eliptical.

Reflection 5: Effectiveness of your communication

(Highlight something you feel you did particularily well in your visualizations or written communication. Discuss an aspect of your visualizations or written communication that would benefit the most from further effort.)

For this project, I learned from my project 1 to add very specific comments and notes in each important step of code so that it is less confusing for me and TA to revise it. Again, it would be more effective if I learn to visualize the 3D scatter plot in less confusing way, since data points are dense and overlaping. One solution I can think of now is to use less dense data points and use different pair of color for each vy (vertical axe) values.