How to solve the binary mass function for m2, using sympy

mass_function_solve_sympy.py

from sympy import *
from astropy.constants import G
import astropy.units as u

G = G.to(u.m**3 / (u.solMass * u.s**2))

f_rhs = lambda P, K, e: (P * K**3 * (1 - e**2)**(3/2)) / (2 * pi * G.value)
f_lhs = lambda m1, m2, i: (m2**3 * sin(i)) / (m1 + m2)**2

# the earth
P = (365.25*u.day).to(u.s).value
K = 0.09
e = 0.
i = pi/2
# around the Sun
m1 = 1

m2 = symbols('m2')
mp = nsolve(f_rhs(P, K, e) - f_lhs(m1, m2, i), 0.1, solver='halley')

(float(mp)*u.solMass).to(u.earthMass)