https://www.youtube.com/watch?v=HEfHFsfGXjs

approximate_pi.py

from math import pi, sqrt, atan2, asin

def count_collisions(ratio):
    masses = [1, ratio]
    vels = [0, -1]
    collisions = 0
    angles = []

    while vels[0] < 0 or vels[1] < vels[0]:
        collisions += 1

        if vels[0] < 0:
            vels[0] *= -1
            angles.append(atan2(sqrt(masses[1])*vels[1], sqrt(masses[0])*vels[0]))
            '''print(vels[0], vels[1])
            print("Angle", angles[-1])
            if len(angles) > 1: print("Difference", angles[-1]-angles[-2], 2/sqrt(masses[1]))'''
        else:
            '''
            v0'-v1'=v1-v0
             -> v0'=v1-v0+v1'
            m0v0'+m1v1'=m0v0+m1v1
             -> m0v1-m0v0+(m0+m1)v1'=m0v0+m1v1
             -> (m0+m1)v1'=2*m0v0+(m1-m0)v1

             -> v1'=v0-v1+v0'
             -> (m0+m1)v0'+m1v0-m1v1=m0v0+m1v1
             -> (m0+m1)v0'=2*m1v1+(m0-m1)v0
            '''
            vels[0], vels[1] = \
                2*masses[1]*vels[1]+(masses[0]-masses[1])*vels[0], \
                2*masses[0]*vels[0]+(masses[1]-masses[0])*vels[1]
            vels[0] /= masses[0]+masses[1]
            vels[1] /= masses[0]+masses[1]

        '''print("Velocities", vels)
        print("Momentum", masses[0]*vels[0]+masses[1]*vels[1])
        print("Energy", 1/2*masses[0]*vels[0]**2+1/2*masses[1]*vels[1]**2)'''

    print("Collisions:", collisions)
    print("Approximation for pi:", collisions/sqrt(ratio))
    # print(pi*2/asin(2*sqrt(masses[0]*masses[1])/(masses[0]+masses[1])))
    print()

count_collisions(2000000000000)

'''
sqrt(m0)*v0'=(m1-m0)/(m0+m1)*sqrt(m0)*v0-2*sqrt(m0)*sqrt(m1)/(m0+m1)*sqrt(m1)*v1
sqrt(m1)*v1'=2*sqrt(m0)*sqrt(m1)/(m0+m1)*sqrt(m0)*v0+(m1-m0)/(m0+m1)*sqrt(m1)*v1

cos(t)=(m1-m0)/(m0+m1)
sin(t)=2*sqrt(m0*m1)/(m0+m1)
'''