#!/usr/bin/env python
# coding: utf-8

# <script async src="https://www.googletagmanager.com/gtag/js?id=UA-59152712-8"></script>
# <script>
#   window.dataLayer = window.dataLayer || [];
#   function gtag(){dataLayer.push(arguments);}
#   gtag('js', new Date());
# 
#   gtag('config', 'UA-59152712-8');
# </script>
# 
# # $M\Omega$, the orbital angular velocity, up to and including 3.5 PN order
# 
# ## Author: Zach Etienne
# 
# ## This notebook uses NRPy to construct the orbital angular velocity up to and including 3.5 post-Newtonian order. 
# 
# **Notebook Status:** <font color='green'><b> Validated </b></font>
# 
# **Validation Notes:** All expressions in this notebook were transcribed twice by hand on separate occasions, and expressions were corrected as needed to ensure consistency with published work. Published work was cross-validated and typo(s) in published work were corrected. In addition, all expressions in this notebook were validated against those in the Mathematica notebook used by [Ramos-Buades, Husa, and Pratten (2018)](https://arxiv.org/abs/1810.00036) (thanks to Toni Ramos-Buades for sharing this!) Finally, this tutorial notebook has been confirmed to be self-consistent with its corresponding NRPy+ module, as documented [below](#code_validation). **Additional validation tests may have been performed, but are as yet, undocumented.**
# 
# 
# ### This notebook exists as the following Python module:
# 1. [PN_MOmega.py](../../edit/NRPyPN/PN_MOmega.py)
# 
# ### This notebook & corresponding Python module depend on the following NRPy+/NRPyPN Python modules:
# 1. [indexedexp.py](../../edit/indexedexp.py): [**documentation+tutorial**](../Tutorial-Indexed_Expressions.ipynb)
# 1. [NRPyPN_shortcuts.py](../../edit/NRPyPN/NRPyPN_shortcuts.py): [**documentation**](NRPyPN_shortcuts.ipynb)

# <a id='toc'></a>
# 
# # Table of Contents
# $$\label{toc}$$
# 
# 1. Part 1: [$M\Omega$](#momega), up to and including 3.5PN order, as derived in [Ramos-Buades, Husa, and Pratten (2018)](https://arxiv.org/abs/1810.00036)
# 1. Part 2: [Validation against second transcription and corresponding Python module](#code_validation)
# 1. Part 3: [Validation against trusted numerical values](#code_validationv2) (i.e., in Table V of [Ramos-Buades, Husa, and Pratten (2018)](https://arxiv.org/abs/1810.00036))
# 1. Part 4: [LaTeX PDF output](#latex_pdf_output): $\LaTeX$ PDF Output

