Potentials

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]:
%matplotlib inline

In [2]:
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)

Now let's add a mesh dataset at a few different times so that we can see how the potentials affect the surfaces of the stars.


In [3]:
b.add_dataset('mesh', times=np.linspace(0,1,11), dataset='mesh01')


Out[3]:
<ParameterSet: 2 parameters | contexts: compute, dataset>

Relevant Parameters

The 'pot' parameter defines the stellar surface to be at the given Roche equipotential - and is given at periastron. By default, the 'pot' parameter is constrained (in contrast to PHOEBE legacy and Wilson-Devinney) and instead the 'rpole' parameter should be used to change the size of a star.


In [4]:
print b['pot@component']


ParameterSet: 2 parameters
*          pot@primary@component: 6.28266165375
*        pot@secondary@component: 6.28266165375

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


Parameter: pot@primary@component
                       Qualifier: pot
                     Description: Potential at periastron
                           Value: 6.28266165375
                  Constrained by: rpole@primary@component, q@binary@component, ecc@binary@component, syncpar@primary@component, sma@binary@component
                      Constrains: None
                      Related to: rpole@primary@component, q@binary@component, ecc@binary@component, syncpar@primary@component, sma@binary@component
                 Only visible if: hierarchy.is_contact_binary:False


In [6]:
print b['pot@primary@constraint']


Constrains (qualifier): pot
Expression in SI (value): rocherpole2potential({rpole@primary@component}, {q@binary@component}, {ecc@binary@component}, {syncpar@primary@component}, {sma@binary@component}, 1)
Current Result (result): 6.28266165375

The 'rpole' parameter is defined as the polar radius of the star at periastron.


In [7]:
print b['rpole']


ParameterSet: 2 parameters
         rpole@primary@component: 1.0 solRad
       rpole@secondary@component: 1.0 solRad

In [8]:
print b['rpole@primary']


Parameter: rpole@primary@component
                       Qualifier: rpole
                     Description: Polar radius at periastron
                           Value: 1.0 solRad
                  Constrained by: 
                      Constrains: pot@primary@component
                      Related to: q@binary@component, ecc@binary@component, syncpar@primary@component, sma@binary@component, pot@primary@component
                 Only visible if: hierarchy.is_contact_binary:False

Critical Potentials and System Checks

If you set a value such that the system becomes non-computable, a logger warning will immediately be raised. This will happen in a detached system, for example, if any of the stars are predicted to overflow at periastron. Since the surface of the star (potential) depends on many parameters (see the constraint above), this can be triggered by changing any of:

  • rpole
  • pot
  • q
  • sma
  • ecc
  • syncpar

In [9]:
b.set_value('rpole@primary@component', 3)


Fri, 10 Feb 2017 15:43 PARAMETERS   WARNING primary is overflowing at periastron (L1=3.75, L2=3.21)
Fri, 10 Feb 2017 15:43 PARAMETERS   WARNING primary is overflowing at periastron (L1=3.75, L2=3.21)

At this time, if you were to call run_compute, an error would be thrown. An error isn't immediately thrown when setting rpole, however, since the overflow can be recitified by changing any of the other relevant parameters. For instance, let's change sma to be large enough to account for this value of rpole and you'll see that the error does not occur again.


In [10]:
b.set_value('sma@binary@component', 5)


Fri, 10 Feb 2017 15:43 PARAMETERS   WARNING primary is overflowing at periastron (L1=3.75, L2=3.21)
Fri, 10 Feb 2017 15:43 PARAMETERS   WARNING primary is overflowing at periastron (L1=3.75, L2=3.21)
Fri, 10 Feb 2017 15:44 PARAMETERS   WARNING primary is overflowing at periastron (L1=3.75, L2=3.21)
Fri, 10 Feb 2017 15:44 PARAMETERS   WARNING primary is overflowing at periastron (L1=3.75, L2=3.21)
Fri, 10 Feb 2017 15:44 PARAMETERS   WARNING primary is overflowing at periastron (L1=3.75, L2=3.21)
Fri, 10 Feb 2017 15:44 PARAMETERS   WARNING primary is overflowing at periastron (L1=3.75, L2=3.21)
Fri, 10 Feb 2017 15:44 PARAMETERS   WARNING primary is overflowing at periastron (L1=3.75, L2=3.21)
Fri, 10 Feb 2017 15:44 PARAMETERS   WARNING primary is overflowing at periastron (L1=3.75, L2=3.21)

In [11]:
b.set_value('sma@binary@component', 10)


Fri, 10 Feb 2017 15:44 PARAMETERS   WARNING primary is overflowing at periastron (L1=3.75, L2=3.21)
Fri, 10 Feb 2017 15:44 PARAMETERS   WARNING primary is overflowing at periastron (L1=3.75, L2=3.21)
Fri, 10 Feb 2017 15:44 PARAMETERS   WARNING primary is overflowing at periastron (L1=3.75, L2=3.21)

These logger warnings are handy when running phoebe interactively, but in a script its also handy to be able to check whether the system is currently computable /before/ running run_compute.

This can be done by calling run_checks which returns a boolean (whether the system passes all checks) and a message (a string describing the first failed check).


In [12]:
passed, message = b.run_checks()
print passed, message


True 

In [13]:
b.set_value('sma@binary@component', 5)


Fri, 10 Feb 2017 15:44 PARAMETERS   WARNING primary is overflowing at periastron (L1=3.75, L2=3.21)
Fri, 10 Feb 2017 15:44 PARAMETERS   WARNING primary is overflowing at periastron (L1=3.75, L2=3.21)
Fri, 10 Feb 2017 15:44 PARAMETERS   WARNING primary is overflowing at periastron (L1=3.75, L2=3.21)
Fri, 10 Feb 2017 15:44 PARAMETERS   WARNING primary is overflowing at periastron (L1=3.75, L2=3.21)
Fri, 10 Feb 2017 15:44 PARAMETERS   WARNING primary is overflowing at periastron (L1=3.75, L2=3.21)

In [14]:
passed, message = b.run_checks()
print passed, message


False primary is overflowing at periastron (L1=3.75, L2=3.21)

Alternatively, you could also manually check with the critical potentials.

The set value (or constrained value) of the equipotential (which is defined at periastron) must be GREATER than the critical potentials at each of the lagrange points in order to not overflow.


In [15]:
q = b.get_value('q@binary@component')
F = b.get_value('syncpar@primary@component')
# instantaneous distance in units of sma at periastron is 1-e
d = 1 - b.get_value('ecc@binary@component')
crit_pots = phoebe.libphoebe.roche_critical_potential(q, F, d, L1=True, L2=True, L3=True)
pot = b.get_value('pot@primary@component')
print "pot: {}\nL1_crit: {}\nL2_crit: {}\nL3_crit: {}".format(pot, crit_pots['L1'], crit_pots['L2'], crit_pots['L3'])


pot: 2.52415959238
L1_crit: 3.75
L2_crit: 3.20679622409
L3_crit: 3.20679622409

Here we can see that our star is overflowing at L2. By changing the sma to a value that passed before, we see that the potential changes to a value that is allowed.


In [16]:
b.set_value('sma@binary@component', 10)


Fri, 10 Feb 2017 15:44 PARAMETERS   WARNING primary is overflowing at periastron (L1=3.75, L2=3.21)
Fri, 10 Feb 2017 15:44 PARAMETERS   WARNING primary is overflowing at periastron (L1=3.75, L2=3.21)
Fri, 10 Feb 2017 15:44 PARAMETERS   WARNING primary is overflowing at periastron (L1=3.75, L2=3.21)

In [17]:
q = b.get_value('q@binary@component')
F = b.get_value('syncpar@primary@component')
# instantaneous distance in units of sma at periastron is 1-e
d = 1 - b.get_value('ecc@binary@component')
crit_pots = phoebe.libphoebe.roche_critical_potential(q, F, d, L1=True, L2=True, L3=True)
pot = b.get_value('pot@primary@component')
print "pot: {}\nL1_crit: {}\nL2_crit: {}\nL3_crit: {}".format(pot, crit_pots['L1'], crit_pots['L2'], crit_pots['L3'])


pot: 4.29115961855
L1_crit: 3.75
L2_crit: 3.20679622409
L3_crit: 3.20679622409

Note that when calling critical potentials for the secondary component, you must invert q.


In [18]:
q = 1./b.get_value('q@binary@component')
F = b.get_value('syncpar@secondary@component')
# instantaneous distance in units of sma at periastron is 1-e
d = 1 - b.get_value('ecc@binary@component')
crit_pots = phoebe.libphoebe.roche_critical_potential(q, F, d, L1=True, L2=True, L3=True)
pot = b.get_value('pot@primary@component')
print "pot: {}\nL1_crit: {}\nL2_crit: {}\nL3_crit: {}".format(pot, crit_pots['L1'], crit_pots['L2'], crit_pots['L3'])


pot: 4.29115961855
L1_crit: 3.75
L2_crit: 3.20679622409
L3_crit: 3.20679622409

Semi-Detached Systems

NOTE: Support for semi-detached systems is introduced in version 2.0.5

Semi-detached systems are implemented by constraining either the potential or polar radius to be at the critical L1 value.

Here since the potential is already constrained by rpole, we'll constrain rpole. This is done by applying the 'critical_rpole' constraint on the 'primary' component.


In [23]:
b.add_constraint('critical_rpole', 'primary')


Out[23]:
<ParameterSet: 1 parameters>

If instead rpole was constrained by pot, then we would use the 'critical_potential' constraint instead.

We can view the constraint on rpole by accessing the constraint:


In [26]:
b['rpole@constraint']


Out[26]:
<ConstraintParameter: {rpole@primary@component} = rochecriticalL12rpole({q@binary@component}, {ecc@binary@component}, {syncpar@primary@component}, {sma@binary@component}, 1) => 3.56130872221 solRad>

Now whenever any of the relevant parameters (q, ecc, syncpar, sma) are changed, the rpole and pot will automatically adjust so that they are at the critical L1 value.

Accessing Synthetic Potentials and Polar Radiii

Potentials and Polar Radii are stored in the synthetic meshes for every time in which a mesh is computed. For circular orbits, these should remain constant over time and be equivalent to the input parameters. In the next section, we'll look at an eccentric case which will show how having access to these values as a function of time can be quite handy.


In [19]:
b.run_compute(irrad_method='none')


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

In [20]:
print b['rpole@primary@model']


ParameterSet: 11 parameters
  0.0@rpole@primary@latest@model: 3.56130872221 solRad
  0.1@rpole@primary@latest@model: 3.56130872221 solRad
  0.2@rpole@primary@latest@model: 3.56130872221 solRad
  0.3@rpole@primary@latest@model: 3.56130872221 solRad
  0.4@rpole@primary@latest@model: 3.56130872221 solRad
  0.5@rpole@primary@latest@model: 3.56130872221 solRad
  0.6@rpole@primary@latest@model: 3.56130872221 solRad
  0.7@rpole@primary@latest@model: 3.56130872221 solRad
  0.8@rpole@primary@latest@model: 3.56130872221 solRad
  0.9@rpole@primary@latest@model: 3.56130872221 solRad
  1.0@rpole@primary@latest@model: 3.56130872221 solRad

In [21]:
print b['rpole@secondary@model']


ParameterSet: 11 parameters
  0.0@rpole@secondary@latest@...: 1.0 solRad
  0.1@rpole@secondary@latest@...: 1.0 solRad
  0.2@rpole@secondary@latest@...: 1.0 solRad
  0.3@rpole@secondary@latest@...: 1.0 solRad
  0.4@rpole@secondary@latest@...: 1.0 solRad
  0.5@rpole@secondary@latest@...: 1.0 solRad
  0.6@rpole@secondary@latest@...: 1.0 solRad
  0.7@rpole@secondary@latest@...: 1.0 solRad
  0.8@rpole@secondary@latest@...: 1.0 solRad
  0.9@rpole@secondary@latest@...: 1.0 solRad
  1.0@rpole@secondary@latest@...: 1.0 solRad

In [22]:
print b['pot@primary@model']


ParameterSet: 11 parameters
    0.0@pot@primary@latest@model: 3.75
    0.1@pot@primary@latest@model: 3.75
    0.2@pot@primary@latest@model: 3.75
    0.3@pot@primary@latest@model: 3.75
    0.4@pot@primary@latest@model: 3.75
    0.5@pot@primary@latest@model: 3.75
    0.6@pot@primary@latest@model: 3.75
    0.7@pot@primary@latest@model: 3.75
    0.8@pot@primary@latest@model: 3.75
    0.9@pot@primary@latest@model: 3.75
    1.0@pot@primary@latest@model: 3.75

In [23]:
print b['pot@secondary@model']


ParameterSet: 11 parameters
  0.0@pot@secondary@latest@model: 10.9950371902
  0.1@pot@secondary@latest@model: 10.9950371902
  0.2@pot@secondary@latest@model: 10.9950371902
  0.3@pot@secondary@latest@model: 10.9950371902
  0.4@pot@secondary@latest@model: 10.9950371902
  0.5@pot@secondary@latest@model: 10.9950371902
  0.6@pot@secondary@latest@model: 10.9950371902
  0.7@pot@secondary@latest@model: 10.9950371902
  0.8@pot@secondary@latest@model: 10.9950371902
  0.9@pot@secondary@latest@model: 10.9950371902
  1.0@pot@secondary@latest@model: 10.9950371902

Eccentricity and Potentials

The parameters for 'pot' and 'rpole' are defined to be at periastron, but because of volume conservation, the actual polar radius (and potential) of a star in an eccentric orbit will change as a function of phase.

Having access to the instantaneous synthetic values of both of the parameters as a function of time in the mesh allows us to see how the radii and potentials are changing.

For more information, see the Eccentricity Tutorial


In [ ]: