From 89cafc58b1158c297c77a3e90845fb045ab5f6e5 Mon Sep 17 00:00:00 2001 From: Jason Wang Date: Tue, 5 May 2020 15:40:04 -0700 Subject: [PATCH 1/2] improved mtot calculation for multiplanet --- orbitize/system.py | 9 +++++++-- tests/test_multiplanet.py | 21 ++++++++++++++------- 2 files changed, 21 insertions(+), 9 deletions(-) diff --git a/orbitize/system.py b/orbitize/system.py index 02844dc4..c24264e7 100644 --- a/orbitize/system.py +++ b/orbitize/system.py @@ -235,8 +235,13 @@ def compute_model(self, params_arr): # mass of secondary bodies are in order from -1-num_bodies until -2 in order. mass = params_arr[-1-self.num_secondary_bodies+(body_num-1)] m0 = params_arr[-1] - mtot = m0 + mass - # TODO: include the masses of other bodies? + # For what mtot to use to calculate central potential, we should use the mass enclosed in a sphere with r <= distance of planet. + # We need to select all planets with sma < this planet. + all_smas = params_arr[0:6*self.num_secondary_bodies:6] + within_orbit = np.where(all_smas <= sma) + all_pl_masses = params_arr[-1-self.num_secondary_bodies:-1] + inside_masses = all_pl_masses[within_orbit] + mtot = np.sum(inside_masses) + m0 else: # if not fitting for secondary mass, then total mass must be stellar mass mass = None diff --git a/tests/test_multiplanet.py b/tests/test_multiplanet.py index 4cbee512..a24dc75a 100644 --- a/tests/test_multiplanet.py +++ b/tests/test_multiplanet.py @@ -21,17 +21,20 @@ def test_compute_model(): mass_b = 0.001 # Msun m0 = 1 # Msun plx = 1 # mas - period_b = np.sqrt(b_params[0]**3/m0) # generate planet c orbital parameters # at 2 au, and starts on the opposite side of the star relative to b c_params = [2, 0, 0, np.pi, 0, 0] mass_c = 0.002 # Msun - period_c = np.sqrt(c_params[0]**3/m0) + + mtot = m0 + mass_b + mass_c + + period_c = np.sqrt(c_params[0]**3/mtot) + period_b = np.sqrt(b_params[0]**3/mtot) epochs = np.linspace(0, period_c*365.25, 100) + tau_ref_epoch # the full period of c, MJD - ra_model, dec_model, vz_model = kepler.calc_orbit(epochs, b_params[0], b_params[1], b_params[2], b_params[3], b_params[4], b_params[5], plx, m0+mass_b, tau_ref_epoch=tau_ref_epoch) + ra_model, dec_model, vz_model = kepler.calc_orbit(epochs, b_params[0], b_params[1], b_params[2], b_params[3], b_params[4], b_params[5], plx, mtot, tau_ref_epoch=tau_ref_epoch) # generate some fake measurements just to feed into system.py to test bookkeeping # just make a 1 planet fit for now @@ -78,21 +81,25 @@ def test_fit_selfconsist(): mass_b = 0.001 # Msun m0 = 1 # Msun plx = 1 # mas - period_b = np.sqrt(b_params[0]**3/m0) # generate planet c orbital parameters # at 2 au, and starts on the opposite side of the star relative to b c_params = [2, 0, 0, np.pi, 0, 0.5] mass_c = 0.002 # Msun - period_c = np.sqrt(c_params[0]**3/m0) + + mtot_c = m0 + mass_b + mass_c + mtot_b = m0 + mass_b + + period_b = np.sqrt(b_params[0]**3/mtot_b) + period_c = np.sqrt(c_params[0]**3/mtot_c) epochs = np.linspace(0, period_c*365.25, 20) + tau_ref_epoch # the full period of c, MJD # comptue Keplerian orbit of b - ra_model_b, dec_model_b, vz_model = kepler.calc_orbit(epochs, b_params[0], b_params[1], b_params[2], b_params[3], b_params[4], b_params[5], plx, m0+mass_b, mass_for_Kamp=m0, tau_ref_epoch=tau_ref_epoch) + ra_model_b, dec_model_b, vz_model = kepler.calc_orbit(epochs, b_params[0], b_params[1], b_params[2], b_params[3], b_params[4], b_params[5], plx, mtot_b, mass_for_Kamp=m0, tau_ref_epoch=tau_ref_epoch) # comptue Keplerian orbit of c - ra_model_c, dec_model_c, vz_model_c = kepler.calc_orbit(epochs, c_params[0], c_params[1], c_params[2], c_params[3], c_params[4], c_params[5], plx, m0+mass_c, tau_ref_epoch=tau_ref_epoch) + ra_model_c, dec_model_c, vz_model_c = kepler.calc_orbit(epochs, c_params[0], c_params[1], c_params[2], c_params[3], c_params[4], c_params[5], plx, mtot_c, tau_ref_epoch=tau_ref_epoch) # perturb b due to c ra_model_b_orig = np.copy(ra_model_b) From bf833f8b99744fed751b701d4731a74f7d96abae Mon Sep 17 00:00:00 2001 From: Jason Wang Date: Mon, 11 May 2020 15:56:42 -0700 Subject: [PATCH 2/2] m1/m0 -> m1/mtot --- orbitize/system.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/orbitize/system.py b/orbitize/system.py index c24264e7..0f05b92f 100644 --- a/orbitize/system.py +++ b/orbitize/system.py @@ -266,7 +266,7 @@ def compute_model(self, params_arr): # vz_i is the ith companion radial velocity if self.fit_secondary_mass: - vz0 = vz_i*-(mass/m0) # calculating stellar velocity due to ith companion + vz0 = vz_i*-(mass/mtot) # calculating stellar velocity due to ith companion total_rv0 = total_rv0 + vz0 # Adding stellar velocity and gamma # for the model points that correspond to this planet's orbit, add the model prediction @@ -298,11 +298,11 @@ def compute_model(self, params_arr): ## NOTE: we are only handling ra/dec and sep/pa right now ## TOOD: integrate RV into this if len(self.radec[other_body_num]) > 0: - radec_perturb[self.radec[other_body_num], 0] += -(mass/m0) * raoff[self.radec[other_body_num]] - radec_perturb[self.radec[other_body_num], 1] += -(mass/m0) * decoff[self.radec[other_body_num]] + radec_perturb[self.radec[other_body_num], 0] += -(mass/mtot) * raoff[self.radec[other_body_num]] + radec_perturb[self.radec[other_body_num], 1] += -(mass/mtot) * decoff[self.radec[other_body_num]] if len(self.seppa[other_body_num]) > 0: - radec_perturb[self.seppa[other_body_num], 0] += -(mass/m0) * raoff[self.seppa[other_body_num]] - radec_perturb[self.seppa[other_body_num], 1] += -(mass/m0) * decoff[self.seppa[other_body_num]] + radec_perturb[self.seppa[other_body_num], 0] += -(mass/mtot) * raoff[self.seppa[other_body_num]] + radec_perturb[self.seppa[other_body_num], 1] += -(mass/mtot) * decoff[self.seppa[other_body_num]] if self.fit_secondary_mass: if len(total_rv0[self.rv[0]]) > 0: