Constraints

Setup

Let's first make sure we have the latest version of PHOEBE 2.0 installed. (You can comment out this line if you don't use pip for your installation or don't want to update to the latest release).


In [ ]:
!pip install -I "phoebe>=2.0,<2.1"

As always, let's do imports and initialize a logger and a new Bundle. See Building a System for more details.


In [1]:
import phoebe
from phoebe import u # units
import numpy as np
import matplotlib.pyplot as plt

logger = phoebe.logger()

b = phoebe.default_binary()


WARNING: Constant u'Gravitational constant' is already has a definition in the u'si' system [astropy.constants.constant]
WARNING: Constant u'Solar mass' is already has a definition in the u'si' system [astropy.constants.constant]
WARNING: Constant u'Solar radius' is already has a definition in the u'si' system [astropy.constants.constant]
WARNING: Constant u'Solar luminosity' is already has a definition in the u'si' system [astropy.constants.constant]
/usr/local/lib/python2.7/dist-packages/astropy/units/quantity.py:782: FutureWarning: comparison to `None` will result in an elementwise object comparison in the future.
  return super(Quantity, self).__eq__(other)

What are Constraints?

Constraints live in their own context of the Bundle, and many are created by default - either when you add a component or when you set the system hierarchy.

Let's look at all the existing constraints for our binary system.


In [2]:
b.filter(context='constraint')


Out[2]:
<ParameterSet: 20 parameters | kinds: star, orbit>

To see what all of these constraints do, see the 'Built-in Constraints' section below.

For now let's look at a single constraint.


In [3]:
b['constraint']['primary']['mass']


Out[3]:
<ConstraintParameter: {mass@primary@component} = (39.478418 * ({sma@binary@component} ** 3.000000)) / ((({period@binary@component} ** 2.000000) * (1.000000 + {q@binary@component})) * 0.000000000066740799999999990427) => 0.998813135806 solMass>

Here we see the equation used to derive the mass of the primary star from its orbit, as well as the current value

If we access the Parameter that it is constraining we can see that it is automatically kept up-to-date.


In [4]:
print b.get_value('mass@primary@component')


0.998813135806

The parameter is aware that it's being constrained and all the relevant linking parameters.


In [5]:
print b['mass@primary@component']


Parameter: mass@primary@component
                       Qualifier: mass
                     Description: Mass
                           Value: 0.998813135806 solMass
                  Constrained by: sma@binary@component, period@binary@component, q@binary@component
                      Constrains: None
                      Related to: sma@binary@component, period@binary@component, q@binary@component

If you change the hierarchy, built-in cross-object constraints (like mass that depends on its parent orbit), will be adjusted to reflect the new hierarchy. See the 'Changing Hierarchies' section below for more details.

Built-in Constraints

There are a number of built-in constraints that will be applied to your system by default. These are all listed below:

asini

This constraint handles computing the projected semi-major axis along the line of sight and can be automatically inverted to solve for either 'asini', 'sma', or 'incl'.


In [6]:
b['asini@constraint']


Out[6]:
<ConstraintParameter: {asini@binary@component} = {sma@binary@component} * (sin({incl@binary@component})) => 5.3 solRad>

esinw, ecosw

These constraints handle computing the projected eccentricity which can be helpful in that they are better representations of the geometry of a light curve and result in symmetric posteriors for near-circular orbits.

Both can be inverted to also automatically solve for 'ecc' or 'per0'.


In [7]:
b['esinw@constraint']


Out[7]:
<ConstraintParameter: {esinw@binary@component} = {ecc@binary@component} * (sin({per0@binary@component})) => 0.0>

In [8]:
b['ecosw@constraint']


Out[8]:
<ConstraintParameter: {ecosw@binary@component} = {ecc@binary@component} * (cos({per0@binary@component})) => 0.0>

t0

This constraint handles converting between different t0 conventions - namely providing a reference time at periastron passage (t0_perpass) and at superior conjunction (t0_supconj).

Currently, this constraint only supports inverting to be solved for 't0_supconj' (ie you cannot automatically invert this constraint to constraint phshift or per0).


In [9]:
b['t0_perpass@constraint']


Out[9]:
<ConstraintParameter: {t0_perpass@binary@component} = {t0_supconj@binary@component} + ((({phshift@binary@component} - 0.250000) + ({per0@binary@component} / 6.283185307179586231995926937088)) * {period@binary@component}) => -0.25 d>

freq

This constraint handles the simple conversion to frequency from period - whether that be rotational or orbital - and does support inversion to solve for 'period'.


In [10]:
b['freq@constraint']


Out[10]:
<ParameterSet: 3 parameters | kinds: star, orbit>

In [11]:
b['freq@binary@constraint']


Out[11]:
<ConstraintParameter: {freq@binary@component} = 6.283185 / {period@binary@component} => 6.283185 rad / d>

In [12]:
b['freq@primary@constraint']


Out[12]:
<ConstraintParameter: {freq@primary@component} = 6.283185 / {period@primary@component} => 6.283185 rad / d>

mass

This constraint handles solving for the mass of a component by obeying Kepler's third law within the parent orbit.

It can be inverted to solve for 'sma' or 'period' (in addition to 'mass'), but not 'q'.


In [13]:
b['mass@constraint']


Out[13]:
<ParameterSet: 2 parameters | components: primary, secondary>

In [14]:
b['mass@primary@constraint']


Out[14]:
<ConstraintParameter: {mass@primary@component} = (39.478418 * ({sma@binary@component} ** 3.000000)) / ((({period@binary@component} ** 2.000000) * (1.000000 + {q@binary@component})) * 0.000000000066740799999999990427) => 0.998813135806 solMass>

component sma

This constraint handles computing the semi-major axis of a component about the center of mass of its parent orbit. Note that this is not the same as the semi-major axis of the parent orbit.

This currently can be inverted to solve for 'sma' of the parent orbit, but not 'q'.


In [15]:
b['sma@constraint']


Out[15]:
<ParameterSet: 2 parameters | components: primary, secondary>

In [16]:
b['sma@primary@constraint']


Out[16]:
<ConstraintParameter: {sma@primary@component} = {sma@binary@component} / (1.000000 + (1.000000 / {q@binary@component})) => 2.65 solRad>

potential

This constraint handles solving for the equipotential of the roche surface from the polar radius.

This currently can be inverted to solve for 'rpole' but not 'sma', 'syncpar', 'ecc', or 'q'.


In [17]:
b['pot@constraint']


Out[17]:
<ParameterSet: 2 parameters | components: primary, secondary>

In [18]:
b['pot@primary@constraint']


Out[18]:
<ConstraintParameter: {pot@primary@component} = rocherpole2potential({rpole@primary@component}, {q@binary@component}, {ecc@binary@component}, {syncpar@primary@component}, {sma@binary@component}, 1) => 6.28266165375>

rotation period

This constraint handles computing the rotation period of a star given its synchronicity parameter (syncpar).

It can be inverted to solve for any of the three parameters 'period' (both rotational and orbital) and 'syncpar'.


In [19]:
b['period@constraint']


Out[19]:
<ParameterSet: 2 parameters | components: primary, secondary>

In [20]:
b['period@primary@constraint']


Out[20]:
<ConstraintParameter: {period@primary@component} = {period@binary@component} / {syncpar@primary@component} => 1.0 d>

inclination (aligned binary)

This constraint handles requiring the inclination of a star to be the same as that of the parent orbit. With the constraint enabled, the binary is required to be "aligned". Disabling this constraint (which currently is not yet supported) will allow for misaligned binaries.

This constraint can be inverted to solve for either 'incl'.


In [21]:
b['incl@constraint']


Out[21]:
<ParameterSet: 2 parameters | components: primary, secondary>

In [22]:
b['incl@primary@constraint']


Out[22]:
<ConstraintParameter: {incl@primary@component} = {incl@binary@component} => 90.0 deg>

Re-Parameterizing

NOTE: this is an experimental feature. When re-parameterizing, please be careful and make sure all results and parameters make sense.

As we've just seen, the mass is a constrained (ie derived) parameter. But let's say that you would rather provide masses for some reason (perhaps that was what was provided in a paper). You can choose to provide mass and instead have one of its related parameters constrained


In [23]:
print b['mass@primary@component'].constrained_by


[<Parameter: sma=5.3 solRad | keys: description, value, quantity, default_unit, limits, visible_if, copy_for>, <Parameter: period=1.0 d | keys: description, value, quantity, default_unit, limits, visible_if, copy_for>, <Parameter: q=1.0 | keys: description, value, quantity, default_unit, limits, visible_if, copy_for>]

In [24]:
print b['value@mass@primary@component'], b['value@mass@secondary@component'], b['value@period@orbit@component']


0.998813135806 0.998813135806 1.0

In [25]:
b.flip_constraint('mass@primary', 'period')


Out[25]:
<ConstraintParameter: {period@binary@component} = ((39.478418 * ({sma@binary@component} ** 3.000000)) / (({mass@primary@component} * (1.000000 + {q@binary@component})) * 0.000000000066740799999999990427)) ** 0.500000 => 1.0 d>

In [26]:
b['mass@primary@component'] = 1.0

In [27]:
print b['value@mass@primary@component'], b['value@mass@secondary@component'], b['value@period@orbit@component']


1.0 1.0 0.999406391718

You'll see that when we set the primary mass, the secondary mass has also changed (because the masses are related through q) and the period has changed (based on resolving the Kepler's third law constraint).

Note that the tags for the constraint are based on those of the constrained parameter, so to switch the parameterization back, we'll have to use a slightly different twig.


In [28]:
print b['constraint']


ParameterSet: 20 parameters
         freq@primary@constraint: 6.283185 / {period@primary@component}
  irrad_frac_lost_bol@primary...: 1.000000 - {irrad_frac_refl_bol@primary@component}
       freq@secondary@constraint: 6.283185 / {period@secondary@component}
  irrad_frac_lost_bol@seconda...: 1.000000 - {irrad_frac_refl_bol@secondary@component}
                asini@constraint: {sma@binary@component} * (sin({incl@binary@component}))
           t0_perpass@constraint: {t0_supconj@binary@component} + ((({phshift@binary@component} - 0.250000) + ({per0@binary@component} / 6.283185307179586231995926937088)) * {period@binary@component})
            mean_anom@constraint: (6.283185 * ({t0@system} - {t0_perpass@binary@component})) / {period@binary@component}
                ecosw@constraint: {ecc@binary@component} * (cos({per0@binary@component}))
                esinw@constraint: {ecc@binary@component} * (sin({per0@binary@component}))
          freq@binary@constraint: 6.283185 / {period@binary@component}
        period@binary@constraint: ((39.478418 * ({sma@binary@component} ** 3.000000)) / (({mass@primary@component} * (1.000000 + {q@binary@component})) * 0.000000000066740799999999990427)) ** 0.500000
          sma@primary@constraint: {sma@binary@component} / (1.000000 + (1.000000 / {q@binary@component}))
       period@primary@constraint: {period@binary@component} / {syncpar@primary@component}
         incl@primary@constraint: {incl@binary@component}
          pot@primary@constraint: rocherpole2potential({rpole@primary@component}, {q@binary@component}, {ecc@binary@component}, {syncpar@primary@component}, {sma@binary@component}, 1)
                 mass@constraint: (39.478418 * ({sma@binary@component} ** 3.000000)) / ((({period@binary@component} ** 2.000000) * (1.000000 + (1.000000 / {q@binary@component}))) * 0.000000000066740799999999990427)
        sma@secondary@constraint: {sma@binary@component} / (1.000000 + {q@binary@component})
     period@secondary@constraint: {period@binary@component} / {syncpar@secondary@component}
       incl@secondary@constraint: {incl@binary@component}
        pot@secondary@constraint: rocherpole2potential({rpole@secondary@component}, {q@binary@component}, {ecc@binary@component}, {syncpar@secondary@component}, {sma@binary@component}, 2)

In [29]:
b['period@constraint@binary']


Out[29]:
<ConstraintParameter: {period@binary@component} = ((39.478418 * ({sma@binary@component} ** 3.000000)) / (({mass@primary@component} * (1.000000 + {q@binary@component})) * 0.000000000066740799999999990427)) ** 0.500000 => 0.999406391718 d>

In [30]:
b['period@constraint@binary'].meta


Out[30]:
OrderedDict([('time', None),
             ('qualifier', 'period'),
             ('history', None),
             ('feature', None),
             ('component', 'binary'),
             ('dataset', None),
             ('constraint', None),
             ('compute', None),
             ('model', None),
             ('fitting', None),
             ('feedback', None),
             ('plugin', None),
             ('kind', 'orbit'),
             ('context', 'constraint'),
             ('twig', 'period@binary@orbit@constraint'),
             ('uniquetwig', 'period@binary@constraint')])

Notice that the qualifier tag has changed from 'mass' to 'period' and the 'component' tag has changed from 'primary' to 'binary' (since sma is in the binary).


In [31]:
b.flip_constraint('period@binary', 'mass')


Out[31]:
<ConstraintParameter: {mass@primary@component} = (39.478418 * ({sma@binary@component} ** 3.000000)) / ((({period@binary@component} ** 2.000000) * (1.000000 + {q@binary@component})) * 0.000000000066740799999999990427) => 1.0 solMass>

Changing Hierarchies

NOTE: this is an experimental feature. Just as when re-parameterizing, please be careful and make sure all results and parameters make sense.

Some of the built-in constraints depend on the system hierarchy, and will automatically adjust to reflect changes to the hierarchy.

For example, the masses depend on the period and semi-major axis of the parent orbit but also depend on the mass-ratio (q) which is defined as the primary mass over secondary mass. For this reason, changing the roles of the primary and secondary components should be reflected in the masses (so long as q remains fixed).

In order to show this example, let's set the mass-ratio to be non-unity.


In [32]:
b.set_value('q', 0.8)

Here the star with component tag 'primary' is actually the primary component in the hierarchy, so should have the LARGER mass (for a q < 1.0).


In [33]:
print "M1: {}, M2: {}".format(b.get_value('mass@primary@component'),
                              b.get_value('mass@secondary@component'))


M1: 1.11111111111, M2: 0.888888888889

Now let's flip the hierarchy so that the star with the 'primary' component tag is actually the secondary component in the system (and so takes the role of numerator in q = M2/M1).

For more information on the syntax for setting hierarchies, see the Building a System Tutorial.


In [34]:
b.set_hierarchy('orbit:binary(star:secondary, star:primary)')

In [35]:
print b.get_value('q')


0.8

In [36]:
print "M1: {}, M2: {}".format(b.get_value('mass@primary@component'),
                              b.get_value('mass@secondary@component'))


M1: 0.888888888889, M2: 1.11111111111

Even though under-the-hood the constraints are being rebuilt from scratch, they will remember if you have flipped them to solve for some other parameter.

To show this, let's flip the constraint for the secondary mass to solve for 'period' and then change the hierarchy back to its original value.


In [37]:
print "M1: {}, M2: {}, period: {}, q: {}".format(b.get_value('mass@primary@component'),
                                                 b.get_value('mass@secondary@component'),
                                                 b.get_value('period@binary@component'),
                                                 b.get_value('q@binary@component'))


M1: 0.888888888889, M2: 1.11111111111, period: 0.999406391718, q: 0.8

In [38]:
b.flip_constraint('mass@secondary@constraint', 'period')


Out[38]:
<ConstraintParameter: {period@binary@component} = ((39.478418 * ({sma@binary@component} ** 3.000000)) / (({mass@secondary@component} * (1.000000 + {q@binary@component})) * 0.000000000066740799999999990427)) ** 0.500000 => 0.999406391718 d>

In [39]:
print "M1: {}, M2: {}, period: {}, q: {}".format(b.get_value('mass@primary@component'),
                                                 b.get_value('mass@secondary@component'),
                                                 b.get_value('period@binary@component'),
                                                 b.get_value('q@binary@component'))


M1: 0.888888888889, M2: 1.11111111111, period: 0.999406391718, q: 0.8

In [40]:
b.set_value('mass@secondary@component', 1.0)

In [41]:
print "M1: {}, M2: {}, period: {}, q: {}".format(b.get_value('mass@primary@component'),
                                                 b.get_value('mass@secondary@component'),
                                                 b.get_value('period@binary@component'),
                                                 b.get_value('q@binary@component'))


M1: 0.8, M2: 1.0, period: 1.05346683532, q: 0.8

Next

Next up: let's add a dataset to our Bundle.


In [ ]: