title: Implementing the electric and magnetic field into a 2d simulation of free particles
body: <p>I'm looking to introduce electrodynamics into my former ideal gas simulation, which is simply a bunch of particles inside a rectangle. </p>
  
  <p>Edit it appears I was using the wrong equations, and that was the cause of my inital problem.</p>
  
  <p>I'm  using the python 3.5 IDLE environment, <code>numpy</code> for the vector math and <code>pygame</code> for the visualization. The equations I'm using are:
   <a href="http://www.wolframalpha.com/input/?i=moving+point+charge+electric+field" rel="nofollow noreferrer">electric field</a> and <a href="http://www.wolframalpha.com/input/?i=moving+point+charge+magnetic+field" rel="nofollow noreferrer">magnetic field</a></p>
  
  <p>So far so good, my code is running. But to realisticly recreate the physics, the fields need to expand with velocity <code>c</code>. Additionally, when computing the force, I´d have to use the data from the retarded time, but it gets rewritten every new frame. Can I store that data or is it to memory heavy?
   What method do you suggest ?</p>
  
  <p>The code uses the Ball class whose instances are the particles. It has a method associated to it which computes <code>E(r)</code> and <code>B(r)</code>. To compute the force acting on the particle <code>i</code>, I add all of them together at <code>i</code>´s  position and compute the force acting on <code>i</code>.</p>
  
  <p>Here's the code, where 'mass' serves as radius of the particle aswell as charge and the z dimension is just added for the cross product.  The relevant code is the one affected by the <code>if electrodynamics == true</code> statement:</p>
  
  <pre><code>import pygame
  import time
  import random as r
  import numpy as np
  
  pygame.init()
  
  display_width = 800
  display_height = 600
  
  electrodynamics = True
  show_field = False
  
  Ball_num = 7
  c = 10
  red = (255,0,0)
  white = (255,255,255)
  gameDisplay = pygame.display.set_mode((display_width,display_height))
  pygame.display.set_caption('bouncy')
  clock = pygame.time.Clock()
  
  def normalized(a):
      b = np.linalg.norm(a)
      if not(b == 0):
          return a/b
      else:
          return np.zeros((1,3), dtype = float)
  
  def tot_EM_field_at_charge(charges, charge):
  
      EM = np.array([[0.,0.,0.],[0.,0.,0.]], dtype = float)
  
      for q in charges:
          EM = EM + q.EM_field(charge.position)
  
      return EM
  
  def Force_on_bally(field, charge):
      c = np.cross(charge.velocity, field[1])
      force = charge.mass*(field[0] + np.cross(charge.velocity, field[1]))
      return force
  
  class arrow(object):
  
      def __init__(self, length, x, y, charges):
  
  
          self.length = length
          self.position = np.array([x,y,0])
  
          self.field = tot_EM_field_at_charge(charges, self)
          self.field_mag= np.linalg.norm(self.field[0])
          if self.field_mag == 0:            
              self.field_direction = np.zeros(3,)
          else:    
              self.field_direction = self.field[0]/self.field_mag
  
          self.position_end()
          self.color()
          self.show()
      def position_end(self):
  
          self.position_2 = self.position + self.field_direction * self.length
  
          return self.position_2
  
      def color(self):
          self.color = self.field_mag
          if self.color &lt; 0.05:
              self.color = (46,120,255)
          elif self.color &lt; 0.1:
              self.color = (147,145,252)
          elif self.color &lt; 0.3:
              self.color = (249,23,28)
          elif self.color &lt; 0.6:
              self.color =(251,139,33) 
          elif self.color &lt; 1:
              self.color = (255,255,127)
          else:
              self.color = (255,255,255)
  
  
          return self.color
      def show(self):
          pygame.draw.line(gameDisplay, self.color,(int(self.position[0]), int(self.position[1])), (int(self.position_2[0]), int(self.position_2[1])))
  
  
  
  
  
  
  
  
  class Ball(object):
      def __init__(self, x, y, m, c, v_x, v_y):
          self.position = np.array([x,y,0], dtype = float)
          self.position_2 = np.array([x,y,0], dtype = float)
          self.velocity = np.array([v_x, v_y,0], dtype = float)
          self.velocity_2 = np.array([v_x, v_y,0], dtype = float)
          self.acceleration = np.array([0.,0.,0.], dtype = float)
          self.mass = m
          self.color = c
  
      def acceleration_compute(self,force):
          a = force/self.mass
          self.acceleration += a
  
      def move(self):
          self.velocity += self.acceleration
          self.position += self.velocity
          self.acceleration *= 0
  
      def show(self):
          pygame.draw.circle(gameDisplay, self.color, [int(self.position[0]), int(self.position[1])], self.mass)
  
      def Edgelord(self):
          if ((self.position[0] + self.velocity[0] &gt;= display_width-self.mass) and self.velocity[0] &gt; 0):
              self.velocity[0] *= -1
              self.position[0] = display_width - self.mass + self.velocity[0]
  
  
          elif ((self.position[0] + self.velocity[0] - self.mass  &lt;= 0) and self.velocity[0] &lt; 0 ):
  
              self.velocity[0] *= -1
              self.position[0] = self.mass + self.velocity[0] 
  
  
          elif ((self.position[1] + self.velocity[1] &gt;= display_height - self.mass) and self.velocity[1] &gt; 0):
  
              self.velocity[1] *= -1
              self.position[1] = display_height - self.mass + self.velocity[1]
  
  
          elif ((self.position[1] + self.velocity[1] - self.mass  &lt;= 0) and self.velocity[1] &lt; 0 ):
  
              self.position[1] = self.mass -self.velocity[1]
              self.velocity[1] *= -1
  
      def EM_field(self, R):
          radius = np.linalg.norm(R - self.position)
          if radius != 0:
              unitradius = (R - self.position)/radius
          else:
              unitradius = np.zeros(3, )
  
          if np.linalg.norm(radius) != 0 and np.dot(unitradius, self.velocity)!=1:
              charge      = self.mass / (1 - np.dot(unitradius, self.velocity) ** 3)
  
  
              if radius &lt; self.mass:
                  radius = self.mass
  
              radius2     = radius ** 2
  
              velocity_in_c = self.velocity/c
  
              oneMinusV2  = 1 - np.dot(velocity_in_c, velocity_in_c)
              uMinusV     = unitradius - velocity_in_c
              aCrossUmV   = np.cross(uMinusV, self.acceleration);
              Eleft       = (oneMinusV2 * (unitradius - velocity_in_c)) / radius2
              Eright      = np.cross(unitradius, aCrossUmV) / (radius*c**2)
              E           = charge * (Eleft - Eright)
              #E = np.zeros(3, )
              B           = np.cross(unitradius/c, ((charge*c**2) * (Eleft - Eright)))
  
              EM_field = np.array([E,B], dtype = float)
          else:
              EM_field = np.zeros((2,3), dtype = float)
  
          return EM_field
  
  ballys = []
  
  for i in range(Ball_num):
      #ballys.insert(i, Ball(r.randrange(300,display_width - 5, 10),r.randrange(200,display_height/2,1)   , r.randrange(5,10,1),(r.randint(1,255),r.randint(1,255),r.randint(1,255)), r.randint(-200,200)/1000, r.randint(-200,200)/1000))
      ballys.insert(i, Ball(200 + i*50, 220 + i*20 , 10,(r.randint(1,255),r.randint(1,255),r.randint(1,255)),0, 0 ))
  #ballys.append( Ball(300 + 50, 300, 10,(r.randint(1,255),r.randint(1,255),r.randint(1,255)),10, 0 )) 
  
  up = np.zeros(3,)
  down = np.zeros(3,)
  right = np.array([0.,0.,0.])
  left = np.array([0.,0.,0.])
  grav = np.array([0.,0.1,0.])
  repulsion = np.array([0.,0.,0.])
  
  crashed = False
  
  while not crashed :
  
      for event in pygame.event.get():
          if event.type == pygame.QUIT:
              pygame.quit()
              quit()
  
          if event.type == pygame.KEYDOWN:
              if event.key == pygame.K_LEFT:
                  left = np.array([-0.1, 0.,0.])
              if event.key == pygame.K_RIGHT:
                  right = np.array([0.1,0.,0.])
              if event.key == pygame.K_DOWN:
                  down = np.array([0.,0.1,0.])
              if event.key == pygame.K_UP:
                  up = np.array([0.,-0.1,0.])
  
          if event.type == pygame.KEYUP:
              if event.key == pygame.K_LEFT or event.key == pygame.K_RIGHT or event.key == pygame.K_UP or event.key == pygame.K_DOWN :
                  right = np.array([0.,0.,0.])
                  left = np.array([0.,0.,0.])
                  down = np.zeros(3,)
                  up = np.zeros(3,)
  
      gameDisplay.fill(white)
  
      if show_field == True:
          for i in range(display_width//20):
              for j in range(display_height//20):
                  arry = arrow(8, 10 + i*20, 10 + j*20 , ballys)
  
  
  
      if electrodynamics == True:
          for bally in ballys:
              bally.acceleration_compute(Force_on_bally(tot_EM_field_at_charge(ballys, bally), bally))
  
      for i, bally in enumerate(ballys):
  
          #if electrodynamics == True:
            #  bally.acceleration_compute(Force_on_bally(tot_EM_field_at_charge(ballys, bally), bally))
          bally.Edgelord()
  
          bally.acceleration_compute(up)
          bally.acceleration_compute(down)
          bally.acceleration_compute(right)
          bally.acceleration_compute(left)
  
          #ballys[i].acceleration_compute(grav * ballys[i].mass)
  
          for bally2 in ballys[i+1:]:
  
              #checks collisions
              if  np.linalg.norm(bally.position - bally2.position) &lt;= bally.mass  + bally2.mass  :
  
                  bally.velocity_2 = (bally.mass * bally.velocity + bally2.mass * bally2.velocity + bally2.mass *(bally2.velocity - bally.velocity))/ (bally.mass + bally2.mass)
                  bally2.velocity_2 = (bally.mass * bally2.velocity + bally.mass * bally.velocity + bally.mass *(bally.velocity - bally2.velocity))/ (bally2.mass + bally.mass)
  
                  #prevents balls getting stuck in each other and assignes new velocitys
                  if not(np.linalg.norm(bally.position + bally.velocity_2  - (bally2.position + bally2.velocity_2) ) &lt;= bally.mass  + bally2.mass):
                          bally.velocity = bally.velocity_2
                          bally2.velocity = bally2.velocity_2
  
          bally.Edgelord()
  
          bally.move()
          bally.show()
  
      pygame.display.update()
      clock.tick(60)
  
  pygame.quit()
  quit()
  </code></pre>