Skip to content

Instantly share code, notes, and snippets.

@j-faria
Created February 28, 2019 00:32
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save j-faria/a2e487041d3d8f5244c65a02ca3bf2f2 to your computer and use it in GitHub Desktop.
Save j-faria/a2e487041d3d8f5244c65a02ca3bf2f2 to your computer and use it in GitHub Desktop.
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