title: Kinetic Theory Model
body: <p>Edit: I've now fixed the problem I asked about. The spheres were leaving the box in the corners, where the if statements (in the while loop shown below) got confused. In the bits of code that reverse the individual components of velocity on contact with walls, some elif statements were used. When elif is used (as far as I can tell) if the sphere exceeds more than one position limit at a time, the program only reverses the velocity component for one of them. This is rectified when replacing elif with simply if. I'm not sure if I quite understand the reason behind this, so hopefully someone cleverer than I will comment such information, but for now, if anyone has the same problem, I hope my limited input helps!</p>
  
  <p>Some context first:</p>
  
  <p>I'm trying to build a model of the kinetic theory of gases in VPython, as a revision exercise for my (Physics) degree. This involves me building a hollow box and putting a bunch of spheres in it, randomly positioned throughout the box. I then need to assign each of the spheres its own random velocity and then use a loop to adjust the position of each sphere with reference to its velocity vector and a time step. </p>
  
  <p>The spheres should also undergo elastic collisions with each wall and all other spheres.</p>
  
  <p>When a sphere meets a wall in the x-direction, its x-velocity component is reversed and similarly in the y and z directions.
  When a sphere meets another sphere, they swap velocities.</p>
  
  <p>Currently, my code works so far as creating the right number of spheres and distributing them randomly and giving each sphere its own random velocity. The spheres also move as they should, except for collisions. The spheres should all stay inside the box as they should bounce off all the walls. They appear to be bouncing off each other, however, occasionally a sphere or two will go straight through the box.</p>
  
  <p>I am extremely new to programming and I don't quite understand what's going on here or why it's happening but I'd be very grateful if someone could help me.</p>
  
  <p>Below is the code I have so far (I've tried to comment what I'm doing at each step):</p>
  
  <pre><code>##########################################################
  # This code is meant to create an empty box and then create
  # a certain number of spheres (num_spheres) that will sit inside
  # the box. Each sphere will then be assigned a random velocity vector.
  # A loop will then adjust the position of each sphere to make them
  # move. The spheres will undergo elastic collisions with the box walls
  # and also with the other spheres in the box.
  
  ##########################################################
  
  from visual import *
  import random as random
  import numpy as np
  
  num_spheres = 15 
  
  fps = 24 #fps of while loop (later)
  
  dt = 1.0/fps #time step
  
  l = 40 #length of box
  
  w = 2 #width of box
  
  radius = 0.5 #radius of spheres
  
  ##########################################################
  # Creating an empty box with sides length/height l, width w
  
  wallR = box(pos = (l/2.0,0,0), size=(w,l,l), color=color.white, opacity=0.25)
  wallL = box(pos = (-l/2.0,0,0), size=(w,l,l), color=color.white, opacity=0.25)
  wallU = box(pos = (0,l/2.0,0), size=(l,w,l), color=color.white, opacity=0.25)
  wallD = box(pos = (0,-l/2.0,0), size=(l,w,l), color=color.white, opacity=0.25)
  wallF = box(pos = (0,0,l/2.0), size=(l,l,w), color=color.white, opacity=0.25)
  wallB = box(pos = (0,0,-l/2.0), size=(l,l,w), color=color.white, opacity=0.25)
  
  #defining a function that creates a list of 'num_spheres' randomly positioned spheres
  def create_spheres(num):
      global l, radius
      particles = []  # Create an empty list
      for i in range(0,num): # Loop i from 0 to num-1
          v = np.random.rand(3)
          particles.append(sphere(pos= (3.0/4.0*l) * (v - 0.5), #pos such that spheres are inside box
                              radius = radius, color=color.red, index=i))
          # each sphere is given an index for ease of referral later
      return particles
  
  #defining a global variable = the array of velocities for the spheres
  velarray = []
  
  #defining a function that gives each sphere a random velocity
  def velocity_spheres(sphere_list):
       global velarray
       for sphere in spheres:
          #making the sign of each velocity component random
          rand = random.randint(0,1) 
          if rand == 1:
              sign = 1
          else:
              sign = -1
  
          mu = 10 #defining an average for normal distribution
          sigma = 0.1 #defining standard deviation of normal distribution
  
          # 3 random numbers form the velocity vector
          vel = vector(sign*random.normalvariate(mu, sigma),sign*random.normalvariate(mu, sigma),
                   sign*random.normalvariate(mu, sigma))
  
          velarray.append(vel)
  
  
  spheres = create_spheres(num_spheres) #creating some spheres
  velocity_spheres(spheres) # invoking the velocity function
  
  while True:
      rate(fps)
  
      for sphere in spheres:
          sphere.pos += velarray[sphere.index]*dt
          #incrementing sphere position by reference to its own velocity vector
  
          if abs(sphere.pos.x) &gt; (l/2.0)-w-radius:
              (velarray[sphere.index])[0] = -(velarray[sphere.index])[0]
              #reversing x-velocity on contact with a side wall
  
          elif abs(sphere.pos.y) &gt; (l/2.0)-w-radius:
              (velarray[sphere.index])[1] = -(velarray[sphere.index])[1]
              #reversing y-velocity on contact with a side wall
  
          elif abs(sphere.pos.z) &gt; (l/2.0)-w-radius:
              (velarray[sphere.index])[2] = -(velarray[sphere.index])[2]
              #reversing z-velocity on contact with a side wall
  
          for sphere2 in spheres: #checking other spheres
              if sphere2 != sphere:
                  #making sure we aren't checking the sphere against itself
  
                  if abs(sphere2.pos-sphere.pos) &lt; (sphere.radius+sphere2.radius):
                      #if the other spheres are touching the sphere we are looking at
  
                      v1 = velarray[sphere.index]
                      #noting the velocity of the first sphere before the collision
  
                      velarray[sphere.index] = velarray[sphere2.index]
                      #giving the first sphere the velocity of the second before the collision
                      velarray[sphere2.index] = v1
                      #giving the second sphere the velocity of the first before the collision
  </code></pre>
  
  <p>Thanks again for any help!</p>