# <a id='momega'></a>
# 
# # Part 1: $M\Omega$, up to and including 3.5PN order, as derived in [Ramos-Buades, Husa, and Pratten (2018)](https://arxiv.org/abs/1810.00036) \[Back to [top](#toc)\]
# $$\label{momega}$$ 
# 
# As described in the [nonspinning Hamiltonian notebook](PN-Hamiltonian-Nonspinning.ipynb), the basic physical system assumes two point particles of mass $m_1$ and $m_2$ with corresponding momentum vectors $\mathbf{P}_1$ and $\mathbf{P}_2$, and displacement vectors $\mathbf{X}_1$ and $\mathbf{X}_2$ with respect to the center of mass. Here we also consider the spin vectors of each point mass $\mathbf{S}_1$ and $\mathbf{S}_2$, respectively.
# 
# To reduce the possibility of copying error, the equation for $M\Omega$ is taken directly from the arXiv LaTeX source code of Eq A1 in [Ramos-Buades, Husa, and Pratten (2018)](https://arxiv.org/abs/1810.00036), and only mildly formatted to (1) improve presentation in Jupyter notebooks, (2) to ensure some degree of consistency in notation across different terms in other NRPyPN notebooks, and (3) to correct any errors. In particular, the boxed negative sign at 2.5PN order ($a_5$ below) was missing in the original equation. We will later show that this negative sign is necessary for consistency with the expression up to 3PN order in [Healy, Lousto, Nakano, and Zlochower (2017)](https://arxiv.org/abs/1702.00872):
# $$
# M\Omega = \frac{1}{r^{3/2}}\left(1 + \sum_{k=2}^7 \frac{a_k}{r^{k/2}}\right),
# $$
# where all terms in boxes should be replaced by 1:
# \begin{align}
# a_2 &= -\left[ \frac{  \left(3 q^2+5 q+3\right)}{2 (q+1)^2}\right] \\
# a_3 &= -\frac{(3 q+4) \chi _{1z}}{4 (q+1)^2 }-  \frac{q (4 q+3) \chi _{2z}}{4 (q+1)^2 } \\
# a_4 &= -\frac{3 q^2 \chi _{2x}^2}{2 (q+1)^2}+ \frac{3 q^2 \chi _{2y}^2}{4 (q+1)^2 }
# +\frac{3 q^2 \chi _{2z}^2}{4 (q+1)^2 }+\frac{24 q^4+103 q^3+164 q^2+103 q+24}{16 (q+1)^4 } \\
# &\quad -\frac{3 \chi _{1x}^2}{2 (q+1)^2 }-\frac{3 q \chi _{1x}\chi _{2x}}{(q+1)^2 }+\frac{3 \chi _{1y}^2}{4 (q+1)^2 }+\frac{3 q \chi _{1y} \chi _{2y}}{2 (q+1)^2 }+\frac{3 \chi _{1z}^2}{4 (q+1)^2 }+\frac{3 q \chi _{1z} \chi _{2z}}{2 (q+1)^2} \\
# a_5 &= \frac{3 \left(13 q^3+34 q^2+30 q+16\right) \chi _{1z}}{16 (q+1)^4}+  \frac{3 q \left(16 q^3+30 q^2+34 q+13\right) \chi _{2z}}{16 (q+1)^4 }\\
# a_6 &=  \frac{\left(155 q^2+180 q+76\right) \chi _{1x}^2}{16 (q+1)^4 \boxed{r^3}}+\frac{q \left(120 q^2+187 q+120\right) \chi _{1x} \chi _{2x}}{8 (q+1)^4 \boxed{r^3}}-\frac{\left(55 q^2+85 q+43\right) \chi _{1y}^2}{8 (q+1)^4 \boxed{r^3}} \\
# &  -\frac{q \left(54 q^2+95 q+54\right) \chi _{1y} \chi _{2y}}{4 (q+1)^4 \boxed{r^3}}-\frac{q \left(96 q^2+127 q+96\right) \chi _{1z} \chi _{2z}}{16 (q+1)^4 \boxed{r^3}}+\frac{q^2 \left(76 q^2+180 q+155\right) \chi _{2x}^2}{16 (q+1)^4 \boxed{r^3}} \\
# &  -\frac{q^2 \left(43 q^2+85 q+55\right) \chi _{2y}^2}{8 (q+1)^4 \boxed{r^3}}-\frac{q^2 (2 q+5) (14 q+27) \chi _{2z}^2}{32 (q+1)^4 \boxed{r^3}} -\frac{(5 q+2) (27 q+14) \chi _{1z}^2}{32 (q+1)^4 \boxed{r^3}} \\
# & +\frac{501 \pi ^2 q (q+1)^4-4 \left(120 q^6+2744 q^5+10049 q^4+14820 q^3+10049 q^2+2744 q+120\right)}{384 (q+1)^6 \boxed{r^3}} \\
# a_7 &= \frac{3 (4 q+1) q^3 \chi _{2 x}^2 \chi _{2 z}}{2 (q+1)^4}-\frac{3 (4 q+1) q^3 \chi _{2 y}^2 \chi _{2 z}}{8 (q+1)^4}-\frac{3 (4 q+1) q^3 \chi _{2 z}^3}{8 (q+1)^4}+\chi _{1x} \left(\frac{9 (2 q+1) q^2 \chi _{2 x} \chi _{2 z}}{4 (q+1)^4}+\frac{9 (q+2) q \chi _{2 x} \chi _{\boxed{?}z}}{4 (q+1)^4}\right) \\
# & +\chi _{1y} \left(\frac{9 q^2 \chi _{2 y} \chi _{1z}}{4 (q+1)^4}+\frac{9 q^2 \chi _{2 y} \chi _{2 z}}{4 (q+1)^4}\right) \\
# & +\chi _{1z} \left(\frac{9 q^2 (2 q+3) \chi _{2 x}^2}{4 (q+1)^4}-\frac{9 q^2 (q+2) \chi _{2 y}^2}{4 (q+1)^4}-\frac{9 q^2 \chi _{2 z}^2}{4 (q+1)^3}-\frac{135 q^5+385 q^4+363 q^3+377 q^2+387 q+168}{32 (q+1)^6}\right) \\
# & -\frac{\left(168 q^5+387 q^4+377 q^3+363 q^2+385 q+135\right) q \chi _{2 z}}{32 (q+1)^6}+\chi _{1x}^2 \left(\frac{3 (q+4) \chi _{1z}}{2 (q+1)^4}+\frac{9 q (3 q+2) \chi _{2 z}}{4 (q+1)^4}\right)\\
# &+\chi _{1y}^2 \left(-\frac{3 (q+4) \chi _{1z}}{8 (q+1)^4}-\frac{9 q (2 q+1) \chi _{2 z}}{4 (q+1)^4}\right)-\frac{9 q \chi _{1z}^2 \chi _{2 z}}{4 (q+1)^3}-\frac{3 (q+4) \chi _{1z}^3}{8 (q+1)^4},
# \end{align}

# Let's divide and conquer, tackling the coefficients one at a time:
# \begin{align}
# a_2 &= -\left[ \frac{  \left(3 q^2+5 q+3\right)}{2 (q+1)^2}\right] \\
# a_3 &= -\frac{(3 q+4) \chi _{1z}}{4 (q+1)^2 }-  \frac{q (4 q+3) \chi _{2z}}{4 (q+1)^2 } \\
# a_4 &= -\frac{3 q^2 \chi _{2x}^2}{2 (q+1)^2}+ \frac{3 q^2 \chi _{2y}^2}{4 (q+1)^2 }
# +\frac{3 q^2 \chi _{2z}^2}{4 (q+1)^2 }+\frac{24 q^4+103 q^3+164 q^2+103 q+24}{16 (q+1)^4 } \\
# &\quad -\frac{3 \chi _{1x}^2}{2 (q+1)^2 }-\frac{3 q \chi _{1x}\chi _{2x}}{(q+1)^2 }+\frac{3 \chi _{1y}^2}{4 (q+1)^2 }+\frac{3 q \chi _{1y} \chi _{2y}}{2 (q+1)^2 }+\frac{3 \chi _{1z}^2}{4 (q+1)^2 }+\frac{3 q \chi _{1z} \chi _{2z}}{2 (q+1)^2}
# \end{align}

# In[1]:


# Step 0: Add NRPy's directory to the path
# https://stackoverflow.com/questions/16780014/import-file-from-parent-directory
import sympy as sp               # SymPy: The Python computer algebra package upon which NRPy+ depends
import indexedexpNRPyPN as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
from NRPyPN_shortcuts import div # NRPyPN: shortcuts for e.g., vector operations

# Step 1: Construct terms a_2, a_3, and a_4, from
#  Eq A2 of Ramos-Buades, Husa, and Pratten (2018)
#    https://arxiv.org/abs/1810.00036
#  These terms have been independently validated
#    against the same terms in Eq 6 of
#  Healy, Lousto, Nakano, and Zlochower (2017)
#    https://arxiv.org/abs/1702.00872
def MOmega__a_2_thru_a_4(m1,m2, chi1x,chi1y,chi1z, chi2x,chi2y,chi2z):
    q = m2/m1 # It is assumed that q >= 1, so m2 >= m1.
    global a_2,a_3,a_4
    a_2 = -((3*q**2+5*q+3)/(2*(q+1)**2))
    a_3 = (-(3*q+4)*chi1z/(4*(q+1)**2) - q*(4*q+3)*chi2z/(4*(q+1)**2))
    a_4 = (-3*q**2*chi2x**2/(2*(q+1)**2)
           +3*q**2*chi2y**2/(4*(q+1)**2)
           +3*q**2*chi2z**2/(4*(q+1)**2)
           +(+24*q**4 + 103*q**3 + 164*q**2 + 103*q + 24)/(16*(q+1)**4)
           -3*chi1x**2/(2*(q+1)**2)
           -3*q*chi1x*chi2x/(q+1)**2
           +3*chi1y**2/(4*(q+1)**2)
           +3*q*chi1y*chi2y/(2*(q+1)**2)
           +3*chi1z**2/(4*(q+1)**2)
           +3*q*chi1z*chi2z/(2*(q+1)**2))


# In[2]:


# Second version, for validation purposes only.
def MOmega__a_2_thru_a_4v2(m1,m2, chi1x,chi1y,chi1z, chi2x,chi2y,chi2z):
    q = m2/m1 # It is assumed that q >= 1, so m2 >= m1.
    global a_2v2,a_3v2,a_4v2
    a_2v2 = - (3*q**2+5*q+3)/(2*(q+1)**2)
    a_3v2 = +(-(3*q+4)*chi1z/(4*(q+1)**2)
                  -q*(4*q+3)*chi2z/(4*(q+1)**2))
    a_4v2 = +(-(3*q**2*chi2x**2)/(2*(q+1)**2) + 3*q**2*chi2y**2/(4*(q+1)**2) + 3*q**2*chi2z**2/(4*(q+1)**2)
                  +(24*q**4+103*q**3+164*q**2+103*q+24)/(16*(q+1)**4) - 3*chi1x**2/(2*(q+1)**2)
                  -3*q*chi1x*chi2x/(q+1)**2 + 3*chi1y**2/(4*(q+1)**2) + 3*q*chi1y*chi2y/(2*(q+1)**2)
                  +3*chi1z**2/(4*(q+1)**2) + 3*q*chi1z*chi2z/(2*(q+1)**2))


# In[3]:


# Third version, directly from Toni Ramos-Buades' Mathematica notebook (thanks Toni!)
def MOmega__a_2_thru_a_4v3(m1,m2, chi1x,chi1y,chi1z, chi2x,chi2y,chi2z):
    q = m2/m1 # It is assumed that q >= 1, so m2 >= m1.
    global a_2v3,a_3v3,a_4v3
    a_2v3 = ( -(3 + 5*q + 3*q**2)/(2*(1 + q)**2) )
    a_3v3 = ( (-4*chi1z - 3*chi1z*q - 3*chi2z*q - 4*chi2z*q**2)/(4*(1 + q)**2) )
    a_4v3 = ( (-3*chi1x**2)/(2*(1 + q)**2) + (3*chi1y**2)/(4*(1 + q)**2) + (3*chi1z**2)/(4*(1 + q)**2) - (3*chi1x*chi2x*q)/(1 + q)**2 + (3*chi1y*chi2y*q)/(2*(1 + q)**2) + (3*chi1z*chi2z*q)/(2*(1 + q)**2) -
       (3*chi2x**2*q**2)/(2*(1 + q)**2) + (3*chi2y**2*q**2)/(4*(1 + q)**2) + (3*chi2z**2*q**2)/(4*(1 + q)**2) + (24 + 103*q + 164*q**2 + 103*q**3 + 24*q**4)/(16*(1 + q)**4) )


# Next, $a_5$ and $a_6$:
# 
# \begin{align}
# a_5 &= \frac{3 \left(13 q^3+34 q^2+30 q+16\right) \chi _{1z}}{16 (q+1)^4}+  \frac{3 q \left(16 q^3+30 q^2+34 q+13\right) \chi _{2z}}{16 (q+1)^4 }\\
# a_6 &=  \frac{\left(155 q^2+180 q+76\right) \chi _{1x}^2}{16 (q+1)^4 \boxed{r^3}}+\frac{q \left(120 q^2+187 q+120\right) \chi _{1x} \chi _{2x}}{8 (q+1)^4 \boxed{r^3}}-\frac{\left(55 q^2+85 q+43\right) \chi _{1y}^2}{8 (q+1)^4 \boxed{r^3}} \\
# &  -\frac{q \left(54 q^2+95 q+54\right) \chi _{1y} \chi _{2y}}{4 (q+1)^4 \boxed{r^3}}-\frac{q \left(96 q^2+127 q+96\right) \chi _{1z} \chi _{2z}}{16 (q+1)^4 \boxed{r^3}}+\frac{q^2 \left(76 q^2+180 q+155\right) \chi _{2x}^2}{16 (q+1)^4 \boxed{r^3}} \\
# &  -\frac{q^2 \left(43 q^2+85 q+55\right) \chi _{2y}^2}{8 (q+1)^4 \boxed{r^3}}-\frac{q^2 (2 q+5) (14 q+27) \chi _{2z}^2}{32 (q+1)^4 \boxed{r^3}} -\frac{(5 q+2) (27 q+14) \chi _{1z}^2}{32 (q+1)^4 \boxed{r^3}} \\
# & +\frac{501 \pi ^2 q (q+1)^4-4 \left(120 q^6+2744 q^5+10049 q^4+14820 q^3+10049 q^2+2744 q+120\right)}{384 (q+1)^6 \boxed{r^3}} \\
# \end{align}

# In[4]:


# Construct terms a_5 and a_6, from
#  Eq A1 of Ramos-Buades, Husa, and Pratten (2018)
#    https://arxiv.org/abs/1810.00036
#  These terms have been independently validated
#    against the same terms in Eq 6 of
#  Healy, Lousto, Nakano, and Zlochower (2017)
#    https://arxiv.org/abs/1702.00872
def MOmega__a_5_thru_a_6(m1,m2, chi1x,chi1y,chi1z, chi2x,chi2y,chi2z):
    q = m2/m1 # It is assumed that q >= 1, so m2 >= m1.
    global a_5,a_6
    a_5 = (+3*  (13*q**3 + 34*q**2 + 30*q + 16)*chi1z/(16*(q+1)**4)
           +3*q*(16*q**3 + 30*q**2 + 34*q + 13)*chi2z/(16*(q+1)**4))
    a_6 = (+(+155*q**2 + 180*q + 76)*chi1x**2/(16*(q+1)**4)
           +q*(+120*q**2 + 187*q + 120)*chi1x*chi2x/(8*(q+1)**4)
           -(+55*q**2 + 85*q + 43)*chi1y**2/(8*(q+1)**4)
           -q*(+54*q**2 + 95*q + 54)*chi1y*chi2y/( 4*(q+1)**4)
           -q*(+96*q**2 +127*q + 96)*chi1z*chi2z/(16*(q+1)**4)
           +q**2*(+76*q**2 + 180*q + 155)*chi2x**2/(16*(q+1)**4)
           -q**2*(+43*q**2 +  85*q +  55)*chi2y**2/( 8*(q+1)**4)
           -q**2*(+2*q+5)*(+14*q+27)*chi2z**2/(32*(q+1)**4)
           -     (+5*q+2)*(+27*q+14)*chi1z**2/(32*(q+1)**4)
           +(+501*sp.pi**2*q*(q+1)**4
             -4*(120*q**6 + 2744*q**5 + 10049*q**4 + 14820*q**3 + 10049*q**2 + 2744*q + 120))/(384*(q+1)**6))


# In[5]:


# Second version, for validation purposes only.
def MOmega__a_5_thru_a_6v2(m1,m2, chi1x,chi1y,chi1z, chi2x,chi2y,chi2z):
    q = m2/m1 # It is assumed that q >= 1, so m2 >= m1.
    pi = sp.pi
    global a_5v2,a_6v2
    a_5v2 = +(+3*  (13*q**3+34*q**2+30*q+16)*chi1z/(16*(q+1)**4)
              +3*q*(16*q**3+30*q**2+34*q+13)*chi2z/(16*(q+1)**4))
    a_6v2 =+(+(155*q**2+180*q+76)*chi1x**2   /(16*(q+1)**4) + q*(120*q**2+187*q+120)*chi1x*chi2x/(8*(q+1)**4)
                  -( 55*q**2+ 85*q+43)*chi1y**2   /( 8*(q+1)**4) - q*( 54*q**2+ 95*q+ 54)*chi1y*chi2y/(4*(q+1)**4)
                  -q   *(96*q**2+127*q+ 96)*chi1z*chi2z/(16*(q+1)**4)
                  +q**2*(76*q**2+180*q+155)*chi2x**2   /(16*(q+1)**4)
                  -q**2*(43*q**2+ 85*q+ 55)*chi2y**2   /( 8*(q+1)**4)
                  -q**2*(2*q+5)*(14*q+27)  *chi2z**2   /(32*(q+1)**4)
                  -     (5*q+2)*(27*q+14)  *chi1z**2   /(32*(q+1)**4)
                  +(501*sp.pi**2*q*(q+1)**4 - 4*(120*q**6+2744*q**5+10049*q**4+14820*q**3+10049*q**2+2744*q+120))
                   /(384*(q+1)**6))


# In[6]:


# Third version, directly from Toni Ramos-Buades' Mathematica notebook (thanks Toni!)
def MOmega__a_5_thru_a_6v3(m1,m2, chi1x,chi1y,chi1z, chi2x,chi2y,chi2z):
    q = m2/m1 # It is assumed that q >= 1, so m2 >= m1.
    Pi = sp.pi
    global a_5v3,a_6v3
    a_5v3 = ( (3*(16*chi1z + 30*chi1z*q + 13*chi2z*q + 34*chi1z*q**2 + 34*chi2z*q**2 + 13*chi1z*q**3 + 30*chi2z*q**3 + 16*chi2z*q**4))/(16*(1 + q)**4) )
    a_6v3 =  ( (167*Pi**2*q)/(128.*(1 + q)**2) - (chi2z**2*q**2*(135 + 124*q + 28*q**2))/(32.*(1 + q)**4) - (chi2y**2*q**2*(55 + 85*q + 43*q**2))/(8.*(1 + q)**4) - (chi1y*chi2y*q*(54 + 95*q + 54*q**2))/(4.*(1 + q)**4) -
       (chi1y**2*(43 + 85*q + 55*q**2))/(8.*(1 + q)**4) + (chi2x**2*q**2*(155 + 180*q + 76*q**2))/(16.*(1 + q)**4) - (chi1z*chi2z*q*(96 + 127*q + 96*q**2))/(16.*(1 + q)**4) +
       (chi1x*chi2x*q*(120 + 187*q + 120*q**2))/(8.*(1 + q)**4) - (chi1z**2*(28 + 124*q + 135*q**2))/(32.*(1 + q)**4) + (chi1x**2*(76 + 180*q + 155*q**2))/(16.*(1 + q)**4) -
       (120 + 2744*q + 10049*q**2 + 14820*q**3 + 10049*q**4 + 2744*q**5 + 120*q**6)/(96.*(1 + q)**6) )


