Package ivs :: Package timeseries :: Module keplerorbit

Module keplerorbit

Fit and evaluate binary (Keplerian) orbits.

Be careful with the inclination angle: if you only want an image projected onto the sky, it is not case-sensitive. Only for the radial velocities this becomes important (the primary and secondary can be mixed).

Example usage: compute the absolute orbits of the two components of mu Cassiopeia, and compute the relative orbit of the second component.

See ps file on http://www.chara.gsu.edu/~gudehus/binary.html and Drummond et al., 1995 (the latter's orbital parameters are used).

Neccesary imports:

>>> from ivs.units.constants import *

Orbital elements, Euler angles and masses of the system

>>> P = 21.753*365 # period in days
>>> e = 0.561 # ecc
>>> omega,Omega,i = 332.7, 47.3, 106.8 # degrees
>>> T0 = 1975.738*365 # year in days
>>> plx = 133.2 # mas
>>> M1,M2 = 0.719,0.168 # Msun
>>> times = np.linspace(T0,T0+0.95*P,100)

Derive semi-major axis of absolute orbits (semi-major axis of relative orbit is a+a2)

>>> fA = M2**3/ (M1+M2)**2 * Msol
>>> a1 = ((P*24*3600)**2 * GG * fA / (4*np.pi**2))**(1/3.)
>>> a1 = a1/au
>>> a2 = a1*M1/M2

Convert the angles to radians

>>> omega,Omega,i = omega/180*np.pi,Omega/180*np.pi,i/180*np.pi

Compute the orbits of the primary and secondary in the plane of the sky

>>> x1,y1,z1 = orbit_on_sky(times,(P,e,a1,T0,omega,Omega,i),component='prim',distance=(1000./133.,'pc'))
>>> x2,y2,z2 = orbit_on_sky(times,(P,e,a2,T0,omega,Omega,i),component='sec',distance=(1000./133.,'pc'))

Plot it:

>>> import pylab as pl
>>> colors = times/365.
>>> p = pl.title('mu Cassiopeia')
>>> p = pl.scatter(x1,y1,c=colors,edgecolors='none',cmap=pl.cm.spectral,label='_nolegend_')
>>> p = pl.scatter(x2,y2,c=colors,edgecolors='none',cmap=pl.cm.spectral,label='_nolegend_')
>>> p = pl.scatter(x2-x1,y2-y1,c=colors,edgecolors='none',cmap=pl.cm.spectral,label='_nolegend_')
>>> p = pl.plot(x1,y1,'b-',label='Primary')
>>> p = pl.plot(x2,y2,'g-',label='Secondary')
>>> p = pl.plot(x2-x1,y2-y1,'r-',label='Secondary relative')
>>> p = pl.colorbar()
>>> p = pl.legend(loc='lower right')
>>> p = pl.grid()
>>> p = pl.xlabel(r'Angular size (arcsec) --> N')
>>> p = pl.ylabel(r'Angular size (arcsec) --> E')

]include figure]]ivs_binary_keplerorbit_muCas.png]

Example usage 2: compute the absolute orbits of the two components of Capella, and compute the relative orbit of the primary component.

See Hummel et al. 1994 (?).

Neccesary imports:

>>> from ivs.units.constants import *

Orbital elements, Euler angles and masses of the system

>>> P = 104.022 # period in days
>>> e = 0.000 #ecc
>>> i = 137.18 # degree
>>> T0 = 2447528.45 # year
>>> omega = 0 # degrees
>>> Omega = 40.8 # degrees
>>> plx = 76.20 # mas
>>> M1 = 2.69 #Msun
>>> M2 = 2.56 #Msun
>>> RV0 = 0.
>>> times = np.linspace(T0,T0+0.95*P,100)

Derive semi-major axis of absolute orbits (semi-major axis of relative orbit is a+a2)

>>> fA = M2**3/ (M1+M2)**2 * Msol
>>> a1 = ((P*24*3600)**2 * GG * fA / (4*np.pi**2))**(1/3.)
>>> a1 = a1/au
>>> a2 = a1*M1/M2

Convert the angles to radians

>>> omega,Omega,i = omega/180*np.pi,Omega/180*np.pi,i/180*np.pi

Compute the orbits of the primary and secondary in the plane of the sky

>>> x1,y1,z1 = orbit_on_sky(times,(P,e,a1,T0,omega,Omega,i),component='prim',distance=(1000./plx,'pc'))
>>> x2,y2,z2 = orbit_on_sky(times,(P,e,a2,T0,omega,Omega,i),component='sec',distance=(1000./plx,'pc'))

Plot it:

>>> import pylab as pl
>>> colors = times/365.
>>> p = pl.figure()
>>> p = pl.title('Capella')
>>> p = pl.scatter(x1,y1,c=colors,edgecolors='none',cmap=pl.cm.spectral,label='_nolegend_')
>>> p = pl.scatter(x2,y2,c=colors,edgecolors='none',cmap=pl.cm.spectral,label='_nolegend_')
>>> p = pl.scatter(x1-x2,y1-y2,c=colors,edgecolors='none',cmap=pl.cm.spectral,label='_nolegend_')
>>> p = pl.plot(x1,y1,'b-',label='Primary')
>>> p = pl.plot(x2,y2,'g-',label='Secondary')
>>> p = pl.plot(x1-x2,y1-y2,'r-',label='Secondary relative')
>>> p = pl.colorbar()
>>> p = pl.legend(loc='lower right')
>>> p = pl.grid()
>>> p = pl.xlabel(r'Angular size DEC (arcsec) --> N')
>>> p = pl.ylabel(r'Angular size RA (arcsec) --> E')

]include figure]]ivs_binary_keplerorbit_capella.png]

Functions

[hide private]

Radial velocities

ndarray

radial_velocity(parameters, times=None, theta=None, itermax=8)
Evaluate Radial velocities due to kepler orbit.

source code

2xarray

orbit_in_plane(times, parameters, component='primary', coordinate_frame='polar')
Construct an orbit in the orbital plane.

source code

2xarray

velocity_in_plane(times, parameters, component='primary', coordinate_frame='polar')
Calculate the velocity in the orbital plane.

source code

project_orbit(r, theta, parameters)
Project an orbit onto the plane of the sky.

source code

orbit_on_sky(times, parameters, distance=None, component='primary')
Construct an orbit projected on the sky.

source code

float,float

true_anomaly(M, e, itermax=8)
Calculation of true and eccentric anomaly in Kepler orbits.

source code

float,float

calculate_phase(T, e, omega, pshift=0)
Compute orbital phase from true anomaly T

source code

float,float

calculate_critical_phases(omega, e, pshift=0)
Compute phase of superior conjunction and periastron passage.

source code

float

eclipse_separation(e, omega)
Calculate the eclipse separation between primary and secondary in a light curve.

source code

float

omega_from_eclipse_separation(separation, e)
Caculate the argument of periastron from the eclipse separation and eccentricity.

source code

Kepler's laws

third_law(M=None, a=None, P=None)
Kepler's third law.

source code

Variables

[hide private]

logger = logging.getLogger('IVS.KEPLER')

Function Details

[hide private]

radial_velocity(parameters, times=None, theta=None, itermax=8)

source code

Evaluate Radial velocities due to kepler orbit.

These parameters define the Keplerian orbit if you give times points (times, days):

Period of the system (days)
time of periastron passage T0 (not x0!) (HJD)
eccentricity
omega, the longitude of periastron gives the orientation of the orbit within its own plane (radians)
the semiamplitude of the velocity curve (km/s)
systemic velocity RV0 (RV of centre of mass of system) (km/s)

These parameters define the Keplerian orbit if you give angles (theta, radians):

Period of the system (days)
eccentricity
a, the semi-major axis of the absolute orbit of the component (au)
omega, the longitude of periastron gives the orientation of the orbit within in its own plane (radians)
inclination of the orbit (radians)
systemic velocity RV0 (RV of centre of mass of system) (km/s)

The periastron passage T0 can be derived via x0 by calculating

T0 = x0/(2pi*Freq) + times[0]

See e.g. p 41,42 of Hilditch, 'An Introduction To Close Binary Stars'

Parameters:

parameters (iterable) - parameters of Keplerian orbit (dependent on input)
times (numpy array) - observation times (days)
theta (numpy array) - orbital angles (radians)
itermax (integer) - number of iterations to find true anomaly

Returns: ndarray

fitted radial velocities (km/s)

orbit_in_plane(times, parameters, component='primary', coordinate_frame='polar')

source code

Construct an orbit in the orbital plane.

Give times in days

Parameters contains:

period (days)
eccentricity
semi-major axis (au)
time of periastron passage T0 (not x0!) (HJD)

Return r (m) and theta (radians)

Parameters:

times (array) - times of observations (days)
parameters (list) - list of parameters (P,e,a,T0)
component (string, one of 'primary' or 'secondary') - component to calculate the orbit of. If it's the secondary, the angles will be shifted by 180 degrees
coordinate_frame (str, ('polar' or 'cartesian')) - type of coordinates

Returns: 2xarray

coord1,coord2

velocity_in_plane(times, parameters, component='primary', coordinate_frame='polar')

source code

Calculate the velocity in the orbital plane.

Parameters:

times (array) - times of observations (days)
parameters (list) - list of parameters (P,e,a,T0)
component (string, one of 'primary' or 'secondary') - component to calculate the orbit of. If it's the secondary, the angles will be shifted by 180 degrees
coordinate_frame (str, ('polar' or 'cartesian')) - type of coordinates

Returns: 2xarray

coord1,coord2

project_orbit(r, theta, parameters)

source code

Project an orbit onto the plane of the sky.

Parameters contains the Euler angles:

omega: the longitude of periastron gives the orientation of the orbit within in its own plane (radians)
Omega: PA of ascending node (radians)
i: inclination (radians), i=pi/2 is edge on

Returns x,y (orbit in plane of the sky) and z

See Hilditch p41 for a sketch of the coordinates. The difference with this is that all angles are inverted. In this approach, North is in the positive X direction, East is in the negative Y direction.

orbit_on_sky(times, parameters, distance=None, component='primary')

source code

Construct an orbit projected on the sky.

Parameters contains:

period (days)
eccentricity
semi-major axis (au)
time of periastron passage T0 (not x0!) (HJD)
omega: the longitude of periastron gives the orientation of the orbit within in its own plane (radians)
Omega: PA of ascending node (radians)
i: inclination (radians), i=pi/2 is edge on

You can give an extra parameter 'distance' as a tuple (value,'unit'). This will be used to convert the distances to angular scale (arcsec).

Else, this function returns the distances in AU.

true_anomaly(M, e, itermax=8)

source code

Calculation of true and eccentric anomaly in Kepler orbits.

M is the phase of the star, e is the eccentricity

See p.39 of Hilditch, 'An Introduction To Close Binary Stars'

Parameters:

M (float) - phase
e (float) - eccentricity
itermax (integer) - maximum number of iterations

Returns: float,float

eccentric anomaly (E), true anomaly (theta)

calculate_phase(T, e, omega, pshift=0)

source code

Compute orbital phase from true anomaly T

Parameters:

T (float) - true anomaly
omega (float) - argument of periastron (radians)
e (float) - eccentricity
pshift (float) - phase shift

Returns: float,float

phase of superior conjunction, phase of periastron passage

calculate_critical_phases(omega, e, pshift=0)

source code

Compute phase of superior conjunction and periastron passage.

Example usage: >>> omega = np.pi/4.0 >>> e = 0.3 >>> print calculate_critical_phases(omega,e) (-0.125, -0.057644612788576133, -0.42054512757020118, -0.19235538721142384, 0.17054512757020118)

Parameters:

omega (float) - argument of periastron (radians)
e (float) - eccentricity
pshift (float) - phase shift

Returns: float,float

phase of superior conjunction, phase of periastron passage

eclipse_separation(e, omega)

source code

Calculate the eclipse separation between primary and secondary in a light curve.

Minimum separation at omega=pi Maximum spearation at omega=0

Parameters:

e (float) - eccentricity
omega (float) - argument of periastron (radians)

Returns: float

separation in phase units (0.5 is half)

omega_from_eclipse_separation(separation, e)

source code

Caculate the argument of periastron from the eclipse separation and eccentricity.

separation in phase units.

Parameters:

separation (float) - separation in phase units (0.5 is half)
e (float) - eccentricity

Returns: float

omega (longitude of periastron) in radians

third_law(M=None, a=None, P=None)

source code

Kepler's third law.

Give two quantities, derived the third.

M = total mass system (solar units) a = semi-major axis (au) P = period (d)

>>> print third_law(M=1.,a=1.)
365.256891359 
>>> print third_law(a=1.,P=365.25)
1.00003773538
>>> print third_law(M=1.,P=365.25)
0.999987421856