IPython Notebook | Python Script
As always, let’s do imports and initialize a logger and a new Bundle. See Building a System for more details.
%matplotlib inline
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)
In order for apsidal motion to be apparent, we need an eccentric system that is precessing.
b['ecc'] = 0.2
Let’s set a very noticeable rate of precession.
b['dperdt'] = 2.0 * u.deg/u.d
We’ll add lc and orb datasets to see how the apsidal motion affects each. We’ll need to sample over several orbits of the binary (which has a period of 3 days, by default).
b.add_dataset('lc', times=np.linspace(0,1,101), dataset='lc01')
b.add_dataset('lc', times=np.linspace(4,5,101), dataset='lc02')
<ParameterSet: 15 parameters | contexts: compute, dataset>
b.add_dataset('orb', times=np.linspace(0,5,401), dataset='orb01')
<ParameterSet: 3 parameters | contexts: compute, dataset>
b.run_compute(irrad_method='none')
<ParameterSet: 18 parameters | kinds: orb, lc>
Let’s plot the orbit from above and highlight the positions of each star at each cycle (times that are multiples of the period).
axs, artists = b['orb01@model'].plot(y='zs', time=[0,1,2,3,4,5])
Now looking at the light curve, we can see that this is resulting in the eclipses moving in phase-space.
axs, artists = b['lc01@model'].plot(time=[0,1])
axs, artists = b['lc02@model'].plot(time=[4,5])
axs, artists = b['lc01@model'].plot(x='phases')
axs, artists = b['lc02@model'].plot(x='phases')