# Finally, $a_7$:
# 
# \begin{align}
# a_7 &= \frac{3 (4 q+1) q^3 \chi _{2 x}^2 \chi _{2 z}}{2 (q+1)^4}-\frac{3 (4 q+1) q^3 \chi _{2 y}^2 \chi _{2 z}}{8 (q+1)^4}-\frac{3 (4 q+1) q^3 \chi _{2 z}^3}{8 (q+1)^4}+\chi _{1x} \left(\frac{9 (2 q+1) q^2 \chi _{2 x} \chi _{2 z}}{4 (q+1)^4}+\frac{9 (q+2) q \chi _{2 x} \chi _{\boxed{?}z}}{4 (q+1)^4}\right) \\
# & +\chi _{1y} \left(\frac{9 q^2 \chi _{2 y} \chi _{1z}}{4 (q+1)^4}+\frac{9 q^2 \chi _{2 y} \chi _{2 z}}{4 (q+1)^4}\right) \\
# & +\chi _{1z} \left(\frac{9 q^2 (2 q+3) \chi _{2 x}^2}{4 (q+1)^4}-\frac{9 q^2 (q+2) \chi _{2 y}^2}{4 (q+1)^4}-\frac{9 q^2 \chi _{2 z}^2}{4 (q+1)^3}-\frac{135 q^5+385 q^4+363 q^3+377 q^2+387 q+168}{32 (q+1)^6}\right) \\
# & -\frac{\left(168 q^5+387 q^4+377 q^3+363 q^2+385 q+135\right) q \chi _{2 z}}{32 (q+1)^6}+\chi _{1x}^2 \left(\frac{3 (q+4) \chi _{1z}}{2 (q+1)^4}+\frac{9 q (3 q+2) \chi _{2 z}}{4 (q+1)^4}\right)\\
# &+\chi _{1y}^2 \left(-\frac{3 (q+4) \chi _{1z}}{8 (q+1)^4}-\frac{9 q (2 q+1) \chi _{2 z}}{4 (q+1)^4}\right)-\frac{9 q \chi _{1z}^2 \chi _{2 z}}{4 (q+1)^3}-\frac{3 (q+4) \chi _{1z}^3}{8 (q+1)^4}
# \end{align}

# In[7]:


# Construct term a_7, from Eq A1 of
#  Ramos-Buades, Husa, and Pratten (2018)
#    https://arxiv.org/abs/1810.00036
def MOmega__a_7(m1,m2, chi1x,chi1y,chi1z, chi2x,chi2y,chi2z):
    q = m2/m1 # It is assumed that q >= 1, so m2 >= m1.
    global a_7
    a_7 = (+3*(4*q+1)*q**3*chi2x**2*chi2z/(2*(q+1)**4)
           -3*(4*q+1)*q**3*chi2y**2*chi2z/(8*(q+1)**4)
           -3*(4*q+1)*q**3*chi2z**3      /(8*(q+1)**4)
           +chi1x*(+9*(2*q+1)*q**2*chi2x*chi2z/(4*(q+1)**4)
                   +9*(1*q+2)*q   *chi2x*chi1z/(4*(q+1)**4))
           +chi1y*(+9*q**2*chi2y*chi1z/(4*(q+1)**4)
                   +9*q**2*chi2y*chi2z/(4*(q+1)**4))
           +chi1z*(+9*q**2*(2*q+3)*chi2x**2/(4*(q+1)**4)
                   -9*q**2*(  q+2)*chi2y**2/(4*(q+1)**4)
                   -9*q**2        *chi2z**2/(4*(q+1)**3)
                   -(135*q**5 + 385*q**4 + 363*q**3 + 377*q**2 + 387*q + 168)/(32*(q+1)**6))
           -(+168*q**5 + 387*q**4 + 377*q**3 + 363*q**2 + 385*q + 135)*q*chi2z/(32*(q+1)**6)
           +chi1x**2*(+3*(q+4)*chi1z/(2*(q+1)**4)
                      +9*q*(3*q+2)*chi2z/(4*(q+1)**4))
           +chi1y**2*(-3*(q+4)*chi1z/(8*(q+1)**4)
                      -9*q*(2*q+1)*chi2z/(4*(q+1)**4))
           -9*q*chi1z**2*chi2z/(4*(q+1)**3)
           -3*(q+4)*chi1z**3/(8*(q+1)**4))


# In[8]:


# Second version, for validation purposes only.
def MOmega__a_7v2(m1,m2, chi1x,chi1y,chi1z, chi2x,chi2y,chi2z):
    q = m2/m1 # It is assumed that q >= 1, so m2 >= m1.
    global a_7v2
    a_7v2  = +(+(3*(4*q+1)*q**3*chi2x**2*chi2z)/(2*(q+1)**4)
                  -(3*(4*q+1)*q**3*chi2y**2*chi2z)/(8*(q+1)**4)
                  -(3*(4*q+1)*q**3*chi2z**3)      /(8*(q+1)**4)
                  +chi1x*(+(9*(2*q+1)*q**2*chi2x*chi2z)/(4*(q+1)**4)
                          +(9*(1*q+2)*q   *chi2x*chi1z)/(4*(q+1)**4))
                  +chi1y*(+(9*q**2*chi2y*chi1z)/(4*(q+1)**4)
                          +(9*q**2*chi2y*chi2z)/(4*(q+1)**4))
                  +chi1z*(+(9*q**2*(2*q+3)*chi2x**2)/(4*(q+1)**4)
                          -(9*q**2*(1*q+2)*chi2y**2)/(4*(q+1)**4)
                          -(9*q**2        *chi2z**2)/(4*(q+1)**3)
                          -(135*q**5+385*q**4+363*q**3+377*q**2+387*q+168)/(32*(q+1)**6))
                  -(168*q**5+387*q**4+377*q**3+363*q**2+385*q+135)*q*chi2z/(32*(q+1)**6)
                  +chi1x**2*(+3*(q+4)*chi1z/(2*(q+1)**4) + 9*q*(3*q+2)*chi2z/(4*(q+1)**4))
                  +chi1y**2*(-3*(q+4)*chi1z/(8*(q+1)**4) - 9*q*(2*q+1)*chi2z/(4*(q+1)**4))
                  -9*q*chi1z**2*chi2z/(4*(q+1)**3) - 3*(q+4)*chi1z**3/(8*(q+1)**4))


# In[9]:


# Third version, directly from Toni Ramos-Buades' Mathematica notebook (thanks Toni!)
def MOmega__a_7v3(m1,m2, chi1x,chi1y,chi1z, chi2x,chi2y,chi2z):
    q = m2/m1 # It is assumed that q >= 1, so m2 >= m1.
    global a_7v3
    a_7v3 = ( (-3*(4 + q)*chi1z**3)/(8*(1 + q)**4) - (q*(135 + 385*q + 363*q**2 + 377*q**3 + 387*q**4 + 168*q**5)*chi2z)/(32*(1 + q)**6) - (9*q*chi1z**2*chi2z)/(4*(1 + q)**3) +
       (3*q**3*(1 + 4*q)*chi2x**2*chi2z)/(2*(1 + q)**4) - (3*q**3*(1 + 4*q)*chi2y**2*chi2z)/(8*(1 + q)**4) - (3*q**3*(1 + 4*q)*chi2z**3)/(8*(1 + q)**4) +
       chi1y**2*((-3*(4 + q)*chi1z)/(8*(1 + q)**4) - (9*q*(1 + 2*q)*chi2z)/(4*(1 + q)**4)) + chi1x**2*((3*(4 + q)*chi1z)/(2*(1 + q)**4) + (9*q*(2 + 3*q)*chi2z)/(4*(1 + q)**4)) +
       chi1x*((9*q*(2 + q)*chi1z*chi2x)/(4*(1 + q)**4) + (9*q**2*(1 + 2*q)*chi2x*chi2z)/(4*(1 + q)**4)) + chi1y*((9*q**2*chi1z*chi2y)/(4*(1 + q)**4) + (9*q**2*chi2y*chi2z)/(4*(1 + q)**4)) +
       chi1z*(-(168 + 387*q + 377*q**2 + 363*q**3 + 385*q**4 + 135*q**5)/(32*(1 + q)**6) + (9*q**2*(3 + 2*q)*chi2x**2)/(4*(1 + q)**4) - (9*q**2*(2 + q)*chi2y**2)/(4*(1 + q)**4) - (9*q**2*chi2z**2)/(4*(1 + q)**3)) )


# Putting it all together, recall that
# 
# $$
# M\Omega = \frac{1}{r^{3/2}}\left(1 + \sum_{k=2}^7 \frac{a_k}{r^{k/2}}\right),
# $$
# where $k/2$ is the post-Newtonian order.

# In[10]:


# Finally, sum the expressions for a_k to construct p_t as prescribed:
# MOmega = 1/r^(3/2) * (1 + \sum_{k=2}^7 (a_k/r^{k/2}))
def f_MOmega(m1,m2, chi1U,chi2U, r):
    a = ixp.zerorank1(DIM=10)
    MOmega__a_2_thru_a_4(m1,m2, chi1U[0],chi1U[1],chi1U[2], chi2U[0],chi2U[1],chi2U[2])
    a[2] = a_2
    a[3] = a_3
    a[4] = a_4
    MOmega__a_5_thru_a_6(m1,m2, chi1U[0],chi1U[1],chi1U[2], chi2U[0],chi2U[1],chi2U[2])
    a[5] = a_5
    a[6] = a_6
    MOmega__a_7(         m1,m2, chi1U[0],chi1U[1],chi1U[2], chi2U[0],chi2U[1],chi2U[2])
    a[7] = a_7
    global MOmega
    MOmega = 1 # Term prior to the sum in parentheses
    for k in range(8):
        MOmega += a[k]/r**div(k,2)
    MOmega *= 1/r**div(3,2)


# In[11]:


# Second version, for validation purposes only.
def f_MOmegav2(m1,m2, chi1U,chi2U, r):
    a = ixp.zerorank1(DIM=10)
    MOmega__a_2_thru_a_4v2(m1,m2, chi1U[0],chi1U[1],chi1U[2], chi2U[0],chi2U[1],chi2U[2])
    a[2] = a_2v2
    a[3] = a_3v2
    a[4] = a_4v2
    MOmega__a_5_thru_a_6v2(m1,m2, chi1U[0],chi1U[1],chi1U[2], chi2U[0],chi2U[1],chi2U[2])
    a[5] = a_5v2
    a[6] = a_6v2
    MOmega__a_7v2(         m1,m2, chi1U[0],chi1U[1],chi1U[2], chi2U[0],chi2U[1],chi2U[2])
    a[7] = a_7v2
    global MOmegav2
    MOmegav2 = 1 # Term prior to the sum in parentheses
    for k in range(8):
        MOmegav2 += a[k]/r**div(k,2)
    MOmegav2 *= 1/r**div(3,2)


# <a id='code_validation'></a>
# 
# # Part 2: Validation against second transcription and corresponding Python module \[Back to [top](#toc)\]
# $$\label{code_validation}$$ 
# 
# As a code validation check, we verify agreement between 
# * the SymPy expressions transcribed from the cited published work on two separate occasions, and
# * the SymPy expressions generated in this notebook, and the corresponding Python module.

# In[12]:


from NRPyPN_shortcuts import m1,m2, chi1U,chi2U, q # NRPyPN: Import needed input variables

def error(varname):
    print("ERROR: When comparing Python module & notebook, "+varname+" was found not to match.")
    sys.exit(1)

MOmega__a_2_thru_a_4v2(m1,m2, chi1U[0],chi1U[1],chi1U[2], chi2U[0],chi2U[1],chi2U[2])
MOmega__a_2_thru_a_4v3(m1,m2, chi1U[0],chi1U[1],chi1U[2], chi2U[0],chi2U[1],chi2U[2])
if sp.simplify(a_2v2 - a_2v3) != 0: error("a_2v2")
if sp.simplify(a_3v2 - a_3v3) != 0: error("a_3v2")
if sp.simplify(a_4v2 - a_4v3) != 0: error("a_4v2")

MOmega__a_5_thru_a_6v2(m1,m2, chi1U[0],chi1U[1],chi1U[2], chi2U[0],chi2U[1],chi2U[2])
MOmega__a_5_thru_a_6v3(m1,m2, chi1U[0],chi1U[1],chi1U[2], chi2U[0],chi2U[1],chi2U[2])
if sp.simplify(a_5v2 - a_5v3) != 0: error("a_5v2")
if sp.simplify(a_6v2 - a_6v3) != 0: error("a_6v2")

MOmega__a_7v2(m1,m2, chi1U[0],chi1U[1],chi1U[2], chi2U[0],chi2U[1],chi2U[2])
MOmega__a_7v3(m1,m2, chi1U[0],chi1U[1],chi1U[2], chi2U[0],chi2U[1],chi2U[2])
if sp.simplify(a_7v2 - a_7v3) != 0: error("a_7v2")

f_MOmega(m1,m2, chi1U,chi2U, q)

# Validation against second transcription of the expressions:
f_MOmegav2(m1,m2, chi1U,chi2U, q)
if sp.simplify(MOmega - MOmegav2) != 0: error("MOmegav2")

# Validation against corresponding Python module:
import PN_MOmega as MOm
MOm.f_MOmega(m1,m2, chi1U,chi2U, q)

if sp.simplify(MOmega - MOm.MOmega) != 0: error("MOm.MOmega")
print("ALL TESTS PASS")


# <a id='latex_pdf_output'></a>
# 
# # Part 3: Output this notebook to $\LaTeX$-formatted PDF file \[Back to [top](#toc)\]
# $$\label{latex_pdf_output}$$
# 
# The following code cell converts this Jupyter notebook into a proper, clickable $\LaTeX$-formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename
# [PN-MOmega.pdf](PN-MOmega.pdf) (Note that clicking on this link may not work; you may need to open the PDF file through another means.)

# In[13]:


import os,sys                    # Standard Python modules for multiplatform OS-level functions
import cmdline_helperNRPyPN as cmd    # NRPy+: Multi-platform Python command-line interface
cmd.output_Jupyter_notebook_to_LaTeXed_PDF("PN-MOmega")