Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
How to solve the binary mass function for m2, using sympy
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.