There's clearly a problem with energy conservation (numerics or other bug), but otherwise it looks very interesting already
https://www.beautifulmathuncensored.de/static/GlowScript/GlowScript.html
from vpython import *
#scene.fullscreen = True
G = 100 #prefactor electric force, code adapted from gravity sim therefore 'G'
B=vector(0,1,10) #constant overall manetic field, should be enough to introduce chaos to 2 particle system (later I'd like to have it caused by the nuclear spin and affect the electron spin, all classical of course)
# cp. https://arxiv.org/abs/chao-dyn/9601009
# Lorentz force is q*cross(v,B), here: charge*cross(velocity,B)
# edit initial conditions here
##########
spheres = [
sphere(pos=vector(2,0,0),radius =.1,color=color.blue,charge=-1,mass=1,velocity=vector(0,-5,0),a = vector(0,0,0)),
#sphere(pos=vector(0,0,0),radius=.3,color=color.red,charge=1,mass=7200,velocity=vector(0,0,0),a=vector(0,0,0)),
#sphere(pos=vector(0,12,0),radius=.08,color=color.green,mass=sqrt(4),velocity=vector(1.2,0,0.6),a=vector(0,0,0),trail=curve(color=color.green)),
sphere(pos=vector(0,0,0),radius=.3,color=color.red,charge=1,mass=7200,velocity=vector(0,0,0),a=vector(0,0,0))
#sphere(pos=vector(0,28,0),radius=.4,color=color.orange,mass=sqrt(80),velocity=vector(0.7,0,0.4),a=vector(0,0,0),trail=curve(color=color.orange)),
#sphere(pos=vector(0,32,0),radius=0.2,color=color.white,mass=-sqrt(10),velocity=vector(1.5,0,0.4),a=vector(0,0,0),trail=curve(color=color.white))
]
def acceleration1on2(sphere2,sphere1):
r = sphere2.pos - sphere1.pos
r_mag = mag(r)
normal_r = norm(r)
g = ((G*sphere1.charge*sphere2.charge)/pow(r_mag,2))/sphere2.mass*normal_r
return g
t = 0
dt = .01
while 1:
rate(60)
for i in spheres:
i.a = vector(0,0,0)
soi = vector(0,0,0)
for j in spheres:
if i!=j:
i.a = i.a + acceleration1on2(i,j)
# only for electron Lorentz force, assume nucleus resting in field (should later be caused by it)
# spheres[0].charge*cross(spheres[0].velocity,B)
lorentz_acc = spheres[0].charge*cross(spheres[0].velocity,B)/spheres[0].mass #acceleration due to B on e-
spheres[0].a += lorentz_acc
for i in spheres:
i.velocity = i.velocity + i.a *dt
i.pos = i.pos+i.velocity*dt
#scene.center=vector(spheres[1].pos.x,spheres[1].pos.y,spheres[1].pos.z)
