# apps170608_collision.binary1_conservationsΒΆ

from irfpy.util import collision as coll
import numpy as np

def main():
d = 0.2e-10    # impact parameter in meter
r = 1e-10    #  To define the maximum impact parameter, or radius of particles
m1 = 2.5
m2 = 4.8
v = 100    # in m/s

do_check(m1, v, m2, d, r)
do_check(m1, v, m2, r * 0.99, r)
do_check(m1, v, m2, r * 0.01, r)

do_check(m1, v, m2 * 10000, d, r)

do_check(m1 * 10000, v, m2, d, r)

do_check(m1 * 10000, v, m2, d, r / 40.)

def do_check(m1, v, m2, d, r):
vx1, vy1, vx2, vy2 = coll.binary_collision_stationary(m1, v, m2, d, r)

print('#1', vx1, vy1)
print('#2', vx2, vy2)

px_org = m1 * v
px_new = m1 * vx1 + m2 * vx2
pxcons = np.abs((px_new - px_org)/px_org) < 1e-8

py_new = m1 * vy1 + m2 * vy2
pycons = np.abs(py_new) < 1e-8

en_org = m1 * v * v / 2
en_new = (m1 * (vx1 ** 2 + vy1 ** 2) + m2 * (vx2 ** 2 + vy2 ** 2)) / 2
encons = np.abs((en_new - en_org)/en_org) < 1e-8

print('Momentum (x)', px_org, '->', px_new)
if not pxcons: print('!!!! MOMENT_X not conserved !!!!')
print('Momentum (y)', 0, '->', py_new)
if not pycons: print('!!!! MOMENT_Y not conserved !!!!')
print('Energy', en_org, '->', en_new)
if not encons: print('!!!! ENERGY not conserved !!!!')

if __name__ == '__main__':
main()