Roche potential¶
Let there be two masses \(m_1\) and \(m_2\) with \(m_1 \ge m_2\) and \(m_1\) centered at \((x,y,z) = (0,0,0)\) and \((x,y,z) = (a,0,0)\), where \(a\) is the separation between the masses. Both objects are located in the x,y plane with x-axis counting distance along the connecting line and z measuring distance perpendicular to the orbital plane.
For synchronous rotation and a circular orbit, the Roche potential reads (e.g., Hilditch 2001 “An introduction to close binary stars”)
with
we write
where \(\Phi_n\) is the dimensionless Roche potential.
Note that moving material in the rotating frame is also subject to Coriolis forces, which are non-conservative and, therefore, cannot be incorporated into the scalar Roche potential.
The module provides functions to obtain the derivative of the dimensionless Roche potential with respect to the dimensionless variables representing the three directions. To turn the result into real, physical numbers, the following conversion is required
Example of usage¶
from __future__ import print_function, division
from PyAstronomy import pyasl
import numpy as np
import matplotlib.pylab as plt
x, y = np.linspace(-1.5, 2, 300), np.linspace(-1.6, 1.6, 300)
xx, yy = np.meshgrid(x, y)
# Coordinates in orbital plain
z = 0
# Mass ratio
q = 0.2
# Get dimensional values of Roche potential
p = pyasl.rochepot_dl(xx, yy, z, q)
# Positions (and potentials) of Lagrange points
l1, l1pot = pyasl.get_lagrange_1(q)
l2, l2pot = pyasl.get_lagrange_2(q)
l3, l3pot = pyasl.get_lagrange_3(q)
l4, l5 = pyasl.get_lagrange_4(), pyasl.get_lagrange_5()
l4pot = pyasl.rochepot_dl(l4[0], l4[1], l4[2], q)
l5pot = pyasl.rochepot_dl(l5[0], l5[1], l5[2], q)
print("Effective (dimensionless) radii of first and second mass")
print("According to the approximation of Eggleton 1983:")
r1eff = pyasl.roche_lobe_radius_eggleton(q, 1)
r2eff = pyasl.roche_lobe_radius_eggleton(q, 2)
print(" Reff1: %5.3f" % r1eff)
print(" Reff2: %5.3f" % r2eff)
print()
print("Roche volume and effective radius from Monte Carlo integration:")
mcvol1 = pyasl.roche_vol_MC(q, 1)
mcvol2 = pyasl.roche_vol_MC(q, 2)
print(" MC Roche lobe volume 1: %6.4f +/- %6.4f" % (mcvol1[0:2]))
print(" MC Roche lobe volume 2: %6.4f +/- %6.4f" % (mcvol2[0:2]))
print(" MC effective radius 1: %6.4f +/- %6.4f" % (mcvol1[2:]))
print(" MC effective radius 2: %6.4f +/- %6.4f" % (mcvol2[2:]))
plt.contour(p, [l5pot*1.02, l3pot, l2pot, l1pot], colors=['g', 'c', 'b', 'r'], extent=[-1.5, 2, -1.6, 1.6])
plt.text(l1, 0, 'L1', horizontalalignment='center')
plt.text(l2, 0, 'L2', horizontalalignment='center')
plt.text(l3, 0, 'L3', horizontalalignment='center')
plt.text(l4[0], l4[1], 'L4', horizontalalignment='center')
plt.text(l5[0], l5[1], 'L5', horizontalalignment='center')
plt.show()
Example: Substellar, polar, and side radii of Earth and hot Jupiter planets¶
from PyAstronomy import pyasl
from PyAstronomy import constants as PC
# Shape of the Earth and its Roche lobe (no rotation of Earth)
pc = PC.PyAConstants()
pc.setSystem("SI")
# Mass ratio (m2/m1)
q = pc.MEarth/pc.MSun
print(f"Earth/Sun mass ratio = {q}")
# Radii of Roche lobe along Earth--Sun connecting line, polar (out-of-plane),
# and side (in plane)
rads = pyasl.get_epradius_ss_polar_side(q, pot=None)
# Convert into km
rads = [r*pc.AU/1e3 for r in rads]
print(f"Roche lobe radii (substellar, polar, side) [1e6 km] = " + \
", ".join(["%g"%(r/1e6) for r in rads]))
reff_earth = pc.REarth/pc.RJ
# Small epsilon because rocky Earth is a small body on a large
# orbit (by Roche standards)
erads = pyasl.get_radius_ss_polar_side(q, 1, reff_earth, eps=1e-10)
# Convert into km
erads = [r*pc.RJ/1e3 for r in erads]
print(f"Earth radii (substellar, polar, side) [km] = " + \
", ".join(["%g"%(r) for r in erads[0:3]]))
print()
print("Roche lobe and planetary shape of typical hot Jupiter")
# Typical hot Jupiter
q = 1e-3
# sma in au
sma = 0.02
# Effective (circular disk) transit radius [RJ]
reff = 1.4
rads = pyasl.get_epradius_ss_polar_side(q, pot=None)
# Convert into km
rads = [r*pc.AU/1e3 for r in rads]
print(f"Roche lobe radii (substellar, polar, side) [1e6 km] = " + \
", ".join(["%g"%(r/1e6) for r in rads]))
jrads = pyasl.get_radius_ss_polar_side(q, sma, reff, eps=1e-10)
print(f"Hot Jupiter radii (substellar, polar, side) [RJ] = " + \
", ".join(["%g"%(r) for r in jrads[0:3]]))
With output
Earth/Sun mass ratio = 3.0032983882201428e-06
Roche lobe radii (substellar, polar, side) [1e6 km] = 1.49152, 0.952662, 0.995462
Earth radii (substellar, polar, side) [km] = 6378.14, 6378.14, 6378.14
Roche lobe and planetary shape of typical hot Jupiter
Roche lobe radii (substellar, polar, side) [1e6 km] = 10.1264, 6.52393, 6.80657
Hot Jupiter radii (substellar, polar, side) [RJ] = 1.51623, 1.38718, 1.41294
Functionality and API¶
Roche lobe potential¶
- PyAstronomy.pyasl.rochepot_dl(x, y, z, q)¶
Dimensionless Roche potential (\(\Phi_n\), synchronous rotation)
More massive component (\(m_1\)) is centered at (x,y,z) = (0,0,0). Less massive component (\(m_2\)) is at (1,0,0). The unit of length is the distance between the objects. Both objects are in the x,y plane (x-axis along the connecting line and z perpendicular to the orbital plane).
- Parameters
- x, y, zfloat or array
Location(s) at which to calculate the potential. Unit of length is the distance between the masses m1 and m2.
- qfloat
Mass ratio (0 <= m2/m1 <= 1)
- Returns
- Potentialfloat or array
The potential at the specified location(s)
- PyAstronomy.pyasl.rochepot(x, y, z, m1, m2, a)¶
Roche potential (\(\Phi\), synchronous rotation)
More massive component (\(m_1\)) is at (x,y,z) = (0,0,0). Less massive component (\(m_2\)) is at (a,0,0).
- Parameters
- x, y, zfloat or array
Position [m]
- afloat
Separation of masses [m]
- m1, m2float
Masses of objects [kg] (m1 >= m2)
- Returns
- Roche potentialfloat or array
At specified locations [m**2/s**2]
Partial derivatives¶
- PyAstronomy.pyasl.ddx_rochepot_dl(x, q, y=0, z=0)¶
Derivative of dimensionless Roche potential along x-axis
- Parameters
- x, y, zfloat or array
Position (default for y and z is zero)
- qfloat
Mass ratio m2/m1
- Returns
- Derivative: float or array
d/dx of dimensionless Roche potential
- PyAstronomy.pyasl.ddy_rochepot_dl(y, q, x=1, z=0)¶
Derivative of dimensionless Roche potential along y-axis
- Parameters
- x, y, zfloat or array
Position (default for x is one and zero for z)
- qfloat
Mass ratio m2/m1
- Returns
- Derivative: float or array
d/dy of dimensionless Roche potential
- PyAstronomy.pyasl.ddz_rochepot_dl(z, q, x=1, y=0)¶
Derivative of dimensionless Roche potential along z-axis
- Parameters
- x, y, zfloat or array
Position (default for x is one and zero for z)
- qfloat
Mass ratio m2/m1
- Returns
- Derivative: float or array
d/dz of dimensionless Roche potential
Roche lobe volume and radius¶
- PyAstronomy.pyasl.roche_lobe_radius_eggleton(q, m)¶
Approximate effective dimensionless Roche lobe radius according to Eggelton 1983 (ApJ 268, 368).
The effective Roche lobe radius is the radius of a sphere with the same volume as the actual equipotential surface defining the Roche lobe.
\[r_L \approx \frac{0.49\,q^{2/3}}{0.6\,q^{2/3} + \ln(1+q^{1/3})}\]- Parameters
- qfloat or array
Mass ratio (m2/m1)
- mint, {1,2}
Calculate radius for mass 1 or 2
- PyAstronomy.pyasl.roche_vol_MC(q, m=2, n=100000, pl=None, eps=0.0001, fullout=True)¶
Calculate (dimensionless) volume of equipotential surface such as the Roche lobe
Uses Monte Carlo (MC) integration
- Parameters
- qfloat
Mass ratio (m2/m1)
- mint, {1,2}
Whether to calculate volume for m1 or m2
- nint, optional
Number of samples to be used in the Monte Carlo integration. Default is 100000. Increase to improve accuracy.
- plfloat, optional
The (dimensionless) potential level bounding the volume. If None (default), the Roche lobe potential (L1) is used.
- epsfloat, optional
Margin used to avoid singularities of Roche potential. Default is 1e-4.
- fulloutboolean, optional
If True (default), provides more output than volume.
- Returns
- V, Verrfloats
Dimensionless volume of Roche lobe and estimate of uncertainty
- Reff, Refferr: floats, optional
Radius of a sphere with the same radius (dimensionless) and estimate of uncertainty. Provided if fullout is True (default).
- PyAstronomy.pyasl.roche_yz_extent(q, m=2, pl=None, eps=0.0001)¶
Extent of equipotential surface in y and z direction
- Returns
- dy, dzfloat
Distance to equipotential level along y and z direction.
Roche limit¶
- PyAstronomy.pyasl.roche_limit_fluid(q, R2)¶
Calculate the Roche limit for a fluid body on circular orbit
An expression for the Roche limit, \(a_R\), was given by the eponymous Édouard Albert Roche in 1849 (see, e.g., Rappaport et al. 2013, ApJL 773 for a modern reference)
\[a_R = 2.44 R_1 \sqrt[3]{\frac{\rho_1}{\rho_2}} = 2.44 R_2 q^{-1/3}\]where \(R_{1,2}\) are the radii of the primary and secondary bodies, \(\rho_{1,2}\) the respective densities, and \(q\) the mass ratio; the second equality relies on spherical shapes.
- Parameters
- qfloat
Mass ratio
- R2float
Radius of secondary body (assuming spherical shape). Units are arbitrary and same as result.
- Returns
- aRfloat
The Roche distance in the units of R2.
Lagrange points¶
- PyAstronomy.pyasl.get_lagrange_1(q, getdlrp=True, eps=0.0001)¶
Get location of first Lagrange point
- Parameters
- qfloat
Mass ratio
- getdlrpboolean, optional
If True (default), also the dimensionless Roche potential at the first Lagrange point is returned.
- epsfloat, optional
The potential diverges at x=0 and x=1. The search for a root of the derivative is performed in the interval [eps, 1-eps].
- Returns
- xL1float
The location of the first Lagrange point (xL1, 0, 0).
- Potentialfloat, optional
The dimensionless potential at the first Lagrange point (default is True).
- PyAstronomy.pyasl.get_lagrange_2(q, getdlrp=True, eps=0.0001)¶
Get location of second Lagrange point
- Parameters
- qfloat
Mass ratio
- getdlrpboolean, optional
If True (default), also the dimensionless Roche potential at the first Lagrange point is returned.
- epsfloat, optional
The potential diverges at x=0 and x=1. The search for a root of the derivative is performed in the interval [1+eps, 1+10*rL], where rL is the estimated (effective) Roche lobe radius according to Eggleton 1983.
- Returns
- xL2float
The location of the second Lagrange point (xL2, 0, 0).
- Potentialfloat, optional
The dimensionless potential at the second Lagrange point (default is True).
- PyAstronomy.pyasl.get_lagrange_3(q, getdlrp=True, eps=0.0001)¶
Get location of third Lagrange point
- Parameters
- qfloat
Mass ratio
- getdlrpboolean, optional
If True (default), also the dimensionless Roche potential at the third Lagrange point is returned.
- epsfloat, optional
The potential diverges at x=0. The search for a root of the derivative is performed in the interval [-10, -eps].
- Returns
- xL3float
The location of the first Lagrange point (xL3, 0, 0).
- Potentialfloat, optional
The dimensionless potential at the third Lagrange point (default is True).
- PyAstronomy.pyasl.get_lagrange_4()¶
Get location of forth Lagrange point
Orbital angular momentum is supposed to point into +z direction. L4 is in direction of motion of m2.
- Returns
- x,y,zfloat
The location of the fourth Lagrange point (xL3, 0, 0).
- PyAstronomy.pyasl.get_lagrange_5(getdlrp=True)¶
Get location of fifth Lagrange point
Orbital angular momentum is supposed to point into +z direction. L5 is behind motion of m2.
- Returns
- x,y,zfloat
The location of the fifth Lagrange point (xL3, 0, 0).
Coriolis force¶
- PyAstronomy.pyasl.coriolisAcceleration(omega, v)¶
Get Coriolis acceleration
In the rotating frame, moving material is subject to the Coriolis force, determined by the circular frequency (vector) of the rotating frame and the velocity vector within it.
\[a_C = -2 \vec{\omega} \times \vec{v}\]- Parameters
- omegafloat or 3-element 1d array
Circular frequency. If float is given, vector assumed to be omega*(0,0,1).
- varray
Either 3-element 1d array with velocity vector in corotating frame or array with shape (nv, 3) with nv velocity vectors (i.e., [::,0] gives all x components of velocity)
- Returns
- Coriolis accelerationarray
Either 3-element array or [vn,3] array, depending in input shape of v (units determined by input).
Radii in the Roche geometry¶
- PyAstronomy.pyasl.get_epradius_ss_polar_side(q, pot=None, eps=1e-06)¶
Get dimensionless radii of equipotential surface for secondary mass in three directions
Assumes that the equipotential surface is enclosed within the Roche lobe
- Parameters
- qfloat
Mass ratio m2/m1
- potfloat, optional
Dimensionless potential for equipotential surface. If None (default), the L1 potential will be adopted.
- epsfloat, optional
Lowest dimensionless radius considered in the search.
- brlboolean, optional
Search below Roche lobe (default is True).
- Returns
- rxfloat
Substellar radius, height along connecting line (dimensionless)
- rzfloat
Polar radius (along z-axis)
- ryfloat
Side radius (along y-axis)
- PyAstronomy.pyasl.get_radius_ss_polar_side(q, sma, reff, eps=1e-06)¶
Calculate substellar, polar, and side radii of secondary body
The function first finds the dimensionless Roche potential corresponding to the equipotential surface corresponding to the effective radius and then determines the radii of that surface along the x, y, and z directions.
The effective radius is the radius of a sphere with the same area as the ellipse spanned by the polar and side radii of the equipotential area.
- Parameters
- qfloat
Mass ratio m2/m1
- smafloat
Semi-major axis [AU]
- refffloat
Effective radius of the body [RJ]
- Returns
- rss, rpol, rsidefloat
Radius of equipotential surface at substellar point (along the connecting x-axis), polar radius (along the z-axis), and radius toward the side (y-axis).
- pequifloat
Dimensionless potential of equipotential surface corresponding to effective radius.