From 7b88c9c9404f947c200a00f2abc72d4b89db4221 Mon Sep 17 00:00:00 2001 From: Will Ceva Date: Tue, 3 Feb 2026 13:39:52 +0100 Subject: [PATCH 01/23] New function get_planet_mass_accurate to get the actual mass of a companion (untested) --- src/kima/pykima/analysis.py | 144 ++++++++++++++++++++++++++++++++++++ 1 file changed, 144 insertions(+) diff --git a/src/kima/pykima/analysis.py b/src/kima/pykima/analysis.py index 2333fef..2b661c0 100644 --- a/src/kima/pykima/analysis.py +++ b/src/kima/pykima/analysis.py @@ -225,6 +225,150 @@ def get_planet_mass(P: Union[float, np.ndarray], K: Union[float, np.ndarray], return m_mj.mean(), m_mj.std(), m_mj else: return (m_mj.mean(), m_mj.std(), m_me.mean(), m_me.std()) + +def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.ndarray], + e: Union[float, np.ndarray], I: Union[float, np.ndarray], star_mass: Union[float, Tuple] = 1.0, + full_output=False): + r""" + Calculate the companion mass, $M_p$, given orbital period `P`, + semi-amplitude `K`, eccentricity `e`, inclination `I`,and stellar mass. If `star_mass` is a + tuple with (estimate, uncertainty), this (Gaussian) uncertainty will be + taken into account in the calculation. Note that this accounts for cases where the + mass of the companion is significant compared to the mass of the star (compared to get_planet_mass). + + Args: + P (Union[float, ndarray]): + orbital period [days] + K (Union[float, ndarray]): + semi-amplitude [m/s] + e (Union[float, ndarray]): + orbital eccentricity + I (Union[float, ndarray]): + inclination [radians] + star_mass (Union[float, Tuple]): + stellar mass, or (mass, uncertainty) [Msun] + + This function returns different results depending on the inputs. + + !!! note "If `P`, `K`, `e`, and `I` are floats and `star_mass` is a float" + + Returns: + m (float): planet mass, in $M_{\rm Jup}$ + m (float): planet mass, in $M_{\rm Earth}$ + + !!! note "If `P`, `K`, `e`, and `I` are floats and `star_mass` is a tuple" + + Returns: + m (tuple): planet mass and uncertainty, in $M_{\rm Jup}$ + m (tuple): planet mass and uncertainty, in $M_{\rm Earth}$ + + !!! note "If `P`, `K`, `e`, and `I` are arrays and `full_output=True`" + + Returns: + m (float): + posterior mean for the planet mass, in $M_{\rm Jup}$ + s (float): + posterior standard deviation for the planet mass, in $M_{\rm Jup}$ + M (array): + posterior samples for the planet mass, in $M_{\rm Jup}$ + + !!! note "If `P`, `K`, `e`, and `I` are arrays and `full_output=False`" + + Returns: + m (float): + posterior mean for the planet mass, in $M_{\rm Jup}$ + s (float): + posterior standard deviation for the planet mass, in $M_{\rm Jup}$ + m (float): + posterior mean for the planet mass, in $M_{\rm Earth}$ + s (float): + posterior standard deviation for the planet mass, in $M_{\rm Earth}$ + """ + from astropy.constants import G + from astropy import units as u + + from sympy.solvers import solve + from sympy import Symbol + + from uncertainties import ufloat + + ms = u.meter / u.second + C = (ms * u.day**(1/3) * u.solMass**(2/3) / (2*np.pi*G)**(1/3)).to(u.solMass).value #have to work in units of solar mass, + #since the units of the primary and secondary mass must be the same + if isinstance(P, float): + C = float(C) + + try: + # calculate for one value of the orbital period + P = float(P) + # then K, e, and I should also be floats + try: + K, e, I = float(K), float(e), float(I) + except TypeError: + raise TypeError("K, e, and I should be floats if P is a float") + + #defining the main coefficient of the cubic equation (comprised of the provided orbital parameter values) + D = P * ( (C * K * np.sqrt(1 - e**2)) / np.sin(I) )**3 + + if isinstance(star_mass, tuple) or isinstance(star_mass, list): + + #first, converting the star_mass to a ufloat to take into account the uncertainty using the uncertainties package + star_mass_ufloat = ufloat(star_mass[0], star_mass[1]) + + #writing out the cubic equation to be solved for the companion mass via sympy + cub_eq = m_sol**3 - (D * m_sol**2) - (D * 2 * star_mass_ufloat * m_sol) - (D * star_mass_ufloat**2) + + #solving the cubic equation for the companion mass + m_sol = solve(cub_eq, Symbol('m_sol')) + + #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTSS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES + m_mj = (m_sol.n * u.solMass).to(u.jupiterMass).value + m_mj_err = (m_sol.s * u.solMass).to(u.jupiterMass).value + + m_me = m_mj * mjup2mearth + m_me_err = m_mj_err * mjup2mearth + return (m_mj, m_mj_err), (m_me, m_me_err) + else: + cub_eq = m_sol**3 - (D * m_sol**2) - (D * 2 * star_mass * m_sol) - (D * star_mass**2) + + #solving the cubic equation for the companion mass + m_sol = solve(cub_eq, Symbol('m_sol')) + + #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTSS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES + m_mj = (m_sol * u.solMass).to(u.jupiterMass).value + + m_me = m_mj * mjup2mearth + return m_mj, m_me + + except TypeError: + # calculate for an array of periods + P = np.atleast_1d(P) + if isinstance(star_mass, tuple) or isinstance(star_mass, list): + # include (Gaussian) uncertainty on the stellar mass + star_mass = star_mass_samples(*star_mass, P.shape[0]) + star_mass = np.repeat(star_mass.reshape(-1, 1), P.shape[1], axis=1) + elif isinstance(star_mass, np.ndarray): + # use the stellar mass as provided + star_mass = np.atleast_1d(star_mass) + + #defining the main coefficient of the cubic equation (comprised of the provided orbital parameter values) + D = P * ( (C * K * np.sqrt(1 - e**2)) / np.sin(I) )**3 + + #writing out the cubic equation to be solved for the companion mass via sympy + cub_eq = m_sol**3 - (D * m_sol**2) - (D * 2 * star_mass * m_sol) - (D * star_mass**2) + + #solving the cubic equation for the companion mass + m_sol = solve(cub_eq, Symbol('m_sol')) + + #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTSS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES + m_mj = (m_sol * u.solMass).to(u.jupiterMass).value + m_me = m_mj * mjup2mearth + + if full_output: + return m_mj.mean(), m_mj.std(), m_mj + else: + return (m_mj.mean(), m_mj.std(), m_me.mean(), m_me.std()) + def get_planet_mass_GAIA(P: Union[float, np.ndarray], a0: Union[float, np.ndarray], parallax: Union[float, np.ndarray], From 67e56b0969b3ce35e72a7bfaa33813a83a5c3d72 Mon Sep 17 00:00:00 2001 From: Will Ceva Date: Tue, 3 Feb 2026 15:18:43 +0100 Subject: [PATCH 02/23] New function works for single prim mass values, but not mu +/- std tuples --- src/kima/pykima/analysis.py | 37 +++++++++++++++++++++++++++---------- 1 file changed, 27 insertions(+), 10 deletions(-) diff --git a/src/kima/pykima/analysis.py b/src/kima/pykima/analysis.py index 2b661c0..04a93c9 100644 --- a/src/kima/pykima/analysis.py +++ b/src/kima/pykima/analysis.py @@ -288,13 +288,16 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda from astropy import units as u from sympy.solvers import solve - from sympy import Symbol + from sympy import Symbol, symbols + from sympy import lambdify from uncertainties import ufloat + import uncertainties.umath as um ms = u.meter / u.second C = (ms * u.day**(1/3) * u.solMass**(2/3) / (2*np.pi*G)**(1/3)).to(u.solMass).value #have to work in units of solar mass, - #since the units of the primary and secondary mass must be the same + #since the units of the primary and secondary mass + # must be the same in the cubic equation if isinstance(P, float): C = float(C) @@ -312,30 +315,43 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda if isinstance(star_mass, tuple) or isinstance(star_mass, list): + raise NotImplementedError('Currently not working - uncertainty propagation issue with sympy and uncertainties packages') + #first, converting the star_mass to a ufloat to take into account the uncertainty using the uncertainties package star_mass_ufloat = ufloat(star_mass[0], star_mass[1]) #writing out the cubic equation to be solved for the companion mass via sympy - cub_eq = m_sol**3 - (D * m_sol**2) - (D * 2 * star_mass_ufloat * m_sol) - (D * star_mass_ufloat**2) + m_sol = Symbol('m_sol') + + #including functionality for error propagation by first treating the stellar mass as a constant, + #so that the secondary mass is returned in terms of the primary mass, allowing for the uncertainty to be propagated + #afterwards + star_mass_const = symbols('star_mass_const') + + cub_eq = m_sol**3 - (D * m_sol**2) - (D * 2 * star_mass_const * m_sol) - (D * star_mass_const**2) #solving the cubic equation for the companion mass - m_sol = solve(cub_eq, Symbol('m_sol')) + m_sol_res_i = solve(cub_eq, m_sol)[0] + + #now propagating the uncertainty by substituting back in the star_mass ufloat for the star_mass constant (doesn't work???) + m_sol_res = lambdify([star_mass_const], m_sol_res_i, modules=um)(star_mass_ufloat) #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTSS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES - m_mj = (m_sol.n * u.solMass).to(u.jupiterMass).value - m_mj_err = (m_sol.s * u.solMass).to(u.jupiterMass).value + m_mj = (m_sol_res.n * u.solMass).to(u.jupiterMass).value + m_mj_err = (m_sol_res.s * u.solMass).to(u.jupiterMass).value m_me = m_mj * mjup2mearth m_me_err = m_mj_err * mjup2mearth return (m_mj, m_mj_err), (m_me, m_me_err) else: + m_sol = Symbol('m_sol') cub_eq = m_sol**3 - (D * m_sol**2) - (D * 2 * star_mass * m_sol) - (D * star_mass**2) #solving the cubic equation for the companion mass - m_sol = solve(cub_eq, Symbol('m_sol')) + m_sol_res = solve(cub_eq, m_sol)[0] #getting only the first root (the only real value) #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTSS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES - m_mj = (m_sol * u.solMass).to(u.jupiterMass).value + m_mj = (m_sol_res * u.solMass).to(u.jupiterMass).value m_me = m_mj * mjup2mearth return m_mj, m_me @@ -355,13 +371,14 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda D = P * ( (C * K * np.sqrt(1 - e**2)) / np.sin(I) )**3 #writing out the cubic equation to be solved for the companion mass via sympy + m_sol = Symbol('m_sol') cub_eq = m_sol**3 - (D * m_sol**2) - (D * 2 * star_mass * m_sol) - (D * star_mass**2) #solving the cubic equation for the companion mass - m_sol = solve(cub_eq, Symbol('m_sol')) + m_sol_res = solve(cub_eq, m_sol) #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTSS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES - m_mj = (m_sol * u.solMass).to(u.jupiterMass).value + m_mj = (m_sol_res * u.solMass).to(u.jupiterMass).value m_me = m_mj * mjup2mearth if full_output: From cd883e302604a158b6bccbc62c83fd8ec42ef8ca Mon Sep 17 00:00:00 2001 From: Will Ceva Date: Tue, 3 Feb 2026 16:34:30 +0100 Subject: [PATCH 03/23] Try calculating accurate companion mass with full posteriors (fail) --- src/kima/pykima/analysis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/kima/pykima/analysis.py b/src/kima/pykima/analysis.py index 04a93c9..c82272f 100644 --- a/src/kima/pykima/analysis.py +++ b/src/kima/pykima/analysis.py @@ -362,7 +362,7 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda if isinstance(star_mass, tuple) or isinstance(star_mass, list): # include (Gaussian) uncertainty on the stellar mass star_mass = star_mass_samples(*star_mass, P.shape[0]) - star_mass = np.repeat(star_mass.reshape(-1, 1), P.shape[1], axis=1) + # star_mass = np.repeat(star_mass.reshape(-1, 1), P.shape[1], axis=1) #commented out for now, since just testing the calculation for single companion fits elif isinstance(star_mass, np.ndarray): # use the stellar mass as provided star_mass = np.atleast_1d(star_mass) From 2df7ac7361e66d4b66a652c5b7cec7d4756298a8 Mon Sep 17 00:00:00 2001 From: Will Ceva Date: Tue, 3 Feb 2026 16:56:00 +0100 Subject: [PATCH 04/23] Calculate accurate companion mass with full posteriors via for loops (horrendously slow!) --- src/kima/pykima/analysis.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/kima/pykima/analysis.py b/src/kima/pykima/analysis.py index c82272f..c713194 100644 --- a/src/kima/pykima/analysis.py +++ b/src/kima/pykima/analysis.py @@ -291,6 +291,8 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda from sympy import Symbol, symbols from sympy import lambdify + from tqdm import tqdm + from uncertainties import ufloat import uncertainties.umath as um @@ -375,7 +377,12 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda cub_eq = m_sol**3 - (D * m_sol**2) - (D * 2 * star_mass * m_sol) - (D * star_mass**2) #solving the cubic equation for the companion mass - m_sol_res = solve(cub_eq, m_sol) + #...in a very inefficient way, since for some reason, solving the above equation but with D and star_mass treated + #as constants, and then substituting in the arrays afterwards, doesn't work (get nans, complex numbers, etc...?) + m_sol_res = [] + for cub_eq_i in tqdm(cub_eq): + m_sol_res_i = solve(cub_eq_i, m_sol) + m_sol_res.append(m_sol_res_i[0]) #getting only the first root (the only real value) #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTSS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES m_mj = (m_sol_res * u.solMass).to(u.jupiterMass).value From ddd5dfc6d15f341de92e8b99b282ad5082cdc8e3 Mon Sep 17 00:00:00 2001 From: Will Ceva Date: Wed, 4 Feb 2026 12:50:38 +0100 Subject: [PATCH 05/23] Found vectorized approach to calculate companion mass with full posteriors --- src/kima/pykima/analysis.py | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/src/kima/pykima/analysis.py b/src/kima/pykima/analysis.py index c713194..3eaa501 100644 --- a/src/kima/pykima/analysis.py +++ b/src/kima/pykima/analysis.py @@ -290,6 +290,7 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda from sympy.solvers import solve from sympy import Symbol, symbols from sympy import lambdify + from sympy import roots from tqdm import tqdm @@ -369,23 +370,26 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda # use the stellar mass as provided star_mass = np.atleast_1d(star_mass) + #defining the main coefficient of the cubic equation (comprised of the provided orbital parameter values) D = P * ( (C * K * np.sqrt(1 - e**2)) / np.sin(I) )**3 #writing out the cubic equation to be solved for the companion mass via sympy - m_sol = Symbol('m_sol') - cub_eq = m_sol**3 - (D * m_sol**2) - (D * 2 * star_mass * m_sol) - (D * star_mass**2) + m_ms = Symbol('m_ms') + D_const, star_mass_const = symbols('D_const star_mass_const') + cub_eq = m_ms**3 - (D_const * m_ms**2) - (D_const * 2 * star_mass_const * m_ms) - (D_const * star_mass_const**2) + + m_ms_res_i = list(roots(cub_eq, m_ms))[1] #getting only the root that will provide a real value when substituting with the actual values + + #solving the cubic equation for the companion mass (have to convert the arrays to complex128 to avoid issues with sympy) + D_c = np.asarray(D, dtype=np.complex128) + star_mass_c = np.asarray(star_mass, dtype=np.complex128) - #solving the cubic equation for the companion mass - #...in a very inefficient way, since for some reason, solving the above equation but with D and star_mass treated - #as constants, and then substituting in the arrays afterwards, doesn't work (get nans, complex numbers, etc...?) - m_sol_res = [] - for cub_eq_i in tqdm(cub_eq): - m_sol_res_i = solve(cub_eq_i, m_sol) - m_sol_res.append(m_sol_res_i[0]) #getting only the first root (the only real value) + #need to compare accuracy/precision of solve vs roots!!! + m_ms_res = np.real(lambdify([D_const, star_mass_const], m_ms_res_i, modules=["numpy"])(D_c, star_mass_c)) #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTSS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES - m_mj = (m_sol_res * u.solMass).to(u.jupiterMass).value + m_mj = (m_ms_res * u.solMass).to(u.jupiterMass).value m_me = m_mj * mjup2mearth if full_output: From 4c53228385f51ba228582d29ad891926398352dd Mon Sep 17 00:00:00 2001 From: Will Ceva Date: Wed, 4 Feb 2026 15:24:03 +0100 Subject: [PATCH 06/23] Fix issue with unreal solutions by stating constants are always positive --- src/kima/pykima/analysis.py | 27 ++++++++------------------- 1 file changed, 8 insertions(+), 19 deletions(-) diff --git a/src/kima/pykima/analysis.py b/src/kima/pykima/analysis.py index 3eaa501..c4e8ca3 100644 --- a/src/kima/pykima/analysis.py +++ b/src/kima/pykima/analysis.py @@ -288,11 +288,7 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda from astropy import units as u from sympy.solvers import solve - from sympy import Symbol, symbols - from sympy import lambdify - from sympy import roots - - from tqdm import tqdm + from sympy import symbols, lambdify, Rational from uncertainties import ufloat import uncertainties.umath as um @@ -370,26 +366,19 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda # use the stellar mass as provided star_mass = np.atleast_1d(star_mass) - - #defining the main coefficient of the cubic equation (comprised of the provided orbital parameter values) - D = P * ( (C * K * np.sqrt(1 - e**2)) / np.sin(I) )**3 + #keeping the equation in a form that includes a fractional exponent, in order to invoke sympy.Rational() explicitly, to avoid issues with complex roots + Di = C * P**(1/3) * K * np.sqrt(1 - e**2) / np.sin(I) - #writing out the cubic equation to be solved for the companion mass via sympy - m_ms = Symbol('m_ms') - D_const, star_mass_const = symbols('D_const star_mass_const') - cub_eq = m_ms**3 - (D_const * m_ms**2) - (D_const * 2 * star_mass_const * m_ms) - (D_const * star_mass_const**2) + Di_const, m_sol_i, star_mass_i_const = symbols('Di_const m_sol_i star_mass_i_const', positive=True) - m_ms_res_i = list(roots(cub_eq, m_ms))[1] #getting only the root that will provide a real value when substituting with the actual values + eq = (m_sol_i / ( (m_sol_i + star_mass_i_const)**Rational(2, 3) ) ) - Di_const - #solving the cubic equation for the companion mass (have to convert the arrays to complex128 to avoid issues with sympy) - D_c = np.asarray(D, dtype=np.complex128) - star_mass_c = np.asarray(star_mass, dtype=np.complex128) + m_sol_res_i = solve(eq, m_sol_i)[0].simplify() #getting only the first root, which is the real solution (and then simplifying it) - #need to compare accuracy/precision of solve vs roots!!! - m_ms_res = np.real(lambdify([D_const, star_mass_const], m_ms_res_i, modules=["numpy"])(D_c, star_mass_c)) + m_sol_res = lambdify([Di_const, star_mass_i_const], m_sol_res_i, "numpy")(Di, star_mass) #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTSS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES - m_mj = (m_ms_res * u.solMass).to(u.jupiterMass).value + m_mj = (m_sol_res * u.solMass).to(u.jupiterMass).value m_me = m_mj * mjup2mearth if full_output: From 0647827e47806e3574a41c48e6b05dd0fa46876d Mon Sep 17 00:00:00 2001 From: Will Ceva Date: Wed, 4 Feb 2026 15:48:05 +0100 Subject: [PATCH 07/23] Move companion mass equation solution to top of function --- src/kima/pykima/analysis.py | 54 +++++++++++++------------------------ 1 file changed, 19 insertions(+), 35 deletions(-) diff --git a/src/kima/pykima/analysis.py b/src/kima/pykima/analysis.py index c4e8ca3..27d8897 100644 --- a/src/kima/pykima/analysis.py +++ b/src/kima/pykima/analysis.py @@ -296,10 +296,17 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda ms = u.meter / u.second C = (ms * u.day**(1/3) * u.solMass**(2/3) / (2*np.pi*G)**(1/3)).to(u.solMass).value #have to work in units of solar mass, #since the units of the primary and secondary mass - # must be the same in the cubic equation + # must be the same in the equation defined below if isinstance(P, float): C = float(C) + #getting the general solution for the companion mass in terms of the orbital parameters and the stellar mass + D_const, m_ms_var, star_mass_const = symbols('D_const m_ms_var star_mass_const', positive=True) + eq = (m_ms_var / ( (m_ms_var + star_mass_const)**Rational(2, 3) ) ) - D_const + + m_ms_sol = solve(eq, m_ms_var)[0].simplify() #getting only the first root, which is the real solution (and then simplifying it) + m_ms_func = lambdify([D_const, star_mass_const], m_ms_sol, "numpy") + try: # calculate for one value of the orbital period P = float(P) @@ -309,9 +316,9 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda except TypeError: raise TypeError("K, e, and I should be floats if P is a float") - #defining the main coefficient of the cubic equation (comprised of the provided orbital parameter values) - D = P * ( (C * K * np.sqrt(1 - e**2)) / np.sin(I) )**3 - + #defining the main coefficient of the mass equation to solve (comprised of the provided orbital parameter values) + D = C * P**(1/3) * K * np.sqrt(1 - e**2) / np.sin(I) + if isinstance(star_mass, tuple) or isinstance(star_mass, list): raise NotImplementedError('Currently not working - uncertainty propagation issue with sympy and uncertainties packages') @@ -319,21 +326,8 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda #first, converting the star_mass to a ufloat to take into account the uncertainty using the uncertainties package star_mass_ufloat = ufloat(star_mass[0], star_mass[1]) - #writing out the cubic equation to be solved for the companion mass via sympy - m_sol = Symbol('m_sol') - - #including functionality for error propagation by first treating the stellar mass as a constant, - #so that the secondary mass is returned in terms of the primary mass, allowing for the uncertainty to be propagated - #afterwards - star_mass_const = symbols('star_mass_const') - - cub_eq = m_sol**3 - (D * m_sol**2) - (D * 2 * star_mass_const * m_sol) - (D * star_mass_const**2) - - #solving the cubic equation for the companion mass - m_sol_res_i = solve(cub_eq, m_sol)[0] - - #now propagating the uncertainty by substituting back in the star_mass ufloat for the star_mass constant (doesn't work???) - m_sol_res = lambdify([star_mass_const], m_sol_res_i, modules=um)(star_mass_ufloat) + # error propagation: + m_mj_err = m_mj * (2/3) * star_mass[1] / star_mass[0] #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTSS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES m_mj = (m_sol_res.n * u.solMass).to(u.jupiterMass).value @@ -343,14 +337,10 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda m_me_err = m_mj_err * mjup2mearth return (m_mj, m_mj_err), (m_me, m_me_err) else: - m_sol = Symbol('m_sol') - cub_eq = m_sol**3 - (D * m_sol**2) - (D * 2 * star_mass * m_sol) - (D * star_mass**2) - - #solving the cubic equation for the companion mass - m_sol_res = solve(cub_eq, m_sol)[0] #getting only the first root (the only real value) + m_ms = m_ms_func(D, star_mass) #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTSS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES - m_mj = (m_sol_res * u.solMass).to(u.jupiterMass).value + m_mj = (m_ms * u.solMass).to(u.jupiterMass).value m_me = m_mj * mjup2mearth return m_mj, m_me @@ -366,19 +356,13 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda # use the stellar mass as provided star_mass = np.atleast_1d(star_mass) - #keeping the equation in a form that includes a fractional exponent, in order to invoke sympy.Rational() explicitly, to avoid issues with complex roots - Di = C * P**(1/3) * K * np.sqrt(1 - e**2) / np.sin(I) - - Di_const, m_sol_i, star_mass_i_const = symbols('Di_const m_sol_i star_mass_i_const', positive=True) - - eq = (m_sol_i / ( (m_sol_i + star_mass_i_const)**Rational(2, 3) ) ) - Di_const - - m_sol_res_i = solve(eq, m_sol_i)[0].simplify() #getting only the first root, which is the real solution (and then simplifying it) + #defining the main coefficient of the mass equation to solve (comprised of the provided orbital parameter values) + D = C * P**(1/3) * K * np.sqrt(1 - e**2) / np.sin(I) - m_sol_res = lambdify([Di_const, star_mass_i_const], m_sol_res_i, "numpy")(Di, star_mass) + m_ms = m_ms_func(D, star_mass) #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTSS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES - m_mj = (m_sol_res * u.solMass).to(u.jupiterMass).value + m_mj = (m_ms * u.solMass).to(u.jupiterMass).value m_me = m_mj * mjup2mearth if full_output: From c9fc5c95ad1adafa99d1c827a11c7c87daacb4b9 Mon Sep 17 00:00:00 2001 From: Will Ceva Date: Wed, 4 Feb 2026 16:15:07 +0100 Subject: [PATCH 08/23] get_planet_mass_accurate now working for all single companion calculations scenarios --- src/kima/pykima/analysis.py | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/kima/pykima/analysis.py b/src/kima/pykima/analysis.py index 27d8897..b29405e 100644 --- a/src/kima/pykima/analysis.py +++ b/src/kima/pykima/analysis.py @@ -288,7 +288,7 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda from astropy import units as u from sympy.solvers import solve - from sympy import symbols, lambdify, Rational + from sympy import symbols, lambdify, Rational, diff from uncertainties import ufloat import uncertainties.umath as um @@ -321,17 +321,16 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda if isinstance(star_mass, tuple) or isinstance(star_mass, list): - raise NotImplementedError('Currently not working - uncertainty propagation issue with sympy and uncertainties packages') - - #first, converting the star_mass to a ufloat to take into account the uncertainty using the uncertainties package - star_mass_ufloat = ufloat(star_mass[0], star_mass[1]) + m_ms = m_ms_func(D, star_mass[0]) # error propagation: - m_mj_err = m_mj * (2/3) * star_mass[1] / star_mass[0] - + star_mass_der_sol = diff(m_ms_sol, star_mass_const).simplify() + star_mass_der = lambdify([D_const, star_mass_const], star_mass_der_sol, "numpy")(D, star_mass[0]) + m_ms_err = star_mass_der * star_mass[1] + #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTSS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES - m_mj = (m_sol_res.n * u.solMass).to(u.jupiterMass).value - m_mj_err = (m_sol_res.s * u.solMass).to(u.jupiterMass).value + m_mj = (m_ms * u.solMass).to(u.jupiterMass).value + m_mj_err = (m_ms_err * u.solMass).to(u.jupiterMass).value m_me = m_mj * mjup2mearth m_me_err = m_mj_err * mjup2mearth From 5a243c162d65240c34c05ad84324030dae8468e1 Mon Sep 17 00:00:00 2001 From: Will Ceva Date: Wed, 4 Feb 2026 16:40:22 +0100 Subject: [PATCH 09/23] Move sine into numerator to avoid divide by zero err when i is 0 --- src/kima/pykima/analysis.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/kima/pykima/analysis.py b/src/kima/pykima/analysis.py index b29405e..39ebd41 100644 --- a/src/kima/pykima/analysis.py +++ b/src/kima/pykima/analysis.py @@ -302,7 +302,7 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda #getting the general solution for the companion mass in terms of the orbital parameters and the stellar mass D_const, m_ms_var, star_mass_const = symbols('D_const m_ms_var star_mass_const', positive=True) - eq = (m_ms_var / ( (m_ms_var + star_mass_const)**Rational(2, 3) ) ) - D_const + eq = (( (m_ms_var + star_mass_const)**Rational(2, 3) ) / m_ms_var ) - D_const m_ms_sol = solve(eq, m_ms_var)[0].simplify() #getting only the first root, which is the real solution (and then simplifying it) m_ms_func = lambdify([D_const, star_mass_const], m_ms_sol, "numpy") @@ -317,7 +317,7 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda raise TypeError("K, e, and I should be floats if P is a float") #defining the main coefficient of the mass equation to solve (comprised of the provided orbital parameter values) - D = C * P**(1/3) * K * np.sqrt(1 - e**2) / np.sin(I) + D = np.sin(I) / ( C * P**(1/3) * K * np.sqrt(1 - e**2) ) if isinstance(star_mass, tuple) or isinstance(star_mass, list): @@ -350,13 +350,13 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda if isinstance(star_mass, tuple) or isinstance(star_mass, list): # include (Gaussian) uncertainty on the stellar mass star_mass = star_mass_samples(*star_mass, P.shape[0]) - # star_mass = np.repeat(star_mass.reshape(-1, 1), P.shape[1], axis=1) #commented out for now, since just testing the calculation for single companion fits + star_mass = np.repeat(star_mass.reshape(-1, 1), P.shape[1], axis=1) #commented out for now, since just testing the calculation for single companion fits elif isinstance(star_mass, np.ndarray): # use the stellar mass as provided star_mass = np.atleast_1d(star_mass) #defining the main coefficient of the mass equation to solve (comprised of the provided orbital parameter values) - D = C * P**(1/3) * K * np.sqrt(1 - e**2) / np.sin(I) + D = np.sin(I) / ( C * P**(1/3) * K * np.sqrt(1 - e**2) ) m_ms = m_ms_func(D, star_mass) From 7c1bd78604900c945d177de3ad6bcf039efb3561 Mon Sep 17 00:00:00 2001 From: Will Ceva Date: Thu, 5 Feb 2026 10:13:52 +0100 Subject: [PATCH 10/23] Iteratively calculate outer companion masses by including inner companion masses in primary mass --- src/kima/pykima/analysis.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/src/kima/pykima/analysis.py b/src/kima/pykima/analysis.py index 39ebd41..373c7cc 100644 --- a/src/kima/pykima/analysis.py +++ b/src/kima/pykima/analysis.py @@ -328,7 +328,7 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda star_mass_der = lambdify([D_const, star_mass_const], star_mass_der_sol, "numpy")(D, star_mass[0]) m_ms_err = star_mass_der * star_mass[1] - #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTSS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES + #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES m_mj = (m_ms * u.solMass).to(u.jupiterMass).value m_mj_err = (m_ms_err * u.solMass).to(u.jupiterMass).value @@ -338,7 +338,7 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda else: m_ms = m_ms_func(D, star_mass) - #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTSS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES + #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES m_mj = (m_ms * u.solMass).to(u.jupiterMass).value m_me = m_mj * mjup2mearth @@ -350,7 +350,6 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda if isinstance(star_mass, tuple) or isinstance(star_mass, list): # include (Gaussian) uncertainty on the stellar mass star_mass = star_mass_samples(*star_mass, P.shape[0]) - star_mass = np.repeat(star_mass.reshape(-1, 1), P.shape[1], axis=1) #commented out for now, since just testing the calculation for single companion fits elif isinstance(star_mass, np.ndarray): # use the stellar mass as provided star_mass = np.atleast_1d(star_mass) @@ -358,9 +357,16 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda #defining the main coefficient of the mass equation to solve (comprised of the provided orbital parameter values) D = np.sin(I) / ( C * P**(1/3) * K * np.sqrt(1 - e**2) ) - m_ms = m_ms_func(D, star_mass) + m_ms = np.empty(P.shape, dtype=float) - #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTSS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES + #implementing a for loop to calculate the mass for the outer companions + #more accurately by including the mass of the inner companion(s) in the star_mass + for comp in range(P.shape[1]): + m_ms_comp = m_ms_func(D[:, comp], star_mass) + m_ms[:, comp] = m_ms_comp + star_mass = np.nansum([star_mass, m_ms_comp], axis=0) + + #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES m_mj = (m_ms * u.solMass).to(u.jupiterMass).value m_me = m_mj * mjup2mearth From 67db58d50b61f9c7ea47dd08d5c1d0139433a9ac Mon Sep 17 00:00:00 2001 From: Will Ceva Date: Thu, 5 Feb 2026 10:44:58 +0100 Subject: [PATCH 11/23] Try to change equation to avoid divide by zero error (does not work) --- src/kima/pykima/analysis.py | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/src/kima/pykima/analysis.py b/src/kima/pykima/analysis.py index 373c7cc..10faf63 100644 --- a/src/kima/pykima/analysis.py +++ b/src/kima/pykima/analysis.py @@ -288,7 +288,7 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda from astropy import units as u from sympy.solvers import solve - from sympy import symbols, lambdify, Rational, diff + from sympy import Symbol, symbols, lambdify, Rational, diff from uncertainties import ufloat import uncertainties.umath as um @@ -301,11 +301,12 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda C = float(C) #getting the general solution for the companion mass in terms of the orbital parameters and the stellar mass - D_const, m_ms_var, star_mass_const = symbols('D_const m_ms_var star_mass_const', positive=True) - eq = (( (m_ms_var + star_mass_const)**Rational(2, 3) ) / m_ms_var ) - D_const + sini_const = Symbol('sini_const', positive=True, real=True) + D_const, m_ms_var, star_mass_const = symbols('D_const m_ms_var star_mass_const', positive=True, real=True) + eq = ( (m_ms_var * sini_const) / ( (m_ms_var + star_mass_const)**Rational(2, 3) ) ) - D_const m_ms_sol = solve(eq, m_ms_var)[0].simplify() #getting only the first root, which is the real solution (and then simplifying it) - m_ms_func = lambdify([D_const, star_mass_const], m_ms_sol, "numpy") + m_ms_func = lambdify([sini_const, D_const, star_mass_const], m_ms_sol, "numpy") try: # calculate for one value of the orbital period @@ -317,15 +318,15 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda raise TypeError("K, e, and I should be floats if P is a float") #defining the main coefficient of the mass equation to solve (comprised of the provided orbital parameter values) - D = np.sin(I) / ( C * P**(1/3) * K * np.sqrt(1 - e**2) ) + D = C * P**(1/3) * K * np.sqrt(1 - e**2) if isinstance(star_mass, tuple) or isinstance(star_mass, list): - m_ms = m_ms_func(D, star_mass[0]) + m_ms = m_ms_func(np.sin(I), D, star_mass[0]) # error propagation: star_mass_der_sol = diff(m_ms_sol, star_mass_const).simplify() - star_mass_der = lambdify([D_const, star_mass_const], star_mass_der_sol, "numpy")(D, star_mass[0]) + star_mass_der = lambdify([sini_const, D_const, star_mass_const], star_mass_der_sol, "numpy")(np.sin(I), D, star_mass[0]) m_ms_err = star_mass_der * star_mass[1] #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES @@ -336,7 +337,7 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda m_me_err = m_mj_err * mjup2mearth return (m_mj, m_mj_err), (m_me, m_me_err) else: - m_ms = m_ms_func(D, star_mass) + m_ms = m_ms_func(np.sin(I), D, star_mass) #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES m_mj = (m_ms * u.solMass).to(u.jupiterMass).value @@ -355,14 +356,14 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda star_mass = np.atleast_1d(star_mass) #defining the main coefficient of the mass equation to solve (comprised of the provided orbital parameter values) - D = np.sin(I) / ( C * P**(1/3) * K * np.sqrt(1 - e**2) ) + D = C * P**(1/3) * K * np.sqrt(1 - e**2) m_ms = np.empty(P.shape, dtype=float) #implementing a for loop to calculate the mass for the outer companions #more accurately by including the mass of the inner companion(s) in the star_mass for comp in range(P.shape[1]): - m_ms_comp = m_ms_func(D[:, comp], star_mass) + m_ms_comp = m_ms_func(np.sin(I[:, comp]), D[:, comp], star_mass) m_ms[:, comp] = m_ms_comp star_mass = np.nansum([star_mass, m_ms_comp], axis=0) From 1f32bf8f9953aec1b7be603c612019fb4a0b0c51 Mon Sep 17 00:00:00 2001 From: Will Ceva Date: Thu, 5 Feb 2026 11:18:44 +0100 Subject: [PATCH 12/23] Switch to sympy predicate nonnegative to try to avoid divide by zero (does not work) --- src/kima/pykima/analysis.py | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/src/kima/pykima/analysis.py b/src/kima/pykima/analysis.py index 10faf63..f4d3bcd 100644 --- a/src/kima/pykima/analysis.py +++ b/src/kima/pykima/analysis.py @@ -288,7 +288,7 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda from astropy import units as u from sympy.solvers import solve - from sympy import Symbol, symbols, lambdify, Rational, diff + from sympy import symbols, lambdify, Rational, diff from uncertainties import ufloat import uncertainties.umath as um @@ -301,12 +301,11 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda C = float(C) #getting the general solution for the companion mass in terms of the orbital parameters and the stellar mass - sini_const = Symbol('sini_const', positive=True, real=True) - D_const, m_ms_var, star_mass_const = symbols('D_const m_ms_var star_mass_const', positive=True, real=True) - eq = ( (m_ms_var * sini_const) / ( (m_ms_var + star_mass_const)**Rational(2, 3) ) ) - D_const + D_const, m_ms_var, star_mass_const = symbols('D_const m_ms_var star_mass_const', nonnegative=True, real=True) + eq = (( (m_ms_var + star_mass_const)**Rational(2, 3) ) / m_ms_var ) - D_const m_ms_sol = solve(eq, m_ms_var)[0].simplify() #getting only the first root, which is the real solution (and then simplifying it) - m_ms_func = lambdify([sini_const, D_const, star_mass_const], m_ms_sol, "numpy") + m_ms_func = lambdify([D_const, star_mass_const], m_ms_sol, "numpy") try: # calculate for one value of the orbital period @@ -318,15 +317,15 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda raise TypeError("K, e, and I should be floats if P is a float") #defining the main coefficient of the mass equation to solve (comprised of the provided orbital parameter values) - D = C * P**(1/3) * K * np.sqrt(1 - e**2) + D = np.sin(I) / ( C * P**(1/3) * K * np.sqrt(1 - e**2) ) if isinstance(star_mass, tuple) or isinstance(star_mass, list): - m_ms = m_ms_func(np.sin(I), D, star_mass[0]) + m_ms = m_ms_func(D, star_mass[0]) # error propagation: star_mass_der_sol = diff(m_ms_sol, star_mass_const).simplify() - star_mass_der = lambdify([sini_const, D_const, star_mass_const], star_mass_der_sol, "numpy")(np.sin(I), D, star_mass[0]) + star_mass_der = lambdify([D_const, star_mass_const], star_mass_der_sol, "numpy")(D, star_mass[0]) m_ms_err = star_mass_der * star_mass[1] #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES @@ -337,7 +336,7 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda m_me_err = m_mj_err * mjup2mearth return (m_mj, m_mj_err), (m_me, m_me_err) else: - m_ms = m_ms_func(np.sin(I), D, star_mass) + m_ms = m_ms_func(D, star_mass) #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES m_mj = (m_ms * u.solMass).to(u.jupiterMass).value @@ -356,14 +355,14 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda star_mass = np.atleast_1d(star_mass) #defining the main coefficient of the mass equation to solve (comprised of the provided orbital parameter values) - D = C * P**(1/3) * K * np.sqrt(1 - e**2) + D = np.sin(I) / ( C * P**(1/3) * K * np.sqrt(1 - e**2) ) m_ms = np.empty(P.shape, dtype=float) #implementing a for loop to calculate the mass for the outer companions #more accurately by including the mass of the inner companion(s) in the star_mass for comp in range(P.shape[1]): - m_ms_comp = m_ms_func(np.sin(I[:, comp]), D[:, comp], star_mass) + m_ms_comp = m_ms_func(D[:, comp], star_mass) m_ms[:, comp] = m_ms_comp star_mass = np.nansum([star_mass, m_ms_comp], axis=0) From 242d18e9bde16dbeb565f54ffc866407e88e9cc4 Mon Sep 17 00:00:00 2001 From: Will Ceva Date: Thu, 5 Feb 2026 11:49:34 +0100 Subject: [PATCH 13/23] Change equation and switch to sympy predicate nonnegative to avoid divide by zero (does not work) --- src/kima/pykima/analysis.py | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/src/kima/pykima/analysis.py b/src/kima/pykima/analysis.py index f4d3bcd..abee806 100644 --- a/src/kima/pykima/analysis.py +++ b/src/kima/pykima/analysis.py @@ -288,7 +288,7 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda from astropy import units as u from sympy.solvers import solve - from sympy import symbols, lambdify, Rational, diff + from sympy import Symbol, symbols, lambdify, Rational, diff from uncertainties import ufloat import uncertainties.umath as um @@ -301,11 +301,12 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda C = float(C) #getting the general solution for the companion mass in terms of the orbital parameters and the stellar mass + sini_const = Symbol('sini_const', nonnegative=True, real=True) D_const, m_ms_var, star_mass_const = symbols('D_const m_ms_var star_mass_const', nonnegative=True, real=True) - eq = (( (m_ms_var + star_mass_const)**Rational(2, 3) ) / m_ms_var ) - D_const + eq = ( (m_ms_var * sini_const) / ( (m_ms_var + star_mass_const)**Rational(2, 3) ) ) - D_const m_ms_sol = solve(eq, m_ms_var)[0].simplify() #getting only the first root, which is the real solution (and then simplifying it) - m_ms_func = lambdify([D_const, star_mass_const], m_ms_sol, "numpy") + m_ms_func = lambdify([sini_const, D_const, star_mass_const], m_ms_sol, "numpy") try: # calculate for one value of the orbital period @@ -317,15 +318,15 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda raise TypeError("K, e, and I should be floats if P is a float") #defining the main coefficient of the mass equation to solve (comprised of the provided orbital parameter values) - D = np.sin(I) / ( C * P**(1/3) * K * np.sqrt(1 - e**2) ) + D = C * P**(1/3) * K * np.sqrt(1 - e**2) if isinstance(star_mass, tuple) or isinstance(star_mass, list): - m_ms = m_ms_func(D, star_mass[0]) + m_ms = m_ms_func(np.sin(I), D, star_mass[0]) # error propagation: star_mass_der_sol = diff(m_ms_sol, star_mass_const).simplify() - star_mass_der = lambdify([D_const, star_mass_const], star_mass_der_sol, "numpy")(D, star_mass[0]) + star_mass_der = lambdify([sini_const, D_const, star_mass_const], star_mass_der_sol, "numpy")(np.sin(I), D, star_mass[0]) m_ms_err = star_mass_der * star_mass[1] #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES @@ -336,7 +337,7 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda m_me_err = m_mj_err * mjup2mearth return (m_mj, m_mj_err), (m_me, m_me_err) else: - m_ms = m_ms_func(D, star_mass) + m_ms = m_ms_func(np.sin(I), D, star_mass) #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES m_mj = (m_ms * u.solMass).to(u.jupiterMass).value @@ -355,14 +356,14 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda star_mass = np.atleast_1d(star_mass) #defining the main coefficient of the mass equation to solve (comprised of the provided orbital parameter values) - D = np.sin(I) / ( C * P**(1/3) * K * np.sqrt(1 - e**2) ) + D = C * P**(1/3) * K * np.sqrt(1 - e**2) m_ms = np.empty(P.shape, dtype=float) #implementing a for loop to calculate the mass for the outer companions #more accurately by including the mass of the inner companion(s) in the star_mass for comp in range(P.shape[1]): - m_ms_comp = m_ms_func(D[:, comp], star_mass) + m_ms_comp = m_ms_func(np.sin(I[:, comp]), D[:, comp], star_mass) m_ms[:, comp] = m_ms_comp star_mass = np.nansum([star_mass, m_ms_comp], axis=0) From 1cbb19be0f48b0946c798c1511b3344f41d27f1b Mon Sep 17 00:00:00 2001 From: Will Ceva Date: Thu, 5 Feb 2026 13:44:49 +0100 Subject: [PATCH 14/23] Original get_planet_mass function also had divide by zero issue - ignoring this error and reverting back to previous commit --- src/kima/pykima/analysis.py | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/src/kima/pykima/analysis.py b/src/kima/pykima/analysis.py index abee806..373c7cc 100644 --- a/src/kima/pykima/analysis.py +++ b/src/kima/pykima/analysis.py @@ -288,7 +288,7 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda from astropy import units as u from sympy.solvers import solve - from sympy import Symbol, symbols, lambdify, Rational, diff + from sympy import symbols, lambdify, Rational, diff from uncertainties import ufloat import uncertainties.umath as um @@ -301,12 +301,11 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda C = float(C) #getting the general solution for the companion mass in terms of the orbital parameters and the stellar mass - sini_const = Symbol('sini_const', nonnegative=True, real=True) - D_const, m_ms_var, star_mass_const = symbols('D_const m_ms_var star_mass_const', nonnegative=True, real=True) - eq = ( (m_ms_var * sini_const) / ( (m_ms_var + star_mass_const)**Rational(2, 3) ) ) - D_const + D_const, m_ms_var, star_mass_const = symbols('D_const m_ms_var star_mass_const', positive=True) + eq = (( (m_ms_var + star_mass_const)**Rational(2, 3) ) / m_ms_var ) - D_const m_ms_sol = solve(eq, m_ms_var)[0].simplify() #getting only the first root, which is the real solution (and then simplifying it) - m_ms_func = lambdify([sini_const, D_const, star_mass_const], m_ms_sol, "numpy") + m_ms_func = lambdify([D_const, star_mass_const], m_ms_sol, "numpy") try: # calculate for one value of the orbital period @@ -318,15 +317,15 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda raise TypeError("K, e, and I should be floats if P is a float") #defining the main coefficient of the mass equation to solve (comprised of the provided orbital parameter values) - D = C * P**(1/3) * K * np.sqrt(1 - e**2) + D = np.sin(I) / ( C * P**(1/3) * K * np.sqrt(1 - e**2) ) if isinstance(star_mass, tuple) or isinstance(star_mass, list): - m_ms = m_ms_func(np.sin(I), D, star_mass[0]) + m_ms = m_ms_func(D, star_mass[0]) # error propagation: star_mass_der_sol = diff(m_ms_sol, star_mass_const).simplify() - star_mass_der = lambdify([sini_const, D_const, star_mass_const], star_mass_der_sol, "numpy")(np.sin(I), D, star_mass[0]) + star_mass_der = lambdify([D_const, star_mass_const], star_mass_der_sol, "numpy")(D, star_mass[0]) m_ms_err = star_mass_der * star_mass[1] #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES @@ -337,7 +336,7 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda m_me_err = m_mj_err * mjup2mearth return (m_mj, m_mj_err), (m_me, m_me_err) else: - m_ms = m_ms_func(np.sin(I), D, star_mass) + m_ms = m_ms_func(D, star_mass) #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES m_mj = (m_ms * u.solMass).to(u.jupiterMass).value @@ -356,14 +355,14 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda star_mass = np.atleast_1d(star_mass) #defining the main coefficient of the mass equation to solve (comprised of the provided orbital parameter values) - D = C * P**(1/3) * K * np.sqrt(1 - e**2) + D = np.sin(I) / ( C * P**(1/3) * K * np.sqrt(1 - e**2) ) m_ms = np.empty(P.shape, dtype=float) #implementing a for loop to calculate the mass for the outer companions #more accurately by including the mass of the inner companion(s) in the star_mass for comp in range(P.shape[1]): - m_ms_comp = m_ms_func(np.sin(I[:, comp]), D[:, comp], star_mass) + m_ms_comp = m_ms_func(D[:, comp], star_mass) m_ms[:, comp] = m_ms_comp star_mass = np.nansum([star_mass, m_ms_comp], axis=0) From 56ff16ae48defbc715a97469a8ce4ac57425c3bd Mon Sep 17 00:00:00 2001 From: Will Ceva Date: Thu, 5 Feb 2026 16:04:51 +0100 Subject: [PATCH 15/23] New functions to get accurate companion semimajor axis (partially tested) --- src/kima/pykima/analysis.py | 143 ++++++++++++++++++++++++++++++++++++ 1 file changed, 143 insertions(+) diff --git a/src/kima/pykima/analysis.py b/src/kima/pykima/analysis.py index 373c7cc..14cbc8b 100644 --- a/src/kima/pykima/analysis.py +++ b/src/kima/pykima/analysis.py @@ -564,6 +564,101 @@ def get_planet_semimajor_axis(P: Union[float, np.ndarray], K=None, return a.mean(), a.std() +def get_planet_semimajor_axis_accurate(P: Union[float, np.ndarray], M_c: Union[float, np.ndarray], + star_mass: Union[float, tuple]=1.0, + full_output=False): + r""" + Calculate the semi-major axis of the companion's orbit given orbital period + `P`, companion mass `M_c`, and stellar mass. + + Args: + P (Union[float, ndarray]): + orbital period [days] + M_c (Union[float, ndarray]): + companion mass [MJup] + star_mass (Union[float, Tuple]): + stellar mass, or (mass, uncertainty) [Msun] + + This function returns different results depending on the inputs. + + !!! note "If `P`, `M_c`, and `star_mass` are floats" + + Returns: + a (float): companion semi-major axis, in AU + + !!! note "If `P` and `M_c` are floats and `star_mass` is a tuple" + + Returns: + a (tuple): semi-major axis and uncertainty, in AU + + !!! note "If `P` and `M_c` are arrays and `full_output=True`" + + Returns: + m_a (float): + posterior mean for the semi-major axis, in AU + s_a (float): + posterior standard deviation for the semi-major axis, in AU + a (array): + posterior samples for the semi-major axis, in AU + + !!! note "If `P` and `M_c` are arrays and `full_output=False`" + + Returns: + m_a (float): + posterior mean for the semi-major axis, in AU + s_a (float): + posterior standard deviation for the semi-major axis, in AU + """ + from astropy.constants import G + from astropy import units as u + + f = (G.to(u.AU**3 / u.solMass / u.day**2) ** (1 / 3)).value + if isinstance(P, float): + f = float(f) + + M_c = M_c * mjup2msun # convert to solar masses + + try: + # calculate for one value of the orbital period + P = float(P) + try: + M_c = float(M_c) + except TypeError: + raise TypeError("M_c should be a float if P is a float") + + if isinstance(star_mass, tuple) or isinstance(star_mass, list): + a = f * (star_mass[0] + M_c)**(1/3) * (P / (2 * np.pi))**(2/3) + # error propagation: + a_err = a * (1/3) * star_mass[1] / (star_mass[0] + M_c) + return a, a_err + else: + a = f * (star_mass + M_c)**(1/3) * (P / (2 * np.pi))**(2/3) + return a + + except TypeError: + # calculate for an array of periods + P = np.atleast_1d(P) + if isinstance(star_mass, tuple) or isinstance(star_mass, list): + # include (Gaussian) uncertainty on the stellar mass + star_mass = star_mass_samples(*star_mass, P.shape[0]) + elif isinstance(star_mass, np.ndarray): + # use the stellar mass as provided + star_mass = np.atleast_1d(star_mass) + + a = np.empty(P.shape, dtype=float) + + #accounting for the masses of the inner companion(s) when calculating + #the semi-major axis of the outer companion(s) + for comp in range(P.shape[1]): + a_comp = f * (star_mass + M_c[:, comp])**(1/3) * (P[:, comp] / (2 * np.pi))**(2/3) + a[:, comp] = a_comp + star_mass = np.nansum([star_mass, M_c[:, comp]], axis=0) + + if full_output: + return a.mean(), a.std(), a + else: + return a.mean(), a.std() + def get_planet_mass_and_semimajor_axis(P, K, e, star_mass=1.0, full_output=False, verbose=False): """ @@ -601,6 +696,54 @@ def get_planet_mass_and_semimajor_axis(P, K, e, star_mass=1.0, return mass, a, star_mass return mass, a +def get_planet_mass_and_semimajor_axis_accurate(P, K, e, I, star_mass=1.0, + full_output=False, verbose=False): + """ + Calculate the companion mass M and the semi-major axis given + orbital period `P`, semi-amplitude `K`, eccentricity `e`, inclination `I`, and stellar mass. + If `star_mass` is a tuple with (estimate, uncertainty), this (Gaussian) + uncertainty will be taken into account in the calculation, and samples from + the star_mass distribution will be returned. + + Units: + P [days] + K [m/s] + e [] + I [radians] + star_mass [Msun] + + Returns: + (M, A) where + M is the output of `get_planet_mass_accurate` + A is the output of `get_planet_semimajor_axis_accurate` + star_mass (np.ndarray): + Gaussian samples if star_mass is a tuple or list + """ + + if verbose: + print('Using star mass = %s solar mass' % star_mass) + + if isinstance(star_mass, tuple) or isinstance(star_mass, list): + # include (Gaussian) uncertainty on the stellar mass + star_mass = star_mass_samples(*star_mass, P.shape[0]) + + mass = get_planet_mass_accurate(P, K, e, I, star_mass, full_output) + + if len(mass) == 2: #this is probably a dumb way to do this + if len(mass[0]) == 2: + M_c = mass[0][0] #ignoring the uncertainty on the planet mass when calculating the semi-major axis, for now + else: + M_c = mass[0] + else: + if full_output: + M_c = mass[-1] + else: + M_c = mass[0] + + a = get_planet_semimajor_axis_accurate(P, M_c, star_mass, full_output) + if isinstance(star_mass, np.ndarray): + return mass, a, star_mass + return mass, a def get_bins(res, start=None, end=None, nbins=100): try: From 2ca108dcbadc433873df01ea28eea49de24cca50 Mon Sep 17 00:00:00 2001 From: Will Ceva Date: Thu, 5 Feb 2026 16:13:42 +0100 Subject: [PATCH 16/23] Change meaning of full_output for new comp mass func and account for nan values in mass and SMA --- src/kima/pykima/analysis.py | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/src/kima/pykima/analysis.py b/src/kima/pykima/analysis.py index 14cbc8b..2e50ba2 100644 --- a/src/kima/pykima/analysis.py +++ b/src/kima/pykima/analysis.py @@ -371,9 +371,9 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda m_me = m_mj * mjup2mearth if full_output: - return m_mj.mean(), m_mj.std(), m_mj + return (np.nanmean(m_mj, axis=0), np.nanstd(m_mj, axis=0), np.nanmean(m_me, axis=0), np.nanstd(m_me, axis=0), m_mj) else: - return (m_mj.mean(), m_mj.std(), m_me.mean(), m_me.std()) + return np.nanmean(m_mj, axis=0), np.nanstd(m_mj, axis=0), m_mj def get_planet_mass_GAIA(P: Union[float, np.ndarray], a0: Union[float, np.ndarray], @@ -655,9 +655,9 @@ def get_planet_semimajor_axis_accurate(P: Union[float, np.ndarray], M_c: Union[f star_mass = np.nansum([star_mass, M_c[:, comp]], axis=0) if full_output: - return a.mean(), a.std(), a + return np.nanmean(a, axis=0), np.nanstd(a, axis=0), a else: - return a.mean(), a.std() + return np.nanmean(a, axis=0), np.nanstd(a, axis=0) def get_planet_mass_and_semimajor_axis(P, K, e, star_mass=1.0, full_output=False, verbose=False): @@ -735,10 +735,7 @@ def get_planet_mass_and_semimajor_axis_accurate(P, K, e, I, star_mass=1.0, else: M_c = mass[0] else: - if full_output: - M_c = mass[-1] - else: - M_c = mass[0] + M_c = mass[-1] a = get_planet_semimajor_axis_accurate(P, M_c, star_mass, full_output) if isinstance(star_mass, np.ndarray): From 8af7df80d52ed8161de211dac5b42812ccfd1b80 Mon Sep 17 00:00:00 2001 From: Will Ceva Date: Thu, 5 Feb 2026 16:39:24 +0100 Subject: [PATCH 17/23] New functions to get accurate comp mass and SMA values are now (basically) fully tested and functional --- src/kima/pykima/analysis.py | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/src/kima/pykima/analysis.py b/src/kima/pykima/analysis.py index 2e50ba2..0c7a44d 100644 --- a/src/kima/pykima/analysis.py +++ b/src/kima/pykima/analysis.py @@ -725,12 +725,19 @@ def get_planet_mass_and_semimajor_axis_accurate(P, K, e, I, star_mass=1.0, if isinstance(star_mass, tuple) or isinstance(star_mass, list): # include (Gaussian) uncertainty on the stellar mass - star_mass = star_mass_samples(*star_mass, P.shape[0]) + star_mass = star_mass_samples(*star_mass, P.shape[0]) #this means that this function + #will break if P is a float and a + #tuple is provided for the star mass, + #and the same applies for the original + #version of this function mass = get_planet_mass_accurate(P, K, e, I, star_mass, full_output) - if len(mass) == 2: #this is probably a dumb way to do this - if len(mass[0]) == 2: + if len(mass) == 2: #this is probably a dumb way to do this - good fix would be: + #move the arrays to the first output position and include the + #planet mass uncertainty in the semi-major axis calculation - + #then the correct mass value would always be in mass[0] + if not isinstance(mass[0], float): M_c = mass[0][0] #ignoring the uncertainty on the planet mass when calculating the semi-major axis, for now else: M_c = mass[0] From e7d92638e5c084f5c4960d29a85519d3d9d2e7c1 Mon Sep 17 00:00:00 2001 From: Will Ceva Date: Fri, 6 Feb 2026 11:43:25 +0100 Subject: [PATCH 18/23] Add companion mass error in error propagation for new sma calculation --- src/kima/pykima/analysis.py | 38 ++++++++++++++++++------------------- 1 file changed, 18 insertions(+), 20 deletions(-) diff --git a/src/kima/pykima/analysis.py b/src/kima/pykima/analysis.py index 0c7a44d..c2b3185 100644 --- a/src/kima/pykima/analysis.py +++ b/src/kima/pykima/analysis.py @@ -371,9 +371,9 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda m_me = m_mj * mjup2mearth if full_output: - return (np.nanmean(m_mj, axis=0), np.nanstd(m_mj, axis=0), np.nanmean(m_me, axis=0), np.nanstd(m_me, axis=0), m_mj) + return (m_mj, np.nanmean(m_mj, axis=0), np.nanstd(m_mj, axis=0), np.nanmean(m_me, axis=0), np.nanstd(m_me, axis=0)) else: - return np.nanmean(m_mj, axis=0), np.nanstd(m_mj, axis=0), m_mj + return m_mj, np.nanmean(m_mj, axis=0), np.nanstd(m_mj, axis=0) def get_planet_mass_GAIA(P: Union[float, np.ndarray], a0: Union[float, np.ndarray], @@ -624,13 +624,22 @@ def get_planet_semimajor_axis_accurate(P: Union[float, np.ndarray], M_c: Union[f try: M_c = float(M_c) except TypeError: - raise TypeError("M_c should be a float if P is a float") + if isinstance(M_c, tuple) or isinstance(M_c, list): + pass + else: + raise TypeError("M_c must be a scalar or tuple/list with shape (2,)") if isinstance(star_mass, tuple) or isinstance(star_mass, list): - a = f * (star_mass[0] + M_c)**(1/3) * (P / (2 * np.pi))**(2/3) - # error propagation: - a_err = a * (1/3) * star_mass[1] / (star_mass[0] + M_c) - return a, a_err + if isinstance(M_c, tuple) or isinstance(M_c, list): + a = f * (star_mass[0] + M_c[0])**(1/3) * (P / (2 * np.pi))**(2/3) + # error propagation: + a_err = np.sqrt(star_mass[1]**2 + M_c[1]**2) * a * (1/3) / (star_mass[0] + M_c[0]) + return a, a_err + else: + a = f * (star_mass[0] + M_c)**(1/3) * (P / (2 * np.pi))**(2/3) + # error propagation: + a_err = a * (1/3) * star_mass[1] / (star_mass[0] + M_c) + return a, a_err else: a = f * (star_mass + M_c)**(1/3) * (P / (2 * np.pi))**(2/3) return a @@ -655,7 +664,7 @@ def get_planet_semimajor_axis_accurate(P: Union[float, np.ndarray], M_c: Union[f star_mass = np.nansum([star_mass, M_c[:, comp]], axis=0) if full_output: - return np.nanmean(a, axis=0), np.nanstd(a, axis=0), a + return a, np.nanmean(a, axis=0), np.nanstd(a, axis=0) else: return np.nanmean(a, axis=0), np.nanstd(a, axis=0) @@ -732,19 +741,8 @@ def get_planet_mass_and_semimajor_axis_accurate(P, K, e, I, star_mass=1.0, #version of this function mass = get_planet_mass_accurate(P, K, e, I, star_mass, full_output) - - if len(mass) == 2: #this is probably a dumb way to do this - good fix would be: - #move the arrays to the first output position and include the - #planet mass uncertainty in the semi-major axis calculation - - #then the correct mass value would always be in mass[0] - if not isinstance(mass[0], float): - M_c = mass[0][0] #ignoring the uncertainty on the planet mass when calculating the semi-major axis, for now - else: - M_c = mass[0] - else: - M_c = mass[-1] - a = get_planet_semimajor_axis_accurate(P, M_c, star_mass, full_output) + a = get_planet_semimajor_axis_accurate(P, mass[0], star_mass, full_output) if isinstance(star_mass, np.ndarray): return mass, a, star_mass return mass, a From 3f9703b74ae9a82d77fe2f0a6b8b0d3e9b9c6b2c Mon Sep 17 00:00:00 2001 From: Will Ceva Date: Fri, 6 Feb 2026 13:14:36 +0100 Subject: [PATCH 19/23] Bug tests and fixes (now works in all parameter combinations) --- src/kima/pykima/analysis.py | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/src/kima/pykima/analysis.py b/src/kima/pykima/analysis.py index c2b3185..3939592 100644 --- a/src/kima/pykima/analysis.py +++ b/src/kima/pykima/analysis.py @@ -353,6 +353,9 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda elif isinstance(star_mass, np.ndarray): # use the stellar mass as provided star_mass = np.atleast_1d(star_mass) + else: + # ensure star_mass is an array of the same shape as P + star_mass = star_mass * np.ones(P.shape[0]) #defining the main coefficient of the mass equation to solve (comprised of the provided orbital parameter values) D = np.sin(I) / ( C * P**(1/3) * K * np.sqrt(1 - e**2) ) @@ -616,16 +619,16 @@ def get_planet_semimajor_axis_accurate(P: Union[float, np.ndarray], M_c: Union[f if isinstance(P, float): f = float(f) - M_c = M_c * mjup2msun # convert to solar masses try: # calculate for one value of the orbital period P = float(P) try: M_c = float(M_c) + M_c = M_c * mjup2msun # convert to solar masses except TypeError: if isinstance(M_c, tuple) or isinstance(M_c, list): - pass + M_c = (M_c[0] * mjup2msun, M_c[1] * mjup2msun) # convert to solar masses else: raise TypeError("M_c must be a scalar or tuple/list with shape (2,)") @@ -653,8 +656,12 @@ def get_planet_semimajor_axis_accurate(P: Union[float, np.ndarray], M_c: Union[f elif isinstance(star_mass, np.ndarray): # use the stellar mass as provided star_mass = np.atleast_1d(star_mass) + else: + # ensure star_mass is an array of the same shape as P + star_mass = star_mass * np.ones(P.shape[0]) a = np.empty(P.shape, dtype=float) + M_c = M_c * mjup2msun # convert to solar masses #accounting for the masses of the inner companion(s) when calculating #the semi-major axis of the outer companion(s) @@ -733,12 +740,11 @@ def get_planet_mass_and_semimajor_axis_accurate(P, K, e, I, star_mass=1.0, print('Using star mass = %s solar mass' % star_mass) if isinstance(star_mass, tuple) or isinstance(star_mass, list): - # include (Gaussian) uncertainty on the stellar mass - star_mass = star_mass_samples(*star_mass, P.shape[0]) #this means that this function - #will break if P is a float and a - #tuple is provided for the star mass, - #and the same applies for the original - #version of this function + if isinstance(P, float): + pass + else: + # include (Gaussian) uncertainty on the stellar mass + star_mass = star_mass_samples(*star_mass, P.shape[0]) mass = get_planet_mass_accurate(P, K, e, I, star_mass, full_output) From 9a14d07bd47ff0a6de01135d141ca7d968ba663c Mon Sep 17 00:00:00 2001 From: Will Ceva Date: Fri, 6 Feb 2026 15:25:40 +0100 Subject: [PATCH 20/23] Switch numerator/denominator of comp mass equation (to always ensure nan results correspond to those where sini = 0) --- src/kima/pykima/analysis.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/kima/pykima/analysis.py b/src/kima/pykima/analysis.py index 3939592..578555f 100644 --- a/src/kima/pykima/analysis.py +++ b/src/kima/pykima/analysis.py @@ -302,7 +302,7 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda #getting the general solution for the companion mass in terms of the orbital parameters and the stellar mass D_const, m_ms_var, star_mass_const = symbols('D_const m_ms_var star_mass_const', positive=True) - eq = (( (m_ms_var + star_mass_const)**Rational(2, 3) ) / m_ms_var ) - D_const + eq = (m_ms_var/( (m_ms_var + star_mass_const)**Rational(2, 3) )) - D_const m_ms_sol = solve(eq, m_ms_var)[0].simplify() #getting only the first root, which is the real solution (and then simplifying it) m_ms_func = lambdify([D_const, star_mass_const], m_ms_sol, "numpy") @@ -317,7 +317,7 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda raise TypeError("K, e, and I should be floats if P is a float") #defining the main coefficient of the mass equation to solve (comprised of the provided orbital parameter values) - D = np.sin(I) / ( C * P**(1/3) * K * np.sqrt(1 - e**2) ) + D = ( C * P**(1/3) * K * np.sqrt(1 - e**2) ) / np.sin(I) if isinstance(star_mass, tuple) or isinstance(star_mass, list): @@ -358,7 +358,7 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda star_mass = star_mass * np.ones(P.shape[0]) #defining the main coefficient of the mass equation to solve (comprised of the provided orbital parameter values) - D = np.sin(I) / ( C * P**(1/3) * K * np.sqrt(1 - e**2) ) + D = ( C * P**(1/3) * K * np.sqrt(1 - e**2) ) / np.sin(I) m_ms = np.empty(P.shape, dtype=float) From 0e5154571d4de27e4c83ca3044aee90da01bf889 Mon Sep 17 00:00:00 2001 From: Will Ceva Date: Fri, 6 Feb 2026 17:19:30 +0100 Subject: [PATCH 21/23] Apply PEP 8 conventions using VSCode plugin --- src/kima/pykima/analysis.py | 245 +++++++++++++++++++++--------------- 1 file changed, 142 insertions(+), 103 deletions(-) diff --git a/src/kima/pykima/analysis.py b/src/kima/pykima/analysis.py index 578555f..21bd7fd 100644 --- a/src/kima/pykima/analysis.py +++ b/src/kima/pykima/analysis.py @@ -13,6 +13,7 @@ from .. import MODELS from .utils import mjup2mearth, mjup2msun, distribution_support + def np_most_probable(results: KimaResults): """ Return the value of Np with the highest posterior probability. @@ -117,7 +118,7 @@ def star_mass_samples(star_mass, star_mass_uncertainty, n): from scipy.stats import truncnorm loc = star_mass scale = star_mass_uncertainty - lower, upper = 0.0, 100.0 # arbitrary upper limit + lower, upper = 0.0, 100.0 # arbitrary upper limit a, b = (lower - loc) / scale, (upper - loc) / scale return truncnorm.rvs(a, b, loc=loc, scale=scale, size=n) @@ -182,7 +183,8 @@ def get_planet_mass(P: Union[float, np.ndarray], K: Union[float, np.ndarray], from astropy.constants import G from astropy import units as u ms = u.meter / u.second - C = (ms * u.day**(1/3) * u.solMass**(2/3) / (2*np.pi*G)**(1/3)).to(u.jupiterMass).value + C = (ms * u.day**(1/3) * u.solMass**(2/3) / + (2*np.pi*G)**(1/3)).to(u.jupiterMass).value if isinstance(P, float): C = float(C) @@ -226,15 +228,19 @@ def get_planet_mass(P: Union[float, np.ndarray], K: Union[float, np.ndarray], else: return (m_mj.mean(), m_mj.std(), m_me.mean(), m_me.std()) + def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.ndarray], - e: Union[float, np.ndarray], I: Union[float, np.ndarray], star_mass: Union[float, Tuple] = 1.0, - full_output=False): + e: Union[float, np.ndarray], Inc: Union[float, np.ndarray], + star_mass: Union[float, Tuple] = 1.0, + full_output=False): r""" - Calculate the companion mass, $M_p$, given orbital period `P`, - semi-amplitude `K`, eccentricity `e`, inclination `I`,and stellar mass. If `star_mass` is a - tuple with (estimate, uncertainty), this (Gaussian) uncertainty will be - taken into account in the calculation. Note that this accounts for cases where the - mass of the companion is significant compared to the mass of the star (compared to get_planet_mass). + Calculate the companion mass, $M_p$, given orbital period `P`, semi + -amplitude `K`, eccentricity `e`, inclination `Inc`,and stellar mass. + If `star_mass` is a tuple with (estimate, uncertainty), this + (Gaussian) uncertainty will be taken into account in the + calculation. Note that this accounts for cases where the mass of + the companion is significant compared to the mass of the star + (compared to get_planet_mass). Args: P (Union[float, ndarray]): @@ -243,46 +249,53 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda semi-amplitude [m/s] e (Union[float, ndarray]): orbital eccentricity - I (Union[float, ndarray]): + Inc (Union[float, ndarray]): inclination [radians] star_mass (Union[float, Tuple]): stellar mass, or (mass, uncertainty) [Msun] This function returns different results depending on the inputs. - !!! note "If `P`, `K`, `e`, and `I` are floats and `star_mass` is a float" + !!! note "If `P`, `K`, `e`, and `Inc` are floats and `star_mass` is + a float" Returns: m (float): planet mass, in $M_{\rm Jup}$ m (float): planet mass, in $M_{\rm Earth}$ - !!! note "If `P`, `K`, `e`, and `I` are floats and `star_mass` is a tuple" + !!! note "If `P`, `K`, `e`, and `Inc` are floats and `star_mass` is + a tuple" Returns: m (tuple): planet mass and uncertainty, in $M_{\rm Jup}$ m (tuple): planet mass and uncertainty, in $M_{\rm Earth}$ - !!! note "If `P`, `K`, `e`, and `I` are arrays and `full_output=True`" + !!! note "If `P`, `K`, `e`, and `Inc` are arrays and + `full_output=True`" Returns: m (float): posterior mean for the planet mass, in $M_{\rm Jup}$ s (float): - posterior standard deviation for the planet mass, in $M_{\rm Jup}$ + posterior standard deviation for the planet mass, in + $M_{\rm Jup}$ M (array): posterior samples for the planet mass, in $M_{\rm Jup}$ - !!! note "If `P`, `K`, `e`, and `I` are arrays and `full_output=False`" + !!! note "If `P`, `K`, `e`, and `Inc` are arrays and + `full_output=False`" Returns: m (float): posterior mean for the planet mass, in $M_{\rm Jup}$ s (float): - posterior standard deviation for the planet mass, in $M_{\rm Jup}$ + posterior standard deviation for the planet mass, in + $M_{\rm Jup}$ m (float): posterior mean for the planet mass, in $M_{\rm Earth}$ s (float): - posterior standard deviation for the planet mass, in $M_{\rm Earth}$ + posterior standard deviation for the planet mass, in + $M_{\rm Earth}$ """ from astropy.constants import G from astropy import units as u @@ -290,45 +303,50 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda from sympy.solvers import solve from sympy import symbols, lambdify, Rational, diff - from uncertainties import ufloat - import uncertainties.umath as um - ms = u.meter / u.second - C = (ms * u.day**(1/3) * u.solMass**(2/3) / (2*np.pi*G)**(1/3)).to(u.solMass).value #have to work in units of solar mass, - #since the units of the primary and secondary mass - # must be the same in the equation defined below + + # have to work in units of solar mass, since the units of the + # primary and secondary mass must be the same in the equation + # defined below + C = (ms * u.day ** (1/3) * u.solMass**(2/3) / + (2*np.pi*G)**(1/3)).to(u.solMass).value + if isinstance(P, float): C = float(C) - #getting the general solution for the companion mass in terms of the orbital parameters and the stellar mass - D_const, m_ms_var, star_mass_const = symbols('D_const m_ms_var star_mass_const', positive=True) - eq = (m_ms_var/( (m_ms_var + star_mass_const)**Rational(2, 3) )) - D_const + # getting the general solution for the companion mass in terms of the orbital parameters and the stellar mass + D_const, m_ms_var, star_mass_const = symbols( + 'D_const m_ms_var star_mass_const', positive=True) + eq = (m_ms_var/((m_ms_var + star_mass_const)**Rational(2, 3))) - D_const - m_ms_sol = solve(eq, m_ms_var)[0].simplify() #getting only the first root, which is the real solution (and then simplifying it) + # getting only the first root, which is the real solution (and then simplifying it) + m_ms_sol = solve(eq, m_ms_var)[0].simplify() m_ms_func = lambdify([D_const, star_mass_const], m_ms_sol, "numpy") try: # calculate for one value of the orbital period P = float(P) - # then K, e, and I should also be floats + + # then K, e, and Inc should also be floats try: - K, e, I = float(K), float(e), float(I) + K, e, Inc = float(K), float(e), float(Inc) except TypeError: - raise TypeError("K, e, and I should be floats if P is a float") - - #defining the main coefficient of the mass equation to solve (comprised of the provided orbital parameter values) - D = ( C * P**(1/3) * K * np.sqrt(1 - e**2) ) / np.sin(I) - + raise TypeError("K, e, and Inc should be floats if P is a float") + + # defining the main coefficient of the mass equation to solve (comprised of the provided orbital parameter values) + D = (C * P**(1/3) * K * np.sqrt(1 - e**2)) / np.sin(Inc) + if isinstance(star_mass, tuple) or isinstance(star_mass, list): - + m_ms = m_ms_func(D, star_mass[0]) # error propagation: star_mass_der_sol = diff(m_ms_sol, star_mass_const).simplify() - star_mass_der = lambdify([D_const, star_mass_const], star_mass_der_sol, "numpy")(D, star_mass[0]) + star_mass_der = lambdify( + [D_const, star_mass_const], star_mass_der_sol, "numpy")(D, star_mass[0]) m_ms_err = star_mass_der * star_mass[1] - - #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES + + # convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES m_mj = (m_ms * u.solMass).to(u.jupiterMass).value m_mj_err = (m_ms_err * u.solMass).to(u.jupiterMass).value @@ -338,7 +356,7 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda else: m_ms = m_ms_func(D, star_mass) - #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES + # convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES m_mj = (m_ms * u.solMass).to(u.jupiterMass).value m_me = m_mj * mjup2mearth @@ -355,21 +373,21 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda star_mass = np.atleast_1d(star_mass) else: # ensure star_mass is an array of the same shape as P - star_mass = star_mass * np.ones(P.shape[0]) + star_mass = star_mass * np.ones(P.shape[0]) - #defining the main coefficient of the mass equation to solve (comprised of the provided orbital parameter values) - D = ( C * P**(1/3) * K * np.sqrt(1 - e**2) ) / np.sin(I) + # defining the main coefficient of the mass equation to solve (comprised of the provided orbital parameter values) + D = (C * P**(1/3) * K * np.sqrt(1 - e**2)) / np.sin(Inc) m_ms = np.empty(P.shape, dtype=float) - #implementing a for loop to calculate the mass for the outer companions - #more accurately by including the mass of the inner companion(s) in the star_mass + # implementing a for loop to calculate the mass for the outer companions + # more accurately by including the mass of the inner companion(s) in the star_mass for comp in range(P.shape[1]): m_ms_comp = m_ms_func(D[:, comp], star_mass) m_ms[:, comp] = m_ms_comp star_mass = np.nansum([star_mass, m_ms_comp], axis=0) - #convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES + # convert to jupiter masses - NOTE, THIS IS INCONSISTENT WITH THE USE OF CONSTANTS FROM utils.py, AS DONE BELOW FOR THE CONVERSION TO EARTH MASSES m_mj = (m_ms * u.solMass).to(u.jupiterMass).value m_me = m_mj * mjup2mearth @@ -378,9 +396,9 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda else: return m_mj, np.nanmean(m_mj, axis=0), np.nanstd(m_mj, axis=0) - -def get_planet_mass_GAIA(P: Union[float, np.ndarray], a0: Union[float, np.ndarray], - parallax: Union[float, np.ndarray], + +def get_planet_mass_GAIA(P: Union[float, np.ndarray], a0: Union[float, np.ndarray], + parallax: Union[float, np.ndarray], star_mass: Union[float, Tuple] = 1.0, full_output=False): r""" Calculate the planet (minimum) mass, $M_p \sin i$, given orbital period `P`, @@ -485,7 +503,7 @@ def get_planet_mass_GAIA(P: Union[float, np.ndarray], a0: Union[float, np.ndarra def get_planet_semimajor_axis(P: Union[float, np.ndarray], K=None, - star_mass: Union[float, tuple]=1.0, + star_mass: Union[float, tuple] = 1.0, full_output=False): r""" Calculate the semi-major axis of the planet's orbit given orbital period @@ -568,8 +586,8 @@ def get_planet_semimajor_axis(P: Union[float, np.ndarray], K=None, def get_planet_semimajor_axis_accurate(P: Union[float, np.ndarray], M_c: Union[float, np.ndarray], - star_mass: Union[float, tuple]=1.0, - full_output=False): + star_mass: Union[float, tuple] = 1.0, + full_output=False): r""" Calculate the semi-major axis of the companion's orbit given orbital period `P`, companion mass `M_c`, and stellar mass. @@ -619,7 +637,6 @@ def get_planet_semimajor_axis_accurate(P: Union[float, np.ndarray], M_c: Union[f if isinstance(P, float): f = float(f) - try: # calculate for one value of the orbital period P = float(P) @@ -628,15 +645,19 @@ def get_planet_semimajor_axis_accurate(P: Union[float, np.ndarray], M_c: Union[f M_c = M_c * mjup2msun # convert to solar masses except TypeError: if isinstance(M_c, tuple) or isinstance(M_c, list): - M_c = (M_c[0] * mjup2msun, M_c[1] * mjup2msun) # convert to solar masses + # convert to solar masses + M_c = (M_c[0] * mjup2msun, M_c[1] * mjup2msun) else: - raise TypeError("M_c must be a scalar or tuple/list with shape (2,)") - + raise TypeError( + "M_c must be a scalar or tuple/list with shape (2,)") + if isinstance(star_mass, tuple) or isinstance(star_mass, list): - if isinstance(M_c, tuple) or isinstance(M_c, list): - a = f * (star_mass[0] + M_c[0])**(1/3) * (P / (2 * np.pi))**(2/3) + if isinstance(M_c, tuple) or isinstance(M_c, list): + a = f * (star_mass[0] + M_c[0])**(1/3) * \ + (P / (2 * np.pi))**(2/3) # error propagation: - a_err = np.sqrt(star_mass[1]**2 + M_c[1]**2) * a * (1/3) / (star_mass[0] + M_c[0]) + a_err = np.sqrt(star_mass[1]**2 + M_c[1]**2) * \ + a * (1/3) / (star_mass[0] + M_c[0]) return a, a_err else: a = f * (star_mass[0] + M_c)**(1/3) * (P / (2 * np.pi))**(2/3) @@ -663,10 +684,11 @@ def get_planet_semimajor_axis_accurate(P: Union[float, np.ndarray], M_c: Union[f a = np.empty(P.shape, dtype=float) M_c = M_c * mjup2msun # convert to solar masses - #accounting for the masses of the inner companion(s) when calculating - #the semi-major axis of the outer companion(s) + # accounting for the masses of the inner companion(s) when calculating + # the semi-major axis of the outer companion(s) for comp in range(P.shape[1]): - a_comp = f * (star_mass + M_c[:, comp])**(1/3) * (P[:, comp] / (2 * np.pi))**(2/3) + a_comp = f * (star_mass + M_c[:, comp])**(1/3) * \ + (P[:, comp] / (2 * np.pi))**(2/3) a[:, comp] = a_comp star_mass = np.nansum([star_mass, M_c[:, comp]], axis=0) @@ -675,6 +697,7 @@ def get_planet_semimajor_axis_accurate(P: Union[float, np.ndarray], M_c: Union[f else: return np.nanmean(a, axis=0), np.nanstd(a, axis=0) + def get_planet_mass_and_semimajor_axis(P, K, e, star_mass=1.0, full_output=False, verbose=False): """ @@ -712,8 +735,9 @@ def get_planet_mass_and_semimajor_axis(P, K, e, star_mass=1.0, return mass, a, star_mass return mass, a + def get_planet_mass_and_semimajor_axis_accurate(P, K, e, I, star_mass=1.0, - full_output=False, verbose=False): + full_output=False, verbose=False): """ Calculate the companion mass M and the semi-major axis given orbital period `P`, semi-amplitude `K`, eccentricity `e`, inclination `I`, and stellar mass. @@ -744,15 +768,16 @@ def get_planet_mass_and_semimajor_axis_accurate(P, K, e, I, star_mass=1.0, pass else: # include (Gaussian) uncertainty on the stellar mass - star_mass = star_mass_samples(*star_mass, P.shape[0]) + star_mass = star_mass_samples(*star_mass, P.shape[0]) mass = get_planet_mass_accurate(P, K, e, I, star_mass, full_output) - + a = get_planet_semimajor_axis_accurate(P, mass[0], star_mass, full_output) if isinstance(star_mass, np.ndarray): return mass, a, star_mass return mass, a + def get_bins(res, start=None, end=None, nbins=100): try: _start, _end = distribution_support(res.priors['Pprior']) @@ -769,18 +794,20 @@ def get_bins(res, start=None, end=None, nbins=100): start = _start if end is None: end = _end - + bins = 10**np.linspace(np.log10(start), np.log10(end), nbins) return bins + def aliases(P): - yearly, solar_dayly, sidereal_dayly, synodic_monthly = 1 / 365.25, 1.0, 1 + 1 / 365.25, 1 / 29 + yearly, solar_dayly, sidereal_dayly, synodic_monthly = 1 / \ + 365.25, 1.0, 1 + 1 / 365.25, 1 / 29 alias_year = [abs(1 / (yearly + i / P)) for i in [1, -1]] alias_solar_day = [abs(1 / (solar_dayly + i / P)) for i in [1, -1]] alias_sidereal_day = [abs(1 / (sidereal_dayly + i / P)) for i in [1, -1]] # alias_synodic_month = [abs(1 / (synodic_monthly + i / P)) for i in [1, -1]] - return alias_year, alias_solar_day, alias_sidereal_day #, alias_synodic_month + return alias_year, alias_solar_day, alias_sidereal_day # , alias_synodic_month # def FIP(results, plot=True, show_peaks=True, oversampling=1, include_aliases=False, # include_known_object=False, include_transiting_planet=False, @@ -790,11 +817,11 @@ def aliases(P): # Args: # res (kima.KimaResults): # The `KimaResults` instance -# plot (bool, optional): +# plot (bool, optional): # Plot the TIP and FIP. Defaults to True. -# show_peaks (bool, optional): +# show_peaks (bool, optional): # Identify and show prominent TIP peaks. Defaults to True. -# oversampling (int, optional): +# oversampling (int, optional): # Oversampling factor for the period binning. Defaults to 1. # Returns: @@ -841,10 +868,10 @@ def aliases(P): # # very memory intensive and slower: # # with np.errstate(divide='ignore'): # # c = np.logical_and( -# # # -# # a[:, None, None] < 1.0 / P, +# # # +# # a[:, None, None] < 1.0 / P, # # 1.0 / P < b[:, None, None] -# # # +# # # # # ).any(axis=2).sum(axis=1) # # TIP = c / results.ESS # # _1 = np.logical_and(a[:, None, None] < 1 / results.posteriors.P, 1 / results.posteriors.P < b[:, None, None]) @@ -878,9 +905,9 @@ def aliases(P): # peaks, _peaks_info = find_peaks(TIP, height=0.9) # for peak in peaks: # axs[0].plot(bins[peak], TIP[peak], 'ro', ms=3) -# axs[0].annotate(f'{bins[peak]:.2f}\nTIP: {TIP[peak]:.{ndigits}f}', -# xy=(bins[peak], TIP[peak]), xytext=(5, 0), -# xycoords='data', textcoords='offset points', +# axs[0].annotate(f'{bins[peak]:.2f}\nTIP: {TIP[peak]:.{ndigits}f}', +# xy=(bins[peak], TIP[peak]), xytext=(5, 0), +# xycoords='data', textcoords='offset points', # ha='left', va='top') # if show_ESS: # ess = f'{results.ESS:,}'.replace(',', r'$\,$') @@ -895,7 +922,7 @@ def aliases(P): # return bins, FIP -def FIP(results, f_center=None, f_width=None, oversampling=5, +def FIP(results, f_center=None, f_width=None, oversampling=5, plot=True, show_peaks=True, peak_indices=None, include_aliases=False, include_known_object=False, include_transiting_planet=False, peak_ndigits=None, show_ESS=True, just_tip=False, ax=None, **kwargs): @@ -914,10 +941,10 @@ def FIP(results, f_center=None, f_width=None, oversampling=5, Tobs = np.ptp(results.data.t) f_width = 1.0 / Tobs f_center = np.arange(f_width, 1.0, f_width / oversampling) - + if f_center is None: f_center = np.arange(f_width, 1.0, f_width / oversampling) - + P = results.posteriors.P.copy() if include_known_object: @@ -954,10 +981,12 @@ def FIP(results, f_center=None, f_width=None, oversampling=5, from scipy.signal import find_peaks if ax is None: if just_tip: - fig, axs = plt.subplots(1, 1, constrained_layout=True, figsize=(6, 3)) + fig, axs = plt.subplots( + 1, 1, constrained_layout=True, figsize=(6, 3)) axs = [axs] else: - fig, axs = plt.subplots(2, 1, constrained_layout=True, sharex=True) + fig, axs = plt.subplots( + 2, 1, constrained_layout=True, sharex=True) else: if just_tip: axs = [ax] @@ -980,9 +1009,9 @@ def FIP(results, f_center=None, f_width=None, oversampling=5, for peak in peaks: axs[0].plot(pcenter[peak], TIP[peak], 'ro', ms=3) - axs[0].annotate(f'{pcenter[peak]:.2f}\nTIP: {TIP[peak]:.{ndigits}f}', - xy=(pcenter[peak], TIP[peak]), xytext=(5, 0), - xycoords='data', textcoords='offset points', + axs[0].annotate(f'{pcenter[peak]:.2f}\nTIP: {TIP[peak]:.{ndigits}f}', + xy=(pcenter[peak], TIP[peak]), xytext=(5, 0), + xycoords='data', textcoords='offset points', ha='left', va='top') if show_ESS: ess = f'{results.ESS:,}'.replace(',', r'$\,$') @@ -990,7 +1019,8 @@ def FIP(results, f_center=None, f_width=None, oversampling=5, axs[0].set(xlabel='Period [days]', ylabel='TIP', ylim=(0, 1.05)) if not just_tip: - axs[1].set(xlabel='Period [days]', ylabel='-log$_{10}$ FIP', ylim=(0, None)) + axs[1].set(xlabel='Period [days]', + ylabel='-log$_{10}$ FIP', ylim=(0, None)) return pcenter, FIP, (fig, axs) @@ -1064,7 +1094,8 @@ def true_within_hdi(results, truths, hdi_prob=0.95, only_periods=False, P = truths[i][0] found = [P in interval for interval in intervals] if any(found): - print(f'Found {P=:<5.4f} d in one of the {hdi_prob} density regions') + print( + f'Found {P=:<5.4f} d in one of the {hdi_prob} density regions') nfound += 1 else: print(f'!! Did not find {P=:<5.4f} d in any density region') @@ -1118,9 +1149,9 @@ def true_within_hdi(results, truths, hdi_prob=0.95, only_periods=False, withins.append(within) - return {'found': nfound, 'hdi': regions, 'withins': withins} + def reorder_P(res, replace=False, passes=1): from .results import posterior_holder @@ -1141,11 +1172,16 @@ def reorder_P(res, replace=False, passes=1): P = new_posterior.P.copy() m, s = np.median(P[mask, 0]), np.std(P[mask, 0]) maybe_wrong_order = np.abs(P[:, 0] - m) > s - new_posterior.P[mask & maybe_wrong_order] = new_posterior.P[mask & maybe_wrong_order][:, ::-1] - new_posterior.K[mask & maybe_wrong_order] = new_posterior.K[mask & maybe_wrong_order][:, ::-1] - new_posterior.e[mask & maybe_wrong_order] = new_posterior.e[mask & maybe_wrong_order][:, ::-1] - new_posterior.w[mask & maybe_wrong_order] = new_posterior.w[mask & maybe_wrong_order][:, ::-1] - new_posterior.φ[mask & maybe_wrong_order] = new_posterior.φ[mask & maybe_wrong_order][:, ::-1] + new_posterior.P[mask & maybe_wrong_order] = new_posterior.P[mask & + maybe_wrong_order][:, ::-1] + new_posterior.K[mask & maybe_wrong_order] = new_posterior.K[mask & + maybe_wrong_order][:, ::-1] + new_posterior.e[mask & maybe_wrong_order] = new_posterior.e[mask & + maybe_wrong_order][:, ::-1] + new_posterior.w[mask & maybe_wrong_order] = new_posterior.w[mask & + maybe_wrong_order][:, ::-1] + new_posterior.φ[mask & maybe_wrong_order] = new_posterior.φ[mask & + maybe_wrong_order][:, ::-1] if replace: res.posteriors = new_posterior @@ -1257,6 +1293,8 @@ def reorder_P2(res, replace=False, passes=1): res.posterior_sample[:, res.indices['planets.e']] = new_posterior.e res.posterior_sample[:, res.indices['planets.w']] = new_posterior.w res.posterior_sample[:, res.indices['planets.φ']] = new_posterior.φ + + def reorder_P5(res, replace=False): from itertools import permutations from copy import deepcopy @@ -1318,7 +1356,7 @@ def sort_planet_samples(res, byP=True, replace=False): new_posterior.P = P[n, sortP].copy() for field in fields: try: - setattr(new_posterior, field, + setattr(new_posterior, field, getattr(res.posteriors, field)[n, sortP].copy()) except AttributeError: pass @@ -1340,7 +1378,8 @@ def sort_planet_samples(res, byP=True, replace=False): except AttributeError: pass try: - res.posterior_sample[:, res.indices['planets.' + field]] = getattr(new_posterior, field) + res.posterior_sample[:, res.indices['planets.' + + field]] = getattr(new_posterior, field) except (AttributeError, KeyError): pass @@ -1355,7 +1394,7 @@ def planet_parameters(results, star_mass=1.0, sample=None, printit=True): if isinstance(star_mass, (tuple, list)): mass_errs = True - format_tuples = lambda data: r'{:10.5f} \pm {:.5f}'.format( + def format_tuples(data): return r'{:10.5f} \pm {:.5f}'.format( data[0], data[1]) if mass_errs else '{:12.5f}'.format(data) if printit: @@ -1717,11 +1756,12 @@ def parameter_clusters(results, n_clusters=None, method='KMeans', if include_ecc: # data = res.posterior_sample[:, res.indices['planets']] - data = np.c_[res.posteriors.P.ravel(), res.posteriors.K.ravel(), res.posteriors.e.ravel()] + data = np.c_[res.posteriors.P.ravel(), res.posteriors.K.ravel(), + res.posteriors.e.ravel()] else: # data = np.c_[res.T, res.A, res.Omega, res.phi] data = np.c_[ - res.posteriors.P.ravel().copy(), + res.posteriors.P.ravel().copy(), res.posteriors.K.ravel().copy() ] @@ -2065,7 +2105,6 @@ def _columns_with_dynamic_range(results): return np.nonzero(dr)[0] - def hpd_grid(sample, alpha=0.05, bw_method=None, roundto=2): """Calculate highest posterior density (HPD) of array for given alpha. The HPD is the minimum width Bayesian credible interval (BCI). The function @@ -2096,7 +2135,7 @@ def hpd_grid(sample, alpha=0.05, bw_method=None, roundto=2): print('done') x = np.linspace(l, u, max(2000, sample.size//5)) y = density.evaluate(x) - #y = density.evaluate(x, l, u) waitting for PR to be accepted + # y = density.evaluate(x, l, u) waitting for PR to be accepted xy_zipped = zip(x, y/np.sum(y)) xy = sorted(xy_zipped, key=lambda x: x[1], reverse=True) xy_cum_sum = 0 @@ -2119,7 +2158,7 @@ def hpd_grid(sample, alpha=0.05, bw_method=None, roundto=2): hpd = list(zip(ite, ite)) modes = [] for value in hpd: - x_hpd = x[(x > value[0]) & (x < value[1])] - y_hpd = y[(x > value[0]) & (x < value[1])] - modes.append(round(x_hpd[np.argmax(y_hpd)], roundto)) - return hpd, x, y, modes \ No newline at end of file + x_hpd = x[(x > value[0]) & (x < value[1])] + y_hpd = y[(x > value[0]) & (x < value[1])] + modes.append(round(x_hpd[np.argmax(y_hpd)], roundto)) + return hpd, x, y, modes From e051bf9bdd1886cfc6fb40d70e731b69bfd900d3 Mon Sep 17 00:00:00 2001 From: cevawm Date: Fri, 6 Mar 2026 11:44:35 +0100 Subject: [PATCH 22/23] Add option to just get msini from new mass and sma functions --- src/kima/pykima/analysis.py | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/src/kima/pykima/analysis.py b/src/kima/pykima/analysis.py index 21bd7fd..257b053 100644 --- a/src/kima/pykima/analysis.py +++ b/src/kima/pykima/analysis.py @@ -230,7 +230,7 @@ def get_planet_mass(P: Union[float, np.ndarray], K: Union[float, np.ndarray], def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.ndarray], - e: Union[float, np.ndarray], Inc: Union[float, np.ndarray], + e: Union[float, np.ndarray], Inc: Union[float=None, np.ndarray=None], star_mass: Union[float, Tuple] = 1.0, full_output=False): r""" @@ -250,7 +250,7 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda e (Union[float, ndarray]): orbital eccentricity Inc (Union[float, ndarray]): - inclination [radians] + inclination [radians]. Defaults to None, in which case the minimum mass will be calculated. star_mass (Union[float, Tuple]): stellar mass, or (mass, uncertainty) [Msun] @@ -329,12 +329,17 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda # then K, e, and Inc should also be floats try: - K, e, Inc = float(K), float(e), float(Inc) + K, e = float(K), float(e) + if Inc is not None: + Inc = float(Inc) except TypeError: raise TypeError("K, e, and Inc should be floats if P is a float") # defining the main coefficient of the mass equation to solve (comprised of the provided orbital parameter values) - D = (C * P**(1/3) * K * np.sqrt(1 - e**2)) / np.sin(Inc) + if Inc is None: + D = (C * P**(1/3) * K * np.sqrt(1 - e**2)) + else: + D = (C * P**(1/3) * K * np.sqrt(1 - e**2)) / np.sin(Inc) if isinstance(star_mass, tuple) or isinstance(star_mass, list): @@ -376,7 +381,10 @@ def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.nda star_mass = star_mass * np.ones(P.shape[0]) # defining the main coefficient of the mass equation to solve (comprised of the provided orbital parameter values) - D = (C * P**(1/3) * K * np.sqrt(1 - e**2)) / np.sin(Inc) + if Inc is None: + D = (C * P**(1/3) * K * np.sqrt(1 - e**2)) + else: + D = (C * P**(1/3) * K * np.sqrt(1 - e**2)) / np.sin(Inc) m_ms = np.empty(P.shape, dtype=float) @@ -736,14 +744,15 @@ def get_planet_mass_and_semimajor_axis(P, K, e, star_mass=1.0, return mass, a -def get_planet_mass_and_semimajor_axis_accurate(P, K, e, I, star_mass=1.0, +def get_planet_mass_and_semimajor_axis_accurate(P, K, e, I=None, star_mass=1.0, full_output=False, verbose=False): """ Calculate the companion mass M and the semi-major axis given orbital period `P`, semi-amplitude `K`, eccentricity `e`, inclination `I`, and stellar mass. If `star_mass` is a tuple with (estimate, uncertainty), this (Gaussian) uncertainty will be taken into account in the calculation, and samples from - the star_mass distribution will be returned. + the star_mass distribution will be returned. If `I` is not provided, then the + minimum mass will be calculated, along with the corresponding semi-major axis. Units: P [days] From 7bdcf83aab383704bd1261b275c3c2aa838d9654 Mon Sep 17 00:00:00 2001 From: cevawm Date: Fri, 6 Mar 2026 13:20:14 +0100 Subject: [PATCH 23/23] Fix inclination argument bug --- src/kima/pykima/analysis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/kima/pykima/analysis.py b/src/kima/pykima/analysis.py index 257b053..fad27ff 100644 --- a/src/kima/pykima/analysis.py +++ b/src/kima/pykima/analysis.py @@ -230,7 +230,7 @@ def get_planet_mass(P: Union[float, np.ndarray], K: Union[float, np.ndarray], def get_planet_mass_accurate(P: Union[float, np.ndarray], K: Union[float, np.ndarray], - e: Union[float, np.ndarray], Inc: Union[float=None, np.ndarray=None], + e: Union[float, np.ndarray], Inc: Union[float, np.ndarray, None], star_mass: Union[float, Tuple] = 1.0, full_output=False): r"""