title: Accelerating the calculations in Python (Simulation of particles in a magnetic field)
body: <p>A have a problem with the speed of the program written in Python. The program is "Simulation of ferromagnetic particles in a magnetic field", and more specifically in the magnetically inert liquid. The program works, but very slowly compared to the same program written in C++, but I have it written in Python, on a project to study.</p>
  
  <p>Overall, the program code is based on loops, with lots of floating point calculations. There is a random number of particles (generated randomly position) which interact with each other under the influence of a magnetic field.</p>
  
  <p>This is the initial position:</p>
  
  <p><a href="http://i.stack.imgur.com/T15Bb.jpg" rel="nofollow">http://i.stack.imgur.com/T15Bb.jpg</a> </p>
  
  <p>And the final:</p>
  
  <p><a href="http://i.stack.imgur.com/0nU5D.jpg" rel="nofollow">http://i.stack.imgur.com/0nU5D.jpg</a></p>
  
  <p>The main loop (in SymMain.py, with k variable) iterations are time steps are calculated coordinates of the particles in them at the time and the forces acting on it (the attractive force and a small repulsive force). In order to speed up, I wanted to use parallel processing to calculate the number of iterations at the same time, but this is not possible, because the calculation in one iteration depend on the calculations in the previous one.</p>
  
  <p>I did not know that Python is so slower than C + +. As an example, the simulation of 529 particles in one time step is calculated (on my computer): </p>
  
  <p>C + + ~ 0.5s</p>
  
  <p>Python ~ 50s</p>
  
  <p>So, how to speed up the program? Please help.</p>
  
  <blockquote>
    <blockquote>
      <blockquote>
        <p>This code is only calculate, not drawing particles as in the pictures. The number of simulate partcile is on SymConst.py, this is nrH*nrL.</p>
      </blockquote>
    </blockquote>
  </blockquote>
  
  <p>SymMain.py</p>
  
  <pre><code>#coding:windows-1250
  
  from os import system
  from SymCalc import *
  from SymParticle import *
  
  
  if __name__ == '__main__':
      App = SymCalc()
      App.MainLoop()
  </code></pre>
  
  <p>SymParticle.py</p>
  
  <pre><code>#coding:windows-1250
  
  from random import randint
  from math import *
  from SymConst import *
  from SymParticle import *
  
  class SymCalc(object):
      def __init__(self):
          # declaration lists containing the properties of the particles
          ParticleList = []
          ParticleListTemp = []
          t = 0.0
  
          # the initial values ​​of particle
          for x in range(0, nParticle):
              ParticleList.append(Particle(x+1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 8e-6))
              ParticleListTemp.append(Particle(x+1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 8e-6))
  
          # generating random coordinates X, Y 
          for x in range(0, nParticle):
              self.Rand(ParticleList[x], x)
  
          # time steps
          for k in range(0, k0):
              print "Time step = {0}".format(k)
  
              # calculation of forces
              for i in range(0, nParticle):
                  for j in range(0, i-1):
                      self.ParticleCalculate(ParticleList[i], ParticleList[j], 1)
                  for j in range(i+1, nParticle):
                      self.ParticleCalculate(ParticleList[i], ParticleList[j], 1)
  
              # display data
              if(k%Constant == 0):
                  for i in range(0, nParticle):
                      self.Print(k, i, ParticleList[i], t)
  
              # changing the position of the particle
              for i in range(0, nParticle):
                  self.ChangeOfPosition (ParticleList[i], dt)
  
                  # reset forces
                  for j in range(0, nParticle):
                      self.ParticleCalculate(ParticleList[i], ParticleList[j], 0)
                      self.ParticleCalculate(ParticleList[i], ParticleListTemp[j], 0)
  
              # next time step
              t += dt
  
  
      # random coordinates of the particles
      def Rand(self, part, lp):
          l = lp%nrL    # współrzędna X pola      
          part.x = (l-nrL/2)*(L0/nrL) + ((randint(0,32767)%100)-50)*(L0/nrL-2*part.r)/99    # X
  
          h = (lp+1-(lp)%nrL)/nrL    # Y      
          part.y = (h-nrH/2)*(H0/nrH) + ((randint(0,32767)%100)-50)*(H0/nrH-2*part.r)/99    # współrzędna Y cząstki
  
  
      # calculating function of the force acting on the particles
      def ParticleCalculate(self, part, part_tmp, p):
          # auxiliary variables
          # r_4, dx, dy, rr, sum, sx, sy, chi_eff, tmp, frepx, frepy, f0
          md = []
  
          if(p == 0):
              part.frx = 0
              part.fry = 0
              part.fx = 0
              part.fy = 0
  
          if(p == 1):
              # versor coordinates connecting the geometrical center of the particle
              dx = part.x - part_tmp.x
              dy = part.y - part_tmp.y
  
              # the distance between the geometric means of the particle
              rr = sqrt(dx*dx + dy*dy)
  
              if(rr &lt; 0.85*(part.r + part_tmp.r)):
                 print "ERROR: Invalid distance between the particles! Simulation aborted..."
  
              if(rr &gt;= 10*part.r):
                  # magnetic dipoles
                  chi_eff = (3.*(MI_P - 1.))/((MI_P - 1) + 3.)
  
                  md.append((4.*MI_0*pi*part.r*part.r*part.r*chi_eff*M_H0)/3.0)
                  md.append((4.*MI_0*pi*part_tmp.r*part_tmp.r*part_tmp.r*chi_eff*M_H0)/3.0)
  
                  tmp = pow(rr,7)
  
                  # first member
                  sum = (5.*(md[0]*part.nx*dx + md[0]*part.ny*dy) * (md[1]*part_tmp.nx*dx + md[1]*part_tmp.ny*dy)) / tmp
                  sx = sum*dx
                  sy = sum*dy
                  tmp = tmp / (rr*rr)
  
                  # second member
                  sum = (md[0]*part.nx*md[1]*part_tmp.nx + md[0]*part.ny*md[1]*part_tmp.ny) / tmp
                  sx -= sum*dx
                  sy -= sum*dy
  
                  # third member
                  sx -= (md[0]*(md[1]*part_tmp.nx*dx + md[1]*part_tmp.ny*dy)*part.nx + md[1]*(md[0]*part.nx*dx + md[0]*part.ny*dy)*part_tmp.nx)/tmp
                  sy -= (md[0]*(md[1]*part_tmp.nx*dx + md[1]*part_tmp.ny*dy)*part.ny + md[1]*(md[0]*part.nx*dx + md[0]*part.ny*dy)*part_tmp.ny)/tmp
  
                  # finally
                  tmp = (-3./(4*pi*MI_0))
                  sx *= tmp
                  sy *= tmp
                  part.fx += sx
                  part.fy += sy
  
                  # short-range repulsive force
                  tmp = pow(rr,15)
                  r_4 = pow((part.r+part.r),4)
  
                  f0 = 3.*fabs(md[0])*fabs(md[1])/(4.*pi*MI_0*r_4)
                  frepx = pow((part.r+part.r),15)*dx/(tmp*rr)*f0
                  frepy = pow((part.r+part.r),15)*dy/(tmp*rr)*f0
  
                  part.frx += frepx;
                  part.fry += frepy;
  
  
      # change the position of the particle
      def ChangeOfPosition(self, part, dt):
          part.ddx = 0
          part.ddy = 0
  
          part.ddx = 1/(6*pi*part.r*eta)*(part.fx+part.frx)*dt
          part.ddy = 1/(6*pi*part.r*eta)*(part.fy+part.fry)*dt
  
          # particles new position value
          part.x += part.ddx
          part.y += part.ddy
  
          if(part.x &lt; -L0/2):
              part.x += L0
          elif(part.x &gt; L0/2):
              part.x -= L0
          elif(part.y &lt; -H0/2):
              part.y += H0
          elif(part.y &gt; H0/2):
              part.y -= H0
  
  
  
      # display data
      def Print(self, k, i, part, t):
          print "-"*50
          print "\nParticle {0}".format(i+1)
          print "The resultant magnetostatic force fx = {0}".format(part.fx - part.frx)
          print "The resultant magnetostatic force fy = {0}\n".format(part.fy - part.fry)
          if(i == nParticle-1):
              print "\n\t\t---t={0}[s]---".format(t)
          print "-"*50
  </code></pre>
  
  <p>SymParticle.py</p>
  
  <pre><code>#coding:windows-1250
  
  class Particle(object):
      # generating a particle properties
      def __init__(self, num, x, y, fx, fy, frx, fry, nx, ny, ddx, ddy, r):
          self.num = num     
          self.x = x     
          self.y = y     
          self.fx = fx     
          self.fy = fy     
          self.frx = frx     
          self.fry = fry     
          self.nx = nx     
          self.ny = ny     
          self.ddx = ddx     
          self.ddy = ddy     
          self.r = r     
  </code></pre>
  
  <p>SymConst.py</p>
  
  <pre><code>#coding:windows-1250
  
  ### Constant
  M_H0 = 3e4
  MI_0 = 12.56e-7     
  MI_P = 2000
  eta = 0.1           
  k0 = 10001
  dt = 1e-6           # time step      
  H0 = 5.95e-4
  L0 = 5.95e-4
  nrH = 4
  nrL = 4
  Constant = 100
  nParticle = nrH*nrL    # number of particle
  </code></pre>