""" Handles collision.
.. note::
YF has old module handling ratherford scattering.
I haven't ported it, but can be done on request.
"""
import numpy as _np
[docs]def binary_collision(mass1, pos1, vel1, mass2, pos2, vel2, collision_length):
r""" Calculate the binary collision
:param mass1: Mass of the particle 1
:param pos1: Postion vector of the particle 1
:param vel1: Velocity vector of the particle 1
:param mass2: Mass of the particle 2
:param pos2: Postion vector of the particle 2
:param vel2: Velocity vector of the particle 2
"""
r1 = _np.array(pos1)
v1 = _np.array(vel1)
r2 = _np.array(pos2)
v2 = _np.array(vel2)
# First, translate to the particle 2 stationary frame
r1_st = r1 - r2
v1_st = v1 - v2
# Call binary_collision_stationary_vector1 to get the velocity in the particle 2 stationary frame
v1_st_new, v2_st_new = binary_collision_stationary_vector1(mass1, r1_st, v1_st, mass2, collision_length)
# Back to the original frame
v1_new = v1_st_new + v2
v2_new = v2_st_new + v2
return v1_new, v2_new
[docs]def binary_collision_stationary_vector1(mass1, pos1, vel1, mass2, collision_length):
r""" Slightly generalized version of :func:`binary_collision_stationary`
Calculate the final velocity after the binary collision.
This function accepts the position and velocity vectors of particle 1.
Particle 2 should be origin, and stationary.
If the negative relative velocity (particles 1 and 2 separates as time goes),
a collision will not happen.
Other conditions are hold as :func:`binary_collision_stationary`.
(e.g. Particle 2 should be at origin, and stationary.)
"""
r1 = _np.array(pos1)
v1 = _np.array(vel1)
r1abs = _np.linalg.norm(r1)
xhat = v1 / _np.linalg.norm(v1)
rhat = r1 / r1abs
x_r_cross = _np.cross(xhat, rhat)
x_r_cross_len = _np.linalg.norm(x_r_cross)
if x_r_cross_len != 0:
zhat = x_r_cross / x_r_cross_len
yhat = _np.cross(zhat, xhat)
impact_parameter = r1.dot(yhat)
speed1 = _np.linalg.norm(v1)
# If the particle is moving away, no collision happens
if rhat.dot(xhat) > 0:
return v1, _np.zeros([3])
vx1, vy1, vx2, vy2 = binary_collision_stationary(mass1, speed1, mass2, impact_parameter, collision_length)
# These values are in the stationary frame.
v1_new = vx1 * xhat + vy1 * yhat
v2_new = vx2 * xhat + vy2 * yhat
return v1_new, v2_new
else: # Impact parameter = 0
if rhat.dot(xhat) > 0:
return v1, _np.zeros([3])
impact_parameter = 0
speed1 = _np.linalg.norm(v1)
vx1, vy1, vx2, vy2 = binary_collision_stationary(mass1, speed1, mass2, impact_parameter, collision_length)
v1_new = vx1 * xhat
v2_new = vx2 * xhat
return v1_new, v2_new
[docs]def binary_collision_stationary(mass1, speed1, mass2, impact_parameter, collision_length):
r""" Cancluate a binary collision in one-particle stationary system.
:param mass1: Mass of moving particle, *m1*
:param speed1: Speed of the moving particle, *v1*. Positive.
:param mass2: Mass of stationary particle, *m2*
:param impact_parameter: The impact parameter, *d*. Positive.
:param collision_length: The radius of the cross section, or *r1+r2* where *r1* is the radius of the moving particke
and *r2* is the radius of the stationary particle.
:returns: *(V1, W1, V2, W2)*
A binary collision is simulated here.
Situation is that a flying (incoming) particle 1 has mass1 with speed 1 (positive y component) in the (negative x, positive y) quadrant.
The stationary particle 2 at the origin will be collided.
The returned is the final velocity (x, y components).
If the impact paramter is more than the collision_length, the particle will not collide.
*Theory*
A binary collision is characterized by the two particles with masses, *m1* and *m2*.
Let's be in a system where the second particle with *m2* is in rest at the origin.
This system change be chosen without losing the generality.
First particle (with *m1*) flies from infinite -x, toward the origin, with a speed of *v1* along x axis.
Just by rotation, we can select a frame that the first particle is in the (-x, +y) quadrant, center in z=0 plane,
with v1>0 with only vx component.
Then, after the collision, the motion of two particles should be confined in x-y plane since no force acting
in the z direction. Assume the resulting velocity vector as *(V1, W1)* for particle *m1*, and
*(V2, W2)* for particle *m2*.
Momentum conservation will give you
.. math::
m_1 v_1 = m_1 V_1 + m_2 V_2
and
.. math::
0 = m_1 W_1 + m_2 W_2
Energy conservation will give you
.. math::
m_1 v_1^2 = m_1 V_1^2 + m_1 W_1^2 + m_2 V_2^2 + m_2 W_2^2
The force only act on the collision, so that the particle *m2* should have the following geometric relation of velocity.
.. math::
W_2 = -\frac{d}{\sqrt{(r_1+r_2)^2 - d^2}} V_2
Solving above four equations to derive the *V1*, *W1*, *V2* and *W2*.
For simplicity, we use the following notations.
.. math::
\mu &= \frac{m_2}{m_1} \\
\delta &= \frac{d}{\sqrt{(r_1+r_2)^2 - d^2}} \\
\alpha &= (\mu+1)(1+\delta^2)
The solutions are
.. math::
V_1 &= v_1 (1 - \frac{2\mu}{\alpha}) \\
W_1 &= \frac{2v_1\mu\delta}{\alpha} \\
V_2 &= \frac{2v_1}{\alpha} \\
W_2 &= -\frac{2v_1\delta}{\alpha} \\
.. todo::
To check the above formulation is correct. In particular the fourth equation should be checked.
How to? Probably we may need another formulation using center of gravity system...
"""
d = impact_parameter
r = collision_length
if d >= r:
return (speed1, 0, 0, 0) # No collision, indicating that the velocity vectors do not change
mu = mass2 / mass1
delta2 = d * d / ((r + d) * (r - d))
delta = _np.sqrt(delta2)
v = speed1
alpha = (1 + mu) * (1 + delta2)
v1 = v * (1 - 2 * mu / alpha)
w1 = 2 * v * mu * delta / alpha
v2 = 2 * v / alpha
w2 = -2 * v * delta / alpha
return (v1, w1, v2, w2)