#1
  1. Contributing User
    Devshed God 1st Plane (5500 - 5999 posts)

    Join Date
    Aug 2011
    Posts
    5,971
    Rep Power
    509

    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 __iadd__(self, other):
            raise NotImplementedError('__iadd__')
    
        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))
    
        def __add__(a,b):
            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)
    
        def bounded(self, radius):
            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()
    [code]Code tags[/code] are essential for python code and Makefiles!
  2. #2
  3. Lord of the Dance
    Devshed Specialist (4000 - 4499 posts)

    Join Date
    Oct 2003
    Posts
    4,162
    Rep Power
    2011
    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
  4. #3
  5. Contributing User
    Devshed God 1st Plane (5500 - 5999 posts)

    Join Date
    Aug 2011
    Posts
    5,971
    Rep Power
    509
    That's just to keep the window open when computation completes. Could just as well have been "press enter to finish" prompt in console.
    [code]Code tags[/code] are essential for python code and Makefiles!

IMN logo majestic logo threadwatch logo seochat tools logo