Triginometric formula for calculating results of 2D elastic collision I am trying to write a code that manages elastic collisions between two circles in 2D space but am failing to locate a working formula. 
For the collision itself the known variables are:


*

*mass of both particles

*initial velocities of both particles

*initial movement angles for both particles

*contact angle


What I try to establish are the final velocities and movement angles of colliding particles.
According to my research, pretty much any source that attempts to solve this problem refers to the below formula available on Wikipedia (link)
\begin{align}
            v'_{1x}&=\frac{v_{1}\cos(\theta_1-\varphi)(m_1-m_2)+2m_2v_{2}\cos(\theta_2-\varphi)}{m_1+m_2}\cos(\varphi)
            \\[0.2em]
            &\quad+v_{1}\sin(\theta_1-\varphi)\cos(\varphi+\frac{\pi}{2})
            \\[0.8em]
            v'_{1y}&=\frac{v_{1}\cos(\theta_1-\varphi)(m_1-m_2)+2m_2v_{2}\cos(\theta_2-\varphi)}{m_1+m_2}\sin(\varphi)
            \\[0.2em]
            &\quad+v_{1}\sin(\theta_1-\varphi)\sin(\varphi+\frac{\pi}{2})
            \end{align}
The formula itself is taken from this website. And has apparently already caused issues to users of this stack exchange previously (a link to a previous question using this formula)
The issue that I have is that replicating this formula in my code does not always provide what appears to be a correct result, but it works as expected in some of the cases. 
I was verifying the expected outcomes against the simulation available here. The results produced by this simmulation differ from the ones produced by the formula, even though the same formula is also displayed on the website hosting the simulation. 
My question is - is this formula actually correct and if it is not what would the correct formula look like? Also, I know it is possible to solve this with vectors, but I am looking for a trigonometric solution.
Thank you

UPD: some examples of issues
Example 1 (appears to operate as expected):
m1 = 100; m2 = 10; v1 = 6; v2 = 1; angle1 = 0; angle2 = 0; contact_angle = 0
results given:


*

*v1x = 5.090909

*v1y = 0

*v2x = 10.09091

*v2y = 0


Example 2 (appears to operate as expected):
m1 = 30; m2 = 20; v1 = 40; v2 = 50; angle1 = 315; angle2 = 135; contact_angle = -45
results given:


*

*v1x = -22.62742

*v1y = 22.62742

*v2x = 41.01219

*v2y = -41.01219


Example 3 (appears to be wrong):
m1 = 10; m2 = 10; v1 = 6; v2 = 1; angle1 = 0; angle2 = 180; contact_angle = -45
results given:


*

*v1x = 2.581109

*v1y = -2.416636

*v2x = 2.486482

*v2y = -2.513894


Surely they cannot both go down post collision?

UPD2: Example 3 was produced by an incorrect code. My reply below includes a version that produces correct results
 A: The formula appears to be correct or at least produces same results as the simmulation on "scincecalculators" (link) that I was using as a reference.
The issue in my Example 3 was due to coding error.
Here is a Python code that produces the correct results (link to code snippet), hopefully somebody will find this useful:
import math

mass1 = 10
mass2 = 10
velocity1 = 6
velocity2 = 1
angle1 = 0
angle2 = 180
contact_angle = -45
pi = math.pi

def final_vx(m1, m2, v1, v2, O1, O2, phi):
  numerator_p1 = v1 * math.cos((O1 - phi) * pi / 180) * (m1 - m2)
  numerator_p2 = 2 * m2 * v2 * math.cos((O2 - phi) * pi / 180)
  numerator = numerator_p1 + numerator_p2
  denominator = m1 + m2
  addition = v1 * math.sin((O1 - phi) * pi / 180) * math.cos((phi * pi / 180) + pi / 2)
  result = (numerator / denominator) * math.cos(phi * pi / 180) + addition

  return(result)


def final_vy(m1, m2, v1, v2, O1, O2, phi):
  numerator_p1 = v1 * math.cos((O1 - phi) * pi / 180) * (m1 - m2)
  numerator_p2 = 2 * m2 * v2 * math.cos((O2 - phi) * pi / 180)
  numerator = numerator_p1 + numerator_p2
  denominator = m1 + m2
  addition = v1 * math.sin((O1 - phi) * pi / 180) * math.sin((phi * pi / 180) + pi / 2)
  result = (numerator / denominator) * math.sin(phi * pi / 180) + addition

  return(result)

v1x = final_vx(mass1, mass2, velocity1, velocity2, angle1, angle2, contact_angle)
v1y = final_vy(mass1, mass2, velocity1, velocity2, angle1, angle2, contact_angle)
v2x = final_vx(mass2, mass1, velocity2, velocity1, angle2, angle1, contact_angle)
v2y = final_vy(mass2, mass1, velocity2, velocity1, angle2, angle1, contact_angle)
print("v1x = " + str(v1x))
print("v1y = " + str(v1y))
print("v2x = " + str(v2x))
print("v2y = " + str(v2y))

