expanding `_get_consistent_inputs` for astrometric + rv parameter sets

[iancze]

For an orbit fit to both double-lined RV and astrometric data of a binary star, you basically fully orient everything in 3D space and determine the individual stellar masses. (There is no information about the radii of the stars). So, we'd like to choose a basic set of (independent) orbital parameters from which every other property can be derived.

Here are a few choices that might work

1) a (physical) period kappa = a1/(a1 + a2) -> Mtot is calculated from (a, period) -> M2 is calculated from Mtot * kappa -> M2 is calculated from Mtot - M2

Alternatively,

2) a (physical) period M2 -> Mtot is calculated from (a, period) -> M1 is calculated from Mtot - M2

or

3) a (physical) period q = M2/M1 -> Mtot is calculated from (a, period) -> M1 is calculated from Mtot and q -> M2 is calculated from Mtot and q

As currently written, if a, period, and m_star are given, _get_consistent_inputs blocks us.

I think option #2 parameterization is going to be the easiest to merge with the current code. So, if we give a, period, and m_planet = M2, then exoplanet calculates rho_star assuming r_star = 1.0.

Then, this flows down the control and m_star is calculated from rho_star and r_star, but, because r_star = 1.0, this isn't the right stellar mass implied by a, period, and m_planet.

# conversion constant
au_to_R_sun = (constants.au / constants.R_sun).value

a_ang = 0.324 # arcsec
parallax = 24.05 # milliarcsec
dpc = 1e3 / parallax
a = a_ang * dpc # au

e = 0.798
i = 96.0 * deg # [rad]
omega = 251.6 * deg - np.pi # Pourbaix reports omega_2, but we want omega_1
Omega = 159.6 * deg
P = 28.8 * 365.25 # days

T0 = Time(1989.92, format="decimalyear")
T0.format = "jd"
T0 = T0.value # [Julian Date]

gamma = 47.97 # km/s; systemic velocity

# kappa = a1 / (a1 + a2)
kappa = 0.45

# calculate Mtot from a, P
Mtot = calc_Mtot(a, P)

M2 = kappa * Mtot
M1 = Mtot - M2

# calculate M1, M2 from kappa and Mtot
print("Masses from Kepler's third law M1:{:.2f}, M2:{:.2f}, Mtot:{:.2f}".format(M1, M2, Mtot))

# n.b. that we include an extra conversion for a, because exoplanet expects a in R_sun
orbit = xo.orbits.KeplerianOrbit(a=a * au_to_R_sun, t_periastron=T0, period=P,
                               incl=i, ecc=e, omega=omega, Omega=Omega, m_planet=M2)

print("m_planet: {:.2f}".format(orbit.m_planet.eval()))
print("m_star: {:.2f}".format(orbit.m_star.eval()))
print("m_tot: {:.2f}".format(orbit.m_total.eval()))

# Output
Correct masses from Kepler's third law M1:1.62, M2:1.33, Mtot:2.95
Calculated masses from _get_consistent_inputs
m_planet: 1.33
m_star: 0.27
m_tot: 1.60

I'm not super familiar with the minimum parameter sets required for transit fitting, or transit + RV fitting, so I'm a little hesitant to dive in and start disrupting control flow through this routine. I would be inclined to change things such that the user specifies a, P, and m_planet and instead of calculating rho_star directly, _get_consistent_inputs calculates m_total and m_star = m_total - m_planet. I could use a little guidance on whether this is advisable package-wide, however.

[dfm]

I also think that I would prefer the second parameter set. The others can easily be calculated in user code.

I'm a bit confused about your code snippet. There could well be a thinko in the way I currently have it implemented. What is the problem that is causing this?

[iancze]

I think it's because of the specification that r_star = 1.0.

Then, this flows down the control and m_star is calculated from rho_star and r_star, but, because r_star = 1.0, this isn't the right stellar mass implied by a, period, and m_planet.

[dfm]

Sure - I saw that. But what should it be instead?

[dfm]

I guess I'm asking: what is the calc_Mtot function above? In principle the current implementation of this function should do what you want in principle and it would be good to work out the bug if there is one.

[dfm]

I think I figured it out and now I understand your suggestion that we use m_total as the parameter. I think that's a good idea, but we'll have to think through the options.

Can you make a list of example test cases like the one above? I'll go through and do the same for transits and RVs. That way we get a reasonable set of unit tests for this function!

[iancze]

Sorry for the delay in getting back (conference talks...). The calc_Mtot function is just calculating Mtot from Kepler's third law. As you mentioned, I think this is also done via _get_consistent_inputs, but since it's currently routed through the stellar density, the r_star issue comes into play.

For the astrometric + RV, I think those three are main test cases I have seen, but option two (a (physical), period, and m_planet) has the best ability to eventually scale to multiple planets. For example, m_planet can still be a vector of multiple planets and m_star is just calculated as m_tot - tt.sum(m_planet), where m_tot is calculated from a and period. So, as long as we can eventually get back out the correct m_star (to use in the RV solution) I think this is a good parameter set to stick with. I'm not even sure it makes sense to start accepting either option 1 or option 3 choices, since these immediately break when you have more than two planets.

One aspect of this whole process that might benefit from a convenience function is converting the positions of the planet and star back into arcseconds. As I've written things, it seems a little cyclical to model using a parameter set of a_ang and parallax, convert to a_phys (in R_sun) to feed to KeplerianOrbit, get_position on the sky plane in R_sun, and then use parallax to convert back to arcseconds. E.g., I'm doing something like this

parallax = pm.Normal("parallax", mu=27.31, sd=0.12) # milliarcsec 
a_ang = pm.Uniform("a_ang", 0.1, 7.0, testval=3.0) # arcsec 

# the distance to the source in parcsecs
dpc = pm.Deterministic("dpc", 1e3 / parallax)

# the semi-major axis in au
a = pm.Deterministic("a", a_ang * dpc) # au

# n.b. that we include an extra conversion for a, because exoplanet expects a in R_sun
orbit = xo.orbits.KeplerianOrbit(a=a * au_to_R_sun, t_periastron=t_periastron, period=P, 
                               incl=incl, ecc=e, omega=omega, Omega=Omega)

rho_phys, theta_model = orbit.get_relative_angles(jds) # the rho, theta model values

# because we've specified a physical value for a, a is now actually in units of R_sun
# So, we'll want to convert back to arcsecs 
rho_model = (rho_phys / au_to_R_sun) / dpc # arcsec

One solution might be to add a parallax parameter to get_position, if it's not None, then the scale is converted from R_sun to arcseconds using the value of parallax.

For the astrometric only case, I think the (a_phys, P) parameterizations work well enough and I'm not sure what else you could fit with instead. Doing this without a parallax (i.e, using a_ang in arcsec as a_phys) works but is a little bit sketchy since it's fooling get_position into doing the above convenience transformation but in a hard-coded and somewhat hidden fashion. What would probably be better in terms of interpretability is to just choose a fake (fixed) parallax, use that to convert a_ang to a_phys to feed to __init__ and then use the same fake parallax to convert get_position back to arcsecs.

[dfm]

Hey @iancze: I just pushed a change to get_consistent_inputs so that it now passes your test. The issue was actually way stupider than we thought... I was just computing the implied density in the wrong units!! That is fixed now so take a look and let me know what you're thinking about parallaxes.

[iancze]

This update works for me! I added the parallax keywords and things seem to be working out. Putting the finishing touches on the astrometric + RV tutorial and I'll open a pull request imminently.