import numpy as np
import random
import math
#first I set up the boundaries of the container
max_w = 100
min_w = -100
max_h = 100
min_h = -100
current_time = 0
#now the big particle
mass = 1 #0
radius = 1 #1
bigx_loc = 0 #2
bigy_loc = 0 #3
bigv_x = 0 #4
bigv_y = 0 #5
big_particle = np.array([mass,radius,bigx_loc,bigy_loc,bigv_x,bigv_y])
The first function populates the area with a bunch of particles. It gives each a random x coordinate, y coordinate, x velocity, and y velocity.
def populate_box():
#this function returns arrays of particles that populate the box
#first we set up empty arrays
pos_x = np.array([])
pos_y = np.array([])
v_x = np.array([])
v_y = np.array([])
#next I want to define the distribution
for i in range(100):
particlex = random.randint(2,10)
if random.getrandbits(1):
particlex = -particlex
pos_x = np.append(pos_x, particlex)
particley = random.randint(2,10)
if random.getrandbits(1):
particley = -particley
pos_y = np.append(pos_y, particley)
#these are test velicities, not how I will actually determine velocities
part_vx = np.random.normal(0,.1)
v_x = np.append(v_x, part_vx)
part_vy = np.random.normal(0,.1)
v_y = np.append(v_y, part_vy)
return pos_x, pos_y, v_x, v_y
The time step function takes all the particles positions, their velocities, the current time, and a variable $\Delta t$ and moves the positions and updates the time of the system based on the velocities and the $\Delta t$. Works forwards as well as backwards! If a particle goes past the boundary, it just it just comes out the other end, like pacman.
test for particle moving x direction only, y direction only, particle moving across positive x boundary, particle moving across negative x boundary, same but for y boundry, particle moving across corner, particles moving differint directions
def time_step(x,y,x_velocity,y_velocity, time, delta_t,big):
"""given particle positions and velocities and current time return new positions and new time
>>> time_step([0.0],[0.0],[3.0],[0.0],0.0,1.0,[0.0,0.0,0.0,0.0,0.0,0.0])
([3.0], [0.0], [3.0], [0.0], 1.0, [0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
>>> time_step([0.0],[0.0],[0.0],[3.0],0.0,2.0,[0.0,0.0,0.0,0.0,0.0,0.0])
([0.0], [6.0], [0.0], [3.0], 2.0, [0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
>>> time_step([9.0],[0.0],[3.0],[0.0],0.0,1.0,[0.0,0.0,0.0,0.0,0.0,0.0])
([-8.0], [0.0], [3.0], [0.0], 1.0, [0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
"""
for i in range(len(x)):
x[i] = x[i] + x_velocity[i]*delta_t
y[i] = y[i] + y_velocity[i]*delta_t
if x[i] > 10:
x[i] = (x[i] -20)
if x[i] <-10:
x[i] = (x[i] + 20)
if y[i] > 10:
y[i] = (y[i] -20)
if y[i] <-10:
y[i] = (y[i] + 20)
big[2] += big[4]*delta_t
big[3] += big[5]*delta_t
if big[2] >10:
big[2] = (big[2] - 20)
if big[2] < -10:
big[2] = big[2] + 20
if big[3] >10:
big[3] = (big[3] - 20)
if big[3] < -10:
big[3] = big[3] + 20
time = time +delta_t
return x, y, x_velocity, y_velocity, time, big
if __name__ == "__main__":
import doctest
doctest.testmod()
This function, when given the position of all particles, checks the distance of each small particle from the big particle via $r = \sqrt{(x_s - x_b)^2 + (y_s - y_b)^2}$ where $x_s$ and $y_s$ is the small particles coordinates and $x_b$ and $y_b$ are the coordinates of the big function. If $r$ is less than the radius of the sphere the function returns the number of whichever particle is inside the big particle. If no particles are inside the big particle the function returns None.
def collision_detect(x, y, big):
'''given positions of particles, detects if collision has occurred and returns collided particle'''
for i in range(len(x)):
distance_from_big = (abs(x[i] - big[2])**2 + abs(y[i] - big[3])**2)
if distance_from_big < big[1]**2:
return i
return None
This next function determines how far back in time one must go in order to get the intersecting particle to the edge of the large particle. The first thing it does is shifts everything into the frame of reference of the big particle.
Once in the frame of reference, we simply need to find the negative time it takes for the particle to move to a position at radius the radius. By the pythagorgean theorem $r^2 = (x +v_x t)^2 + (y + v_y t)^2$ This can be rearranged into $0 = (v_x^2 + v_y^2)t^2 + (2v_x x +2 v_y y)t + (x^2 + y^2 - r^2) $. From this $t$ can be found from the quadratic equation.
Problem with this function. Determinant sometimes negative? Should I use absolute velocity? But I think this means there is another problem
def pointofcollision_finder(offendingparticle, x,y,x_velocity,y_velocity,big):
'''
>>> pointofcollision_finder(0, [0.0], [0.0], [1.0], [0.0], [0.0, 1.0, 0.0, 0.0, 0.0, 0.0])
-1.0
'''
# given a bunch of stuff, return the value of delta t the system has to move in order to find the point of contact
#first I will shift everything into frame of reference of the big particle to make this easier.
shifted_vx = (x_velocity[offendingparticle] - big[4])
shifted_vy = (y_velocity[offendingparticle] - big[5])
shifted_x = (x[offendingparticle] - big[2])
shifted_y = (y[offendingparticle] - big[3])
#okay, now I need to figure out how to calculate delta t such that the position is at radius
#uh I can find that from just the quadratic equation
A = (shifted_vx**2 + shifted_vy**2)
B = (2*shifted_vx*shifted_x + 2*shifted_vy*shifted_y)
C = (shifted_x**2 + shifted_y**2 - big[1]**2)
D = (B**2 - 4*A*C)
if D< 0:
print(shifted_vx, shifted_vy, shifted_x, shifted_y, C)
else:
time_delta = ((-B - math.sqrt(B**2 - 4*A*C)) / 2*A)
return time_delta
if __name__ == "__main__":
import doctest
doctest.testmod()
Well that didn't work. I am going to find a new, better, way of determining backward time
def backwardtime(offendingparticle, x,y,x_velocity,y_velocity,big, timeforward):
shiftvx = x_velocity[offendingparticle] - big[4]
shiftvy = y_velocity[offendingparticle] - big[5]
l = True
while l:
r_sq = (x[offendingparticle] - big[2])**2 + (y[offendingparticle] - big[3])**2
timeforward = timeforward/2
if r_sq < big[1]**2:
x[offendingparticle] = x[offendingparticle] - shiftvx * timeforward
y[offendingparticle] = y[offendingparticle] - shiftvy * timeforward
if r_sq > big[1] **2:
x[offendingparticle] = x[offendingparticle] + shiftvx * timeforward
y[offendingparticle] = y[offendingparticle] + shiftvy * timeforward
if r_sq <= big[1]**2 + .001 and r_sq >= big[1]**2:
l = False
negtime = -(timeforward)
return negtime
third time the charm?
def backtime(offendingparticle,x,y,x_velocity,y_velocity,big):
l = 0
while True:
if ((x[offendingparticle] - big[2])**2 + (y[offendingparticle] - big[3])**2) < big[1]**2:
l = l +.01
x[offendingparticle] = x[offendingparticle] - x_velocity[offendingparticle] * .01
y[offendingparticle] = y[offendingparticle] - y_velocity[offendingparticle] * .01
big[2] = big[2] - big[4] *.01
big[3] = big[3] - big[5] * .01
else:
return -l
Once the particles are at the point of collision (when a particle as at the radius of the large particle), this function calculates the new velocities as a result of collision. First, the function calculates the velocity of the center of mass. $v_{cm} = \frac{m_1v_1 + m_2v_2}{m_1 + m_2}$
Next, this is subtracted from the velocity of the big particle and the small particle in order to shift into the frame of reference of the center of mass. Then the shifted velocities are reflected (the collision), and finally shifted back into the normal frame of reference where the new velocities are output to rewrite the old.
Problem with this function. The velocities when output are ridiculously high. Often get overflow error.
def velocity_exchanger(big, x_velocity, y_velocity, offendingparticle):
#given big particle info and offending particle info, adjust new velocities
#first thing that needs to be done is shift into frame of reference for center of mass
mass = 1 #this is mass of small particles
cm_vx = (big[0]*big[4] + x_velocity[offendingparticle]*mass) / (big[0] + mass)
cm_vy = (big[0]*big[5] + y_velocity[offendingparticle]*mass) / (big[0] + mass)
#next we need to shiift into that frame of reference
shift_big_vx = big[4] - cm_vx
shift_big_vy = big[5] - cm_vy
shift_small_vx = x_velocity[offendingparticle] - cm_vx
shift_small_vy = y_velocity[offendingparticle] - cm_vy
#reflect from the actual collision
new_big_vx = -shift_big_vx
new_big_vy = -shift_big_vy
new_small_vx = - shift_small_vx
new_small_vy = -shift_small_vy
#now shift back into standard frame of reference
big[4] = new_big_vx + cm_vx
big[5] = new_big_vy + cm_vy
x_velocity[offendingparticle] = new_small_vx +cm_vx
y_velocity[offendingparticle] = new_small_vy +cm_vy
return big, x_velocity, y_velocity
print(velocity_exchanger([1.0, 0.0, 0.0, 0.0, 0.0, 0.0 ,0.0], [1], [0], 0))
Finally we put all the functions together. First we populate the area. Then we do a time step, followed by a collision check. If the collision check indicates no collision, we go back to the time step and continue on. If the collision check indicates a collision, we use pointofcollision_finder to find how far back in time we need to move the system in order to for the particle that has intersected the large particle to not do that. Then, with this we do a backwards time step. After this step we collision check again. If there is still a particle inside the big particle, we go back to the point of collision finder function. If this collision check indicates no collision, then the velocities are updated. At this point it all goes back to the first time step. This process continues until the time hits a certain value
#test of everything
record = np.array([])
a,b,c,d = populate_box()
time = 0
while time < 500:
last_collision = None
a,b,c,d, time,big_particle = time_step(a,b,c,d,time,1,big_particle)
f = collision_detect(a,b,big_particle)
while f:
print(f)
print("offending coordinates are", a[f], b[f])
print("big coordinates are", big_particle[2], big_particle[3], big_particle[4], big_particle[5])
print("offending velocity are", c[f], d[f])
neg_time = backwardtime(f, a,b,c,d,big_particle, 1) #backtime(f,a,b,c,d,big_particle)#pointofcollision_finder(f,a,b,c,d,big_particle)
print("neg time is", neg_time)
a,b,c,d,time,big_particle = time_step(a,b,c,d,time,neg_time,big_particle)
last_collision = f
print("new coordinates are", a[f], b[f])
f = collision_detect(a,b,big_particle)
print("new bad particle", f)
if f == None:
big_particle,c,d = velocity_exchanger(big_particle,c,d,last_collision)
break
$g(v_x) = \big ( \frac{m}{2 \pi k T} \big)^{\frac{1}{2}} e^{-m v_x^2 / 2 K T}$
(-0.874655236119)**2 + (-0.455879692584)**2
(-0.874656667837)**2 +(-0.455881657088)**2