Skip to content

Commit

Permalink
Merge pull request #170 from sblunt/multiplanet_fix
Browse files Browse the repository at this point in the history 
improved mtot calculation for multiplanet
  • Loading branch information
@semaphoreP
semaphoreP committed May 14, 2020
2 parents 7d3292e + bf833f8 commit ac2c453
Show file tree
Hide file tree
Showing 2 changed files with 26 additions and 14 deletions.
19 changes: 12 additions & 7 deletions orbitize/system.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -261,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
Expand Down Expand Up @@ -293,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:
Expand Down
21 changes: 14 additions & 7 deletions tests/test_multiplanet.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down

0 comments on commit ac2c453

Please sign in to comment.