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()