Skip to content

Instantly share code, notes, and snippets.

@petigura

petigura/Precession-Terms.ipynb

Created October 31, 2018 19:39
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save petigura/4e6aaa8eaf3b8a299c46a86dcae18c9e to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save petigura/4e6aaa8eaf3b8a299c46a86dcae18c9e to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
Precession-Terms.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Hey, K.\n",
"\n",
"I'm making slow, but steady progress on K2-19. The system of two sub-saturns near the 3:2 commensurability with 8 and 12 day periods. There is also a 1 Re planet at 2.5 days. My REBOUND analysis of the system (without including GR or the inner planet) strongly suggested that the planets are in a *secular*, but not *mean motion resonance*, with aligned apses. \n",
"\n",
"When we chatted with Sean, you recommended checking the effects of GR and J2 precession from the inner planet. I could add these into my REBOUND simulations, but I'd like to gain some intuition here. I computed the apsidal precession frequencies GR and J2 terms and got\n",
"\n",
"- GR 10^-4 rad/year\n",
"- J2 3 x 10^-4 rad/year\n",
"\n",
"However, the secular oscilations of planets c and d occur at ~ 1 rad/year. My intution says that GR and J2 precession only matter if they can differentially precess K2-19b/c out of their secular resonance. Since the secular cycles are so quick, my guess is that these other terms don't matter. \n",
"\n",
"Does this sound resonable? Would you mind looking over my calculations in the ipython notebook \n",
"\n",
"Thanks!\n",
"\n",
"Erik\n",
"\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Populating the interactive namespace from numpy and matplotlib\n"
]
}
],
"source": [
"%pylab inline\n",
"from astropy import constants as c\n",
"from astropy import units as u"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### GR precession\n",
"\n",
"(pulled from Yee et al. 2018)\n",
"\n",
"$\\dot{\\omega} = 3 n \\frac{G M_\\star}{ac^2}$"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" \n",
"for planet b\n",
"omegadot_gr is 9.47156781521e-05 1 / yr\n",
"or one cycle per 66337.3311554 yr\n"
]
}
],
"source": [
"# GR precession\n",
"Mstar = 0.8 * c.M_sun\n",
"per_b = 8*u.day\n",
"\n",
"def omegadot_gr(per):\n",
" n = 2 * np.pi / per\n",
" a = (c.G* Mstar * per_b**2 / 4 / np.pi**2)**(1/3.)\n",
" _omegadot_gr = (3 * n * c.G * Mstar / a / c.c**2).to(u.yr**-1)\n",
" return _omegadot_gr\n",
" \n",
"print \"\"\n",
"print \"for planet b\"\n",
"print \"omegadot_gr is \",omegadot_gr(7.9*u.day)\n",
"print \"or one cycle per\", 2 * pi / omegadot_gr(7.9*u.day)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### J2 precession\n",
"\n",
"Approximate planet d as a wire, equivalent to assuming that the central body has a size $a_d$ and a J2 moment of \n",
"\n",
"$\\frac{1}{2}\\frac{m_d}{M_\\star}$\n",
"\n",
"Object orbits will precess at\n",
"\n",
"$\\dot{\\omega} = 3 n J_2 \\left(\\frac{a_d}{a}\\right)^2$"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" \n",
"for planet b\n",
"omegadot_j2 is 0.000346923811727 1 / yr\n",
"or one cycle per 18111.1387999 yr\n"
]
}
],
"source": [
"m_d = 1 * c.M_earth\n",
"j2 = 0.5 * m_d / Mstar\n",
"\n",
"per_d = 2.5 * u.day\n",
"a_d = (c.G* Mstar * per_d**2 / 4 / pi**2)**(1/3.)\n",
"\n",
"def omegadot_j2(per):\n",
" n = 2 * np.pi / per\n",
" a = (c.G* Mstar * per_b**2 / 4 / np.pi**2)**(1/3.)\n",
" _omegadot_j2 = ((3 * n * j2) * (a_d/a)**2).to(u.yr**-1)\n",
" return _omegadot_j2\n",
"\n",
"\n",
"print \"\"\n",
"print \"for planet b\"\n",
"print \"omegadot_j2 is \", omegadot_j2(7.9*u.day)\n",
"print \"or one cycle per\", 2 * pi / omegadot_j2(7.9*u.day)"
]
},
{
"cell_type": "code",
"execution_count": 54,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Difference in GR precession rates 3.23611900353e-05 1 / yr\n",
"Difference in J2 precession rates 0.00011853230234 1 / yr\n"
]
}
],
"source": [
"print \"Difference in GR precession rates\",omegadot_gr(7.9*u.day) -omegadot_gr(12*u.day)\n",
"print \"Difference in J2 precession rates\", omegadot_j2(7.9*u.day) -omegadot_j2(12*u.day)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# both timescales are much slower than secular oscilations."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.13"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment