1. #### Ion collision program

Code to simulate and draw motion of ions. The main function is setup to collide 3 identical ions having equal speed.
The velocities are designed to stop one of the particles. The purpose was to answer a question on quora. Posted here might help with recent forum questions about tkinter
The program also shows how to use tkinter without the mainloop.

Code:
```#import Numbers  # not available with my installation
import sys
import math
import pprint

class Vector(list):

'''
python3 -m doctest -v thisfile.py

doctest  command python3 -m doctest thisfile.py
>>> a = Vector(3,4)
>>> a.length == 5
True
>>> b = Vector(7, 11)
>>> Vector(10, 15) == (a + b)
True
>>> (a * -3) == Vector(-9, -12)
True
>>> (-1 * a) == Vector(-3, -4)
True
>>> a / 2 == Vector(1.5, 2)
True
>>> a.norm.length
1.0
>>> -a
[-3, -4]
'''

def __init__(self, *args):
super().__init__(args)

def dot(a, b):
'dot product'
assert len(a) == len(b)
return sum(x*y for (x,y,) in zip(a, b))

def __imul__(self, other):
raise NotImplementedError('__imul__')

@property
def norm(self):
'Euclidean norm'
return self / self.length

@property
def length(self):
return math.sqrt(self.dot(self))

return Vector(*[x+y for (x,y,) in zip(a,b,)])

def __sub__(a,b):
return Vector(*[x-y for (x,y,) in zip(a,b,)])

def __mul__(self, c):
'multiplication by scalar'
return Vector(*[c*x for x in self])

__rmul__ = __mul__

def __truediv__(self, c):
'division by scalar'
return self * (1.0 / c)

def __neg__(self):
return self * (-1)

class Ion:

'''
>>> a = Ion(position = Vector(0),
...         velocity = Vector(1),
...         charge = 3,
...         mass = 5)
...
>>> a.T
2.5
>>> a.p == Vector(5)
True
>>> a.update(8)  # no force
>>> a.T
2.5
>>> a.update(1, Vector(10))
>>> print(a.velocity, a.position)
[3.0] [11.0]
>>>
>>>
>>> # demonstrate that the ODE solver is good enough for this example
>>>
>>> # expect acceleration of 10 m/s**2 starting at 10 m
>>> # to hit ground at sqrt(2) seconds
>>> a = Ion(Vector(10),Vector(0))
>>> dt = 0.0001
>>> t = 0
>>> while t < math.sqrt(2):
...     a.update(dt, Vector(-10))
...     t += dt
...
>>> print(round(t, 4), list(map(round, a.position, (4,))))
1.4143 [-0.0019]
>>> # kinetic energy should equal potential energy of mgh
>>> print(round(a.T,1), round(1*10*10.0, 1))
100.0 100.0
'''

#    def __init__(self, position:Vector, velocity:Vector, charge = numbers.Number):
def __init__(self, position:Vector, velocity:Vector, charge = 1, mass = 1):
assert len(position) == len(velocity)
self.position = position
self.velocity = velocity
self.charge = charge
self.mass = mass

@property
def speed(self):
return self.velocity.length

def update(self, time_step, force: Vector = None):
self.velocity = self.velocity + self.acceleration(force) * time_step
self.position = self.position + self.velocity * time_step

def acceleration(self, force: Vector):
if not force:
force = self.zero
else:
assert len(force) == len(self.position)
return force/self.mass

def to(self, other):
return other.position - self.position

def momentum(self):
return self.velocity * self.mass

def kinetic_energy(self):
return self.speed ** 2 * self.mass / 2

T = property(fget = kinetic_energy)

p = property(fget = momentum)

@property
def zero(self):
return Vector(*([0]*len(self.position)))

def __str__(self):
return '{}'.format(self.position)

__repr__ = __str__

class System(list):

'''
system of charged particles
'''

def __init__(self, *args):
super().__init__(args)

def update(self, time_step):
zero = self[0].zero
forces = {}
for (i, atom,) in enumerate(self):
to = atom.to
q = atom.charge
for j in range(i+1, len(self)):
b = self[j]
atob = to(b)
magnitude = (q * b.charge) / (atob.length ** 2)
forces[(j,i)] = f = magnitude * atob.norm
forces[(i,j)] = - f
force = sum((forces[key] for key in forces if key[0] == i), zero)
atom.update(time_step, force)

def potential_energy(self):
zero = self[0].zero
potential = 0
for (i, atom,) in enumerate(self):
to = atom.to
q = atom.charge
for j in range(i+1, len(self)):
b = self[j]
atob = to(b)
potential += (q * b.charge) / atob.length
return potential

def momentum(self):
return sum((a.p for a in self), self[0].zero)

def kinetic_energy(self):
return sum(a.T for a in self)

T = property(fget = kinetic_energy)

V = property(fget = potential_energy)

p = property(fget = momentum)

@property
def energy(self):
return self.T + self.V

def __str__(self):
return pprint.pformat(self)

for (i, atom,) in enumerate(self):
to = atom.to
if any(radius < to(self[j]).length for j in range(i+1, len(self))):
return False
return True

if __name__ ==  '__main__':

import tkinter

def initializer(y = -1):
while True:
a = Ion(Vector(-1,0),Vector(1,0))
b = Ion(Vector(1,0),Vector(-1,0))
c = Ion(Vector(0,y),Vector(0,1))
yield System(a,b,c)

def test():
system = initializer()
a = next(system)
print(a)
a.update(1)
print(a)
print(next(system))
a = next(system)
dt = 1e-4
time = 0
p0 = a.p
while a.bounded(2):
if (time % 0.1) < (dt * 2):
print(time)
print(a)
a.update(dt)
time += dt
pn = a.p
# momentum conservation
print(p0, pn) # [0, 1] [0.0, 0.9999999999999962]

def bb(a:Vector, w = 500, b = 8):
'scale and create a bounding box'
(cx, cy) =  a * (w * 0.80 / 2)
return (w/2 + cx - b, w/2 - cy - b, w/2 + cx + b, w/2 - cy + b,)

def main():
root = tkinter.Tk()
root.title('3 body simulation in 2D')
w = 600
canvas = tkinter.Canvas(root, width=w, height=w, background='black')
canvas.pack()
system = initializer(-0.5)
a = next(system)
depictions = []
oval = canvas.create_oval
text = canvas.create_text(w/2,36,font=('fixed', 24), fill='yellow', text='yellow speed: {:4.2f}'.format(a[-1].speed))
for (atom, color,) in zip(a, 'blue green yellow'.split()):
depictions.append(oval(*bb(atom.position, w), fill=color))
root.update()
dt = 1e-4
time = 0
print(a.energy, 'system energy')
print('energy, kinetic and potential')
while a.bounded(3):
a.update(dt)
for (d, atom,) in zip(depictions, a):
canvas.coords(d, *bb(atom.position, w))
canvas.itemconfig(text, text='yellow speed: {:4.2f}'.format(a[-1].speed))
root.update()
if (time % (300 * dt)) < dt:
T = a.T
V = a.V
print(('{:8.3f}'*3).format(T+V, T, V))
sys.stdout.flush()
time += dt
root.mainloop()
print(a.energy, 'system energy')

main()```
2. Interesting program, although I still try to figure out how it works.

But I noticed you do have a mainloop call near the end/bottom of script
3. That's just to keep the window open when computation completes. Could just as well have been "press enter to finish" prompt in console.