What Could be Causing Discrepancies in My Python Gas Simulation?

[leibo]

Hello.

I wrote a simple gas simulation in python, with pygame. In the simulation the atoms collide elastic collisions with each other and with the walls. Although the simulation does run and the collisions look realistic, when I try to calculate physical parameters like p,T,V,n and compare them to the gas law's there is no match. For example, when the radius is doubled the mean free path becomes smaller, but not in a factor of 2^2=4. the ratio PV/NT is not exactly constant and can be changed by a factor of 10, etc. I wonder what can be the problem? thank you.

this is the code:

from __future__ import division
import pygame, sys,pygame.mixer
import math
import pylab
import numpy as np
import time
import random
from pygame.locals import*
pygame.init()

"""basic variables"""
clock = pygame.time.Clock()
font1 = pygame.font.Font(None, 17)
edgeX,edgeY = 10,10
boxWidth,boxHeight = 950,680
size = width,height = 1000,700
centerX,centerY=int(width/2),int(height/2)
screen = pygame.display.set_mode(size)
x,y = width/2, height/2
radius = 3
number_of_atoms =200
position = []
v = []
color = []

"""set initial positions and velocities"""

for j in range(number_of_atoms):
    position.append([400+random.randint(-250,320),325+random.randint(-285,375),0])
    v.append([random.uniform(-5,5),random.uniform(-5,5)])
    m.append(1)
    color.append((0,0,255))
color[number_of_atoms-1]=((250,0,0))
for i in range(4):
 v.append([0,0])

def colission(a,v,b): # a - list of all atoms, v - list of all velocites (atom j has in position position[j] and travel with velocity v[j]), b - the atom we check for collisions
 for ball in a:
  if ball is not b:
   distanceX,distanceY = (b[0]-ball[0]),(b[1]-ball[1])
   distance = distanceX*distanceX + distanceY*distanceY
   if distance < 36: #radius=3,(2radius)**2=36
    ma,mb=1,1 #the mass of all atoms is 1
    i,j = a.index(ball),a.index(b)
    change = [-0.15*distanceX,-0.15*distanceY]

    Vax,Vay = v[i][0],v[i][1]
    Vbx,Vby = v[j][0],v[j][1]

    Angle =math.atan2(distanceY,distanceX) # math.acos((distanceX/distance))
    cos,sin = math.cos(Angle),math.sin(Angle)

    #changing coordiantes:
    vg,vh=[Vag,Vah],[Vbg,Vbh]=[Vax*cos+Vay*sin,Vay*cos-Vax*sin],[Vbx*cos+Vby*sin,Vby*cos-Vbx*sin]

    #after collision:
    Va2,Vb2=(2*mb*Vbg)/(ma+mb), (2*ma*Vag)/(ma+mb) #because ma=mb the term (ma-mb)*Vag or (mb-ma)*Vbg  cancels out

    v1,v2 = [Va2,Vah] , [Vb2,Vbh]
    #changing coordiantes again:
    sin = -sin
    v1,v2 =  [(Va2*cos+Vah*sin),Vah*cos-Va2*sin], [v2[0]*cos+v2[1]*sin,v2[1]*cos-v2[0]*sin]
    v[j],v[i]= v2,v1 #update velocites
    a[i][0],a[i][1]=a[i][0] + change[0], a[i][1] + change[1] #treating overlaps
    return()

"""the loop"""

while 1:
    screen.fill((0,0,0))

    mouseX,mouseY=pygame.mouse.get_pos()
    for event in pygame.event.get():
     if event.type == pygame.QUIT:
        pygame.quit()
        sys.exit()
     elif event.type == MOUSEBUTTONDOWN:
        if mouseY<edgeY+50 and mouseY>edgeY-50 and mouseX<edgeX+boxWidth and mouseX>edgeX-50:
         z=1
     elif event.type == MOUSEBUTTONUP:
        z=0
    if z==1:
        dy = mouseY-edgeY
        boxHeight = boxHeight-dy
        edgeY=mouseY

    pygame.draw.rect(screen, (200,100,0), (edgeX,edgeY,boxWidth,boxHeight),1)
    for event in pygame.event.get():
     if event.type == pygame.QUIT:
        sys.exit()

    for j in range(len(position[:number_of_atoms])):
        colission(position,v,position[j]) #check for collisions
        position[j][0]+=v[j][0] 
        position[j][1]+=v[j][1]

        """checking for collisions with the walls"""
        if position[j][0]+radius>edgeX+boxWidth:
         
         v[j][0]=-1*(v[j][0])
         position[j]=[(-radius+edgeX+boxWidth)-2,position[j][1],0]

        if position[j][0] - radius<edgeX:
         
         v[j][0]=-1*(v[j][0])
         position[j]=[radius+edgeX+2,position[j][1],0]

        if position[j][1]+radius>edgeY+boxHeight:
        
         v[j][1]=-1*(v[j][1])
         position[j]=[position[j][0],(-radius+edgeY+boxHeight)-2,0]

        if position[j][1]-radius<edgeY:
       
         v[j][1]=-1*(v[j][1])
         position[j]=[position[j][0],edgeY+2+radius,0]

    clock.tick(40)
    for i in range(len(position[:number_of_atoms])):
     circle = pygame.draw.circle(screen,color[i], (int(position[i][0]),int(position[i][1])),radius)
    pygame.display.flip()

 

[uart]

Hello.

I wrote a simple gas simulation in python, with pygame. In the simulation the atoms collide elastic collisions with each other and with the walls. Although the simulation does run and the collisions look realistic, when I try to calculate physical parameters like p,T,V,n and compare them to the gas law's there is no match. For example, when the radius is doubled the mean free path becomes smaller, but not in a factor of 2^2=4. the ratio PV/NT is not exactly constant and can be changed by a factor of 10, etc. I wonder what can be the problem? thank you.

Nice simulation leibo. It looks like it would be good for demonstration purposes.

I think one of the issues that you've got there is the fact that it's 2D instead of 3D like a real (or ideal) gas. I'm not an expert on this but just doing some quick calculations I think that in the 2D model you should still get PA = const * NT, (where here I've defined P as the force per unit length for the 2D case). Note also that area (A) replaces volume (V) in this case.

As for the mean free path, that's another calculation that you're going to have to make adjustments for in the 2D case versus the usual 3D case. Again I could be wrong, but my calculations suggest that the MFP should be proportional to 1/R (instead of 1/R^2) in this 2D case.

I couldn't see anywhere in your code where you were calculating any quantities such as temperature (rms velocity) or pressure. Since you say you're getting inconstant values for these quantities then I assume you must be calculating them somewhere. Was this code removed to make it more compact to post here?

BTW. I had to make a few small additions to your code to get it to run on my computer (maybe different versions of Python or Pygame). Just simple uninitialized variables. I had to add m=[] and z=0 to the initialization section. Not initializing "m" looks like an oversight, but "z" looks likes it's supposed to be initialized to either 0 or 1 somewhere within an if ... elif ... elif ... structure. I didn't fully trace the code to find where, but it must be getting through that block somehow without setting "z", so that also caused the program to crash on my computer. Anyway, putting "z=0" in the initialization section fixed the issue.

 

[gsal]

Here is a gas simulation using visual python.

gas-program

 

[leibo]

Nice simulation leibo. It looks like it would be good for demonstration purposes.

I think one of the issues that you've got there is the fact that it's 2D instead of 3D like a real (or ideal) gas. I'm not an expert on this but just doing some quick calculations I think that in the 2D model you should still get PA = const * NT, (where here I've defined P as the force per unit length for the 2D case). Note also that area (A) replaces volume (V) in this case.

As for the mean free path, that's another calculation that you're going to have to make adjustments for in the 2D case versus the usual 3D case. Again I could be wrong, but my calculations suggest that the MFP should be proportional to 1/R (instead of 1/R^2) in this 2D case.

I couldn't see anywhere in your code where you were calculating any quantities such as temperature (rms velocity) or pressure. Since you say you're getting inconstant values for these quantities then I assume you must be calculating them somewhere. Was this code removed to make it more compact to post here?

BTW. I had to make a few small additions to your code to get it to run on my computer (maybe different versions of Python or Pygame). Just simple uninitialized variables. I had to add m=[] and z=0 to the initialization section. Not initializing "m" looks like an oversight, but "z" looks likes it's supposed to be initialized to either 0 or 1 somewhere within an if ... elif ... elif ... structure. I didn't fully trace the code to find where, but it must be getting through that block somehow without setting "z", so that also caused the program to crash on my computer. Anyway, putting "z=0" in the initialization section fixed the issue.

Thank you very much for your answer. I took into consideration the adjustments for 2D, yet the data I get is not enough accurate. I did not include the measuring code for simplicity. Since the temperature is proportional to the kinetic energy (to be precise - it is proportional to the RMS velocity squared, but the difference is not big), I measured the kinetic energy by adding the following code to the while loop:

EK = 0.5*sum([velocity[0]**2 for velocity in v[:number_of_balls]]) + 0.5*sum([velocity[1]**2 for velocity in v[:number_of_balls]])

(I summed the velocities only up to up to number_of_balls since the four last velocities are the "velocities of the walls "- the purpose of adding that is energy-conservation and momentum-conservation tests I made to check whether the collision function indeed represents elastic collisions. the momentum is fully conserved, but for some reason the kinetic energy is not fully conserved and oscillates a little - maybe this is the problem?)

the area A is boxWidth*boxHeight.

the pressure was measured as follows: with every impact of an atom with the right wall, the atom transfers to the wall momentum of 2mVx. since m=1 the total momentum transferred is dp=2V[j][0]. after every 300 iterations, I divided it by the length of the wall and by 300 (the "time", since F=dp/dt). In the code I added a "pressure" list (initially empty list) and an integer k that counts the number of iterations. I added these lines to the part in the code dealing with collisions with the right wall:

 dp+=2*v[j][0]
if k%300 == 0:
 p = dp/300*boxHeight
 print p
 dp = 0

for measuring the mean free path I created an empty list called "MFP" and added this line to the collision-function:

if (i  == number_of_balls-1 or j  == number_of_balls-1) and (i not in MFP and j not in MFP):
        MFP.append(position[i])

(I checked the mean free path of the atom with index number_of_balls-1).
from this list of coordinates of collisions one can easily find the average MFP:

MFP = np.array(MFP)
     d = []
     for i in range(1,len(MFP)):
      d.append(MFP[i]-MFP[i-1]) #d=[x,y] when x = horizontal distance between 2 collision, y =vertical distance between 2 collision
     distance = []
     for diffrence in d:
      distance.append(math.sqrt((diffrence[0])**2 + (diffrence[1])**2))
     print "MFP=",np.average(distance)
     MFP=MFP.tolist()

I am sorry for the missing z. This z is a part of the section in the code that allows you to scroll down the upper wall of the box with the mouse and reducing the area of the box, during the running of the simulation.

For better results one can take the average of the values of P,T printed during the running of the program but yet the change is not very significant. By the way, if you run the simulation for some time, letting the atoms to spread the kinetic energy between them by collisions, the velocity distribution is pretty much like Maxwell–Boltzmann distribution (as opposed to the initial velocity distribution).

this is the full code:

from __future__ import division
import pygame, sys,pygame.mixer
import math
import pylab
import numpy as np
import time
import random
from pygame.locals import*
pygame.init()

"""basic variables"""
clock = pygame.time.Clock()
font1 = pygame.font.Font(None, 17)
edgeX,edgeY = 0,0
boxWidth,boxHeight = 1000,700
size = width,height = 1000,700
centerX,centerY=int(width/2),int(height/2)
screen = pygame.display.set_mode(size)
x,y = width/2, height/2
radius = 3
number_of_atoms =230
position = []
v = []
color = []
m = []
z = 0
dp = 0
k = 0
MFP = []
"""set initial positions and velocities"""

for j in range(number_of_atoms):
    position.append([edgeX+random.randint(0,boxWidth),edgeY+random.randint(0,boxHeight),0])
    v.append([random.uniform(-5,5),random.uniform(-5,5)])
    m.append(1)
    color.append((0,0,255))
color[number_of_atoms-1]=((250,0,0))
for i in range(4):
 v.append([0,0])

def colission(a,v,b): # a - list of all atoms, v - list of all velocites (atom j is in position position[j] and travel with velocity v[j]), b - the atom we check for collisions
 for ball in a:
  if ball is not b:
   distanceX,distanceY = (b[0]-ball[0]),(b[1]-ball[1])
   distance = distanceX*distanceX + distanceY*distanceY
   if distance < 36: #radius=3,(2radius)**2=36
    ma,mb=1,1 #the mass of all atoms is 1
    i,j = a.index(ball),a.index(b)
    if (i  == number_of_balls-1 or j  == number_of_balls-1) and (i not in MFP and j not in MFP):
        MFP.append(position[i])
    change = [-0.15*distanceX,-0.15*distanceY]

    Vax,Vay = v[i][0],v[i][1]
    Vbx,Vby = v[j][0],v[j][1]

    Angle =math.atan2(distanceY,distanceX) # math.acos((distanceX/distance))
    cos,sin = math.cos(Angle),math.sin(Angle)

    #changing coordiantes:
    vg,vh=[Vag,Vah],[Vbg,Vbh]=[Vax*cos+Vay*sin,Vay*cos-Vax*sin],[Vbx*cos+Vby*sin,Vby*cos-Vbx*sin]

    #after collision:
    Va2,Vb2=(2*mb*Vbg)/(ma+mb), (2*ma*Vag)/(ma+mb) #because ma=mb the term (ma-mb)*Vag or (mb-ma)*Vbg  cancels out

    v1,v2 = [Va2,Vah] , [Vb2,Vbh]
    #changing coordiantes again:
    sin = -sin
    v1,v2 =  [(Va2*cos+Vah*sin),Vah*cos-Va2*sin], [v2[0]*cos+v2[1]*sin,v2[1]*cos-v2[0]*sin]
    v[j],v[i]= v2,v1 #update velocites
    a[i][0],a[i][1]=a[i][0] + change[0], a[i][1] + change[1] #treating overlaps
    return()

"""the loop"""

while 1:
    k+=1
    screen.fill((0,0,0))

    mouseX,mouseY=pygame.mouse.get_pos()
    for event in pygame.event.get():
     if event.type == pygame.QUIT:
        pygame.quit()
        sys.exit()
     elif event.type == MOUSEBUTTONDOWN:
        if mouseY<edgeY+50 and mouseY>edgeY-50 and mouseX<edgeX+boxWidth and mouseX>edgeX-50:
         z=1
     elif event.type == MOUSEBUTTONUP:
        z=0
    if z==1:
        dy = mouseY-edgeY
        boxHeight = boxHeight-dy
        edgeY=mouseY

    pygame.draw.rect(screen, (200,100,0), (edgeX,edgeY,boxWidth,boxHeight),1)
    for event in pygame.event.get():
     if event.type == pygame.QUIT:
        sys.exit()

    #checking for the temperature. Can be any other integer rather than 300
    if k%300 == 0:
      EK = 0.5*sum([velocity[0]**2 for velocity in v[:number_of_balls]]) + 0.5*sum([velocity[1]**2 for velocity in v[:number_of_balls]])
      print "T=",EK

    # checking the pressure:
      p = dp/(300*boxHeight)
      print "P=",p
      dp = 0

    for j in range(len(position[:number_of_atoms])):
        colission(position,v,position[j]) #check for collisions
        position[j][0]+=v[j][0]
        position[j][1]+=v[j][1]

        """checking for collisions with the walls"""

        if position[j][0]+radius>edgeX+boxWidth:
        # calculating pressure
         dp+=2*v[j][0]

         v[j][0]=-1*(v[j][0])
         position[j]=[(-radius+edgeX+boxWidth)-2,position[j][1],0]

        if position[j][0] - radius<edgeX:

         v[j][0]=-1*(v[j][0])
         position[j]=[radius+edgeX+2,position[j][1],0]

        if position[j][1]+radius>edgeY+boxHeight:

         v[j][1]=-1*(v[j][1])
         position[j]=[position[j][0],(-radius+edgeY+boxHeight)-2,0]

        if position[j][1]-radius<edgeY:

         v[j][1]=-1*(v[j][1])
         position[j]=[position[j][0],edgeY+2+radius,0]
    #checking foe mean free path
    if k%1000 == 0:
     MFP = np.array(MFP)
     d = []
     for i in range(1,len(MFP)):
      d.append(MFP[i]-MFP[i-1]) #d=[x,y] when x = horizontal distance between 2 collision, y =vertical distance between 2 collision
     distance = []
     for diffrence in d:
      distance.append(math.sqrt((diffrence[0])**2 + (diffrence[1])**2))
     print "MFP=",np.average(distance)
     MFP=MFP.tolist()

    clock.tick(40)
    for i in range(len(position[:number_of_atoms])):
     circle = pygame.draw.circle(screen,color[i], (int(position[i][0]),int(position[i][1])),radius)
    pygame.display.flip()

Here is a gas simulation using visual python.

gas-program

I know this program (for some reason it make my IDE crush) but I want to know what the problem with my program (-:

please forgive me for my english, I am not a native speaker.

 

[leibo]

This is the most updated version. now the kinetic energy (temperature) is calculated with Vrms. The ratio PA/NT is almost constant and stands on about 1.03 (when you scroll down the upper wall , reducing the area to about half of the original area, the pressure is doubled (-:)

from __future__ import division
import pygame, sys,pygame.mixer
import math
import pylab
import numpy as np
import time
import random
from pygame.locals import*
pygame.init()

"""basic variables"""
clock = pygame.time.Clock()
edgeX,edgeY = 0,0
boxWidth,boxHeight = 1000,700
size = width,height = 1000,700
centerX,centerY=int(width/2),int(height/2)
screen = pygame.display.set_mode(size)
x,y = width/2, height/2
radius = 3
number_of_atoms =230
position = []
v = []
# this is velocity distribution close to Maxwell–Boltzmann distribution
# v = [[-2.1287244853552325, 1.129085317619489], [-1.4575898925432307, -1.2711725536523415], [-2.2097302277479645, 2.7202710712771894], [0.5886788765323574, -3.795609370434213], [6.1789538542622955, -4.033423565067819], [-2.0440228542602634, 1.091295730259362], [1.6132212358775329, -0.5907913294550777], [-0.2516175054065206, -0.7569982911659228], [-0.767834098944953, -1.2465024578108552], [2.9211589754032574, 1.3399207404716227], [1.827119555517413, 0.4861092239445045], [1.4719966556288928, -1.8110668022930578], [0.014573174901874097, -1.6394058668971683], [-0.4434283584239646, 0.623965993550651], [3.0252604143752726, 2.208635209669537], [2.3417487342314587, 3.808728601455674], [-3.471920682852678, -2.87246410407716], [0.5985499889182406, 1.353504417528499], [-2.416855783145403, 1.292420605968518], [-1.443024546146792, -0.3755611375280395], [-5.2312698750465625, 0.3508321653804425], [-3.2885591872235773, 0.2360824989220558], [4.9583454096399056, 3.635330201978962], [-0.17681495748929366, -1.9062349065461373], [5.052487371301375, -0.7498305180096287], [-1.091483827298553, -1.3134606347484017], [-4.569292937508718, 2.2657925196920705], [0.024862400530593964, 1.038245416716403], [-2.7612212253634536, -2.4362409251158947], [3.751202862544148, -1.2884358964475193], [-5.544759276849922, 1.633206866449836], [-1.4929803152561631, 7.427959261771754], [0.6910384380395272, -2.3210682127246605], [4.966179676366409, 8.14277050918934], [-2.939921609018506, 4.247572878924255], [3.46280491334585, -0.22832898919642597], [3.7568113482449323, -3.9818124252045584], [-5.81603207616, -1.8953289928657058], [1.8922119487626696, 5.647835085379633], [1.5865111481144616, 2.440716397288515], [-2.559241988449605, 6.259315242017289], [1.560095547900185, -2.43209068738827], [0.2948421732901573, 0.30410261965283103], [-1.3397082002559542, -3.6253780052346816], [-3.605609964313315, 0.3837913068513812], [0.5357050593902128, 5.3593553783442704], [-1.792412888558205, 0.862244087962059], [1.162155955027183, -2.9130269821889017], [-3.4850237012054675, -7.674228647409126], [4.656126575705954, 4.584905781321144], [-0.7571297335955761, -0.5059533871455572], [-1.2719565404060598, -5.869202492280124], [-1.0217834106552188, -0.7805395470606652], [-0.8630572822069764, 3.0094803400894166], [4.200359299186116, -0.578968615806449], [-0.7869975423781066, 4.19567092164654], [2.633443039109919, 2.2318655698837317], [3.602440137034528, 2.294935016747216], [3.4289832767298516, 0.4544258244837671], [4.513819006019582, -1.2295456367701252], [-0.07632380913320802, 1.2072539498543464], [-1.8068756201328224, 2.74995848527786], [-2.976227145535782, -6.563635270523498], [0.5026244572731093, 3.4693018984139234], [2.996202980131871, 4.699392928649715], [-3.2899406516732057, 1.3134134276847456], [1.718079815818962, -0.1683187895044323], [-0.13518393017429897, 5.54700587906911], [1.8790135713840248, 2.0266924690419112], [3.5647748704942472, -0.02369178903435687], [-2.3926060753474827, 3.657019905033295], [6.4969728497115256, -1.0367301642558044], [-0.023465299643890325, 2.3272129564146233], [-0.6673448250322047, 3.716099643539505], [-4.188916072667548, 1.8798793523488007], [1.6985531474315374, -4.256524967353668], [-3.6760110444317227, 0.450815452016781], [0.8593856519469103, 6.548276513256727], [-0.6874305714315574, 3.924788489008713], [-1.6280123105285504, 4.3272921365645685], [1.9172239398967272, 2.256850677753119], [0.15416495499730978, -2.431193698081671], [1.5053797733536585, -0.23695421632733982], [-2.9137062125131745, -1.8971562459288949], [2.8541364655469073, -0.456823611080871], [-0.09661346709533003, 2.9051650821793364], [0.5072015945158301, -2.2672681386822875], [0.2545044950869553, 4.621789703077026], [-2.243872975995369, -0.8370408949452286], [-0.8941426424681549, -7.172708681656136], [0.20169032842291273, 2.4814670906899496], [-1.8964009017655787, -3.8470469473042566], [0.1977415306381295, 0.9954546081169863], [1.2177433260585913, -3.516139330055328], [1.4584436116027801, -1.5102918861304389], [-4.688466851545976, -0.18279688006604267], [1.1066366006626174, 1.1094869103094414], [-5.898434347641582, -1.2255529598860404], [-1.1645661115261157, 1.0758904080219807], [0.19564921451363354, 0.942295215044142], [-0.958888900885171, 1.858851545058297], [0.86195445606525, 0.0560058442935831], [-1.8794681665160713, -0.5882204653385881], [-2.441063036455345, 0.732607361354868], [2.7193057366101, 1.9389437680157378], [-1.5313904527595503, -1.5367326619406922], [-4.128124335095066, 1.0843340884023764], [-1.7533476670721027, -2.208124824940983], [-1.6950872812235964, -2.277644406649925], [4.063102097910239, -0.6585604177104856], [1.3842578839330968, -2.9221362400820685], [0.47579102203538015, -0.8945850746569246], [0.8453488171596295, 2.052097091864293], [-2.6443681964978722, -0.2057732479640893], [1.5566398411524949, -0.9330184493337836], [-2.9263004794535545, 5.026437883117941], [0.80645416145117, 3.758828723524164], [3.284929090923728, -3.6227000661580258], [-3.447836059212689, -3.647310284015568], [-1.625662257384021, 3.072709815013726], [0.7602769424162704, 4.794024860706607], [-2.2847585561416968, -0.517625690763838], [0.6856737997305061, 3.1792462708961686], [7.434487717013161, 1.9439019283071575], [-2.388321304737894, -1.899416395264461], [0.6252479500513946, 3.369328418717423], [-1.004824462893196, 5.761855273203045], [3.411500943108601, -2.1786472181528787], [0.895632244853948, -1.9381317634703237], [0.8557814888026711, -3.3262602285629748], [-0.23577428717797821, -3.4613930312339027], [2.2984500432385984, 0.2401807812854364], [0.9929487664755621, -0.7603012580415649], [3.9735309321486705, 4.517661732296059], [-0.06734888812098228, -2.750813451991558], [3.1752095910261064, 0.9315615053102624], [-0.026244444512734, -2.728425205461482], [-2.3089666688277997, 8.512834150279348], [-0.9798046514458889, 3.318884096997367], [1.2691900121950108, 2.3290013721534812], [-0.3295679748583097, -1.3852726642934186], [5.065812645471419, -3.181915702213119], [0.9967362401567141, -3.4265072249079513], [-4.6615177233539065, 4.191631712184568], [-3.6756680914878666, 4.365885101440132], [-1.6044671511492943, 2.4026543610193816], [-3.4870059930214214, 1.4797639765420691], [-1.4061910304529295, 2.5507576873009636], [-0.49588655199546927, -0.7718544057641696], [0.31812248876170823, 3.8808675439605653], [0.8297651194240518, 1.1034164139441218], [0.955685077915639, -1.5157001184057073], [1.8606548559009104, -2.9794615857147617], [-7.024947741623986, 1.522061884819359], [2.731956552179851, -3.06780587632241], [0.41716496504854017, -1.691338524499796], [-4.576483535748092, -0.3337037330114587], [-8.676712082318103, -0.19111682462656088], [-0.843767089393642, -1.3400023244380612], [0.9218786743708282, 1.0126582359750793], [-1.0143251878634154, 1.188754004214144], [0.2806049374277204, -4.430459397598726], [-3.3704337274762484, -1.691201925438966], [0.40900347579195095, -0.5204500005128022], [0.2190568893603247, -0.20904850973844363], [-1.1245110690995426, 1.096377357434669], [-0.2972479882775855, 1.4440529160443347], [3.6352137258859685, 0.6692525456857008], [-0.83448505499437, -3.478006844512695], [-5.975529230813983, -5.446779021576986], [1.9480598050090485, -1.459692330763727], [-1.5793475649127857, 1.8953604331776843], [-1.051268733938766, -1.5563922957477605], [2.7055153671479326, 3.5996417861276298], [-1.6945531224440167, -2.2492968341614965], [-6.396378137877722, 4.7435536584544336], [6.115917787317734, 2.211656076002754], [-4.098109410858985, 4.734428742541122], [8.025240069344465, 5.19664002926344], [1.5711890449539099, 7.225795964291294], [1.4568608542881973, 0.5214991953781499], [0.10644490460908843, 0.2787937294079582], [4.555357276896172, 0.8261472271107984], [-4.03032262225062, -0.23569698523079152], [-0.1440683864790067, -0.0371783175326322], [-0.08869216199854768, 1.9774893189749678], [-0.4939406287573891, -5.062897437084835], [3.286740574147082, -0.2442005475869966], [3.1460674470508687, -4.08913639326025], [-4.195191850920548, 3.619360827313388], [-0.06533450721116435, -0.872911940151411], [-1.1120274875897798, 1.6992612421171176], [1.1595191152242585, 0.8163290947290789], [2.562523036898191, -2.428415569295686], [1.0525434696343794, -1.6587306494640617], [3.949853363491484, -3.448051574307901], [-4.898777869231424, -0.5068961237823801], [3.804126963250541, -3.277498044839827], [1.6778773851147775, -2.6711386526611856], [-2.900680891033825, -1.0387419128075335], [-4.848605045148634, 1.8898063518521453], [-2.502675276990972, -2.0535025006104544], [0.5499822935629051, 2.9646298126471278], [-7.338907640997098, 0.0409109530153291], [1.1332902114965717, -0.5468540934107438], [0.16782928688640625, 0.7904767586653285], [1.2582474073794312, 3.0036748546765084], [-0.060684638140671154, -5.053844672495908], [0.9567293470016475, 1.1601246604631166], [-2.8149165015429323, -1.715933191498082], [3.503878199350167, -1.2699697672984964], [2.6067848162696183, -6.063286021185606], [3.7337356915825195, 2.9595867065597177], [-1.3250323138536988, 4.262724531276803], [0.6644042107692495, -0.8352146276981922], [-1.118101489819965, -1.7698967696742673], [2.6978724771433944, 2.6726209521864996], [3.479962638287302, -0.7386315864145241], [-3.2997577760130223, -2.700179680457352], [2.158909231190396, -2.411481523612787], [1.8359060462975636, -2.5338795853408027], [-2.8168490942868356, -1.6877315173363834], [5.927835811346032, 2.4666987802697182], [-1.08505989013212, 2.4804614069323674], [0.7831132294798451, -0.47206045038222133], [0.5537611071882644, -1.9306913340835448], [2.705853755959905, 1.7762306202558578], [-1.2475770658729126, -0.24937838151687836], [1.609192961725411, -6.182121058659931], [0.7438527245790993, 0.4988183863315878], [0, 0], [0, 0], [0, 0], [0, 0]]
color = []
m = []
z = 0
dp = 0
k = 0
MFP = []
"""set initial positions and velocities"""

for j in range(number_of_atoms):
    position.append([edgeX+random.randint(0,boxWidth),edgeY+random.randint(0,boxHeight),0])
    v.append([random.uniform(-5,5),random.uniform(-5,5)])
    m.append(1)
    color.append((0,0,255))
color[number_of_atoms-1]=((250,0,0))
for i in range(4):
 v.append([0,0])

def colission(a,v,b): # a - list of all atoms, v - list of all velocites (atom j is in position position[j] and travel with velocity v[j]), b - the atom we check for collisions
 for ball in a:
  if ball is not b:
   distanceX,distanceY = (b[0]-ball[0]),(b[1]-ball[1])
   distance = distanceX*distanceX + distanceY*distanceY
   if distance < 36: #radius=3,(2radius)**2=36
    ma,mb=1,1 #the mass of all atoms is 1
    i,j = a.index(ball),a.index(b)
    if (i  == number_of_atoms-1 or j  == number_of_atoms-1) and (i not in MFP and j not in MFP):
        MFP.append(position[i])
    change = [-0.15*distanceX,-0.15*distanceY]

    Vax,Vay = v[i][0],v[i][1]
    Vbx,Vby = v[j][0],v[j][1]

    Angle =math.atan2(distanceY,distanceX) # math.acos((distanceX/distance))
    cos,sin = math.cos(Angle),math.sin(Angle)

    #changing coordiantes:
    vg,vh=[Vag,Vah],[Vbg,Vbh]=[Vax*cos+Vay*sin,Vay*cos-Vax*sin],[Vbx*cos+Vby*sin,Vby*cos-Vbx*sin]

    #after collision:
    Va2,Vb2=(2*mb*Vbg)/(ma+mb), (2*ma*Vag)/(ma+mb) #because ma=mb the term (ma-mb)*Vag or (mb-ma)*Vbg  cancels out

    v1,v2 = [Va2,Vah] , [Vb2,Vbh]
    #changing coordiantes again:
    sin = -sin
    v1,v2 =  [(Va2*cos+Vah*sin),Vah*cos-Va2*sin], [v2[0]*cos+v2[1]*sin,v2[1]*cos-v2[0]*sin]
    v[j],v[i]= v2,v1 #update velocites
    a[i][0],a[i][1]=a[i][0] + change[0], a[i][1] + change[1] #treating overlaps
    return()

"""the loop"""

while 1:
    k+=1
    screen.fill((0,0,0))

    mouseX,mouseY=pygame.mouse.get_pos()
    for event in pygame.event.get():
     if event.type == pygame.QUIT:
        pygame.quit()
        sys.exit()
     elif event.type == MOUSEBUTTONDOWN:
        if mouseY<edgeY+50 and mouseY>edgeY-50 and mouseX<edgeX+boxWidth and mouseX>edgeX-50:
         z=1
     elif event.type == MOUSEBUTTONUP:
        z=0
    if z==1:
        dy = mouseY-edgeY
        boxHeight = boxHeight-dy
        edgeY=mouseY

    pygame.draw.rect(screen, (200,100,0), (edgeX,edgeY,boxWidth,boxHeight),1)
    for event in pygame.event.get():
     if event.type == pygame.QUIT:
        sys.exit()

    #checking for the temperature. Can be any other integer rather than 300
    if k%300 == 0:
      v_rms = math.sqrt(sum(velocity[0]**2 + velocity[1]**2 for velocity in v[:number_of_atoms])/number_of_atoms) #Vrms = sqrt(sum(Vx**2+Vy**2)/num of atoms)
      EK = 0.5*1*(v_rms)**2
      print "T=",EK

    # checking the pressure:
      p = dp/(300*boxHeight)
      print "P=",p
      dp = 0

    for j in range(len(position[:number_of_atoms])):
        colission(position,v,position[j]) #check for collisions
        position[j][0]+=v[j][0]
        position[j][1]+=v[j][1]

        """checking for collisions with the walls"""

        if position[j][0]+radius>edgeX+boxWidth:
        # calculating pressure
         dp+=2*v[j][0]

         v[j][0]=-1*(v[j][0])
         position[j]=[(-radius+edgeX+boxWidth)-2,position[j][1],0]

        if position[j][0] - radius<edgeX:

         v[j][0]=-1*(v[j][0])
         position[j]=[radius+edgeX+2,position[j][1],0]

        if position[j][1]+radius>edgeY+boxHeight:

         v[j][1]=-1*(v[j][1])
         position[j]=[position[j][0],(-radius+edgeY+boxHeight)-2,0]

        if position[j][1]-radius<edgeY:

         v[j][1]=-1*(v[j][1])
         position[j]=[position[j][0],edgeY+2+radius,0]
    #checking foe mean free path
    if k%1000 == 0:
     MFP = np.array(MFP[len(MFP)-100:len(MFP)-1])
     d = []
     for i in range(1,len(MFP)):
      d.append(MFP[i]-MFP[i-1]) #d=[x,y] when x = horizontal distance between 2 collision, y =vertical distance between 2 collision
     distance = []
     for difference in d:
      distance.append(math.sqrt((difference[0])**2 + (difference[1])**2))
     print "MFP=",np.average(distance)
     MFP=MFP.tolist()

    clock.tick(40)
    for i in range(len(position[:number_of_atoms])):
     circle = pygame.draw.circle(screen,color[i], (int(position[i][0]),int(position[i][1])),radius)
    pygame.display.flip()