This is an old revision of the document!

Copy and paste the following code into Glowscript to model the ideal gas law.

GlowScript 2.8 VPython
## Constants
R=8.314 # Gas constant
N_Avo=6.02E23 # Avogodro's constant
L = 0.1 # Our container is a cube with L=0.1m on each side
Vol=L*L*L # Volume of our cube is 1 Liter
Atomic_radius = 0.001 # wildly exaggerated size of a gas atom

container = box(pos=vec(0,0,0),size=vec(L,L,L), color=color.white, opacity=0.1) # Create a container for our gas atoms

N_atoms = 300 # Number of gas atoms in your simulation
n=N_atoms/N_Avo # number of moles of gas atoms

Molar_mass = 4E-3 # molar mass of helium in kg/mol

T = 300 # temperature in Kelvin

#v_avg = sqrt(8*R*T/(pi*Molar_mass)) # average velocity of gas atoms
v_rms = sqrt(3*R*T/(Molar_mass)) # root mean squared velocity of gas atoms
Velocity=v_rms #Which velocity would you like to assign your gas atoms: v_avg or v_rms?

Theoretical_Pressure=n*R*T/Vol # Calculate the theoretical pressure based on the ideal gas law

## The lines below create a graph for the heigh vs time. Try creating a graph for the velocity and acceleration
Grph1 = graph(title='Pressure vs Time', xtitle='Time (s)', ytitle='Pressure (atm)', fast=False, ymin=0, ymax=5*Theoretical_Pressure) #initialize our graphs. Useful boundaries: ymin=0, ymax=5*Theoretical_Pressure
ExperimentalPressureGraph = gcurve(color=color.red, label='Experimental_Pressure') #Make a graph for measured pressure
TheoreticalPressureGraph = gcurve(color=color.blue, label='Theoretical_Pressure') #Make a graph for theoretical pressure
    
t=0 # initialize the time variable
dt = 5E-9 #time-step interval
pressure=0 # initialize the pressure tracker
pcount=0 # This counter will track when a particle-wall collision adds to the pressure

## The following lines create all the particles
ListOfParticles = [] # Empty list of particles
for i in range(0,N_atoms): #Loop over all of the atoms
    newparticle = sphere(pos=vector(L*random()-L/2,L*random()-L/2,L*random()-L/2),radius=Atomic_radius, color=color.green,opacity=0.7, visible = True) #Create a spherical particle at a random postion
    newparticle.velocity = vector(Velocity*sin(pi*random())*cos(2*pi*random()),Velocity*sin(pi*random())*sin(2*pi*random()),Velocity*cos(pi*random())) #Randomize velocity
    newparticle.mass = Molar_mass/N_Avo #Assign particles the atomic mass (in kg/atom) of your gas
    ListOfParticles.append(newparticle) #Append the particle to our list of particles

## The following lines simulate random motion and particle-wall collisions
while True: #Run infinitely
    rate(9999) #Determines how fast the code runs
    for particle in ListOfParticles: #Loop over all partciles
        particle.pos = particle.pos + particle.velocity*dt #Update position of every particle based on its velocity
        
        #Check for particle-wall collisions in x,y,z. If there is a collision, reverse the velocity direction and add some pressure to the total.
        if abs(particle.pos.x) >= L/2:
            particle.velocity.x = - particle.velocity.x 
            pressure+=particle.mass*abs(particle.velocity.x)/(L*L*dt)
            pcount+=1 
        if abs(particle.pos.y) >= L/2:
            particle.velocity.y = - particle.velocity.y  
            pressure+=particle.mass*abs(particle.velocity.y)/(L*L*dt)
            pcount+=1 
        if abs(particle.pos.z) >= L/2:
            particle.velocity.z = - particle.velocity.z 
            pressure+=particle.mass*abs(particle.velocity.z)/(L*L*dt)
            pcount+=1 
    
    #The following lines calculate the experimental pressure measured due to particle-wall collisions
    if pcount>=10: #After at least 10 particle-wall collisions occur
        Experimental_Pressure=pressure/pcount # Average the total pressure over the number of particle-wall collisions
        Experimental_Pressure=Experimental_Pressure/101.325 # Convert from kPa to atm
        ExperimentalPressureGraph.plot(t,Experimental_Pressure) #Graph the experimental pressure
        pressure=0 # Reset the pressure tracker
        pcount=0 # Reset the pressure counter

    TheoreticalPressureGraph.plot(t,Theoretical_Pressure) #Graph theoretical pressure

    t=t+dt #move onto the next time-step