from amuse.lab import *

def critical_velocity_four_bodies(m11, m12, m21, m22, a1, a2):
  mu     = (m11+m12)*(m21+m22)/(m11+m12+m21+m22)
  v_sqrt = (constants.G/mu) * ((m11*m12)/a1 + (m21*m22)/a2)
  v_crit = (v_sqrt).sqrt()
  return v_crit

m11 = 137 | units.MSun
m12 = 129 | units.MSun
m21 = 36.6 | units.MSun
m22 = 10 | units.MSun
a1 = 2 | units.au
a2 = 0.2 | units.au

print(critical_velocity_four_bodies(m11, m12, m21, m22, a1, a2).in_(units.kms))