This is the central location where all variables should be defined, and any relationships between them should be given. Having all definitions collected in one file is useful because other files can reference this one, so there is no need for duplication, and less room for mistakes. In particular, the relationships between variables are defined only here, so we don't need to reimplement those relationships. And if you do ever find a mistake, you only have to fix it in one place, then just re-run the other notebooks.
There are two main classes of variables:
The distinction is only required for code output, to ensure that everything gets calculated correctly. The PN equations you write down and manipulate can be in terms of any of these variables.
The fundamental variables that go into PN equations are things like the mass, spins $\chi_1$, and $\chi_2$, orbital angular velocity $\hat{\ell}$, and unit separation vector $\hat{n}$. We also include the tidal-coupling parameters in this list. Also, note that only $M_1$ is included. This is because the total mass is always assumed to be 1, so $M_2 = 1-M_1$.
The derived variables are further distinguished by whether they will need to be recalculated at each time step or not. For example, though we define the spins fundamentally as $\chi_1$ and $\chi_2$, we can also define derived spins $S$ and $\Sigma$, which need to be recalculated if the system is precessing. On the other hand, the masses are constant and fundamentally defined by $M_1$, so $M_2$ and $\nu$ only need to be calculated from that information once.
In [1]:
# Make sure division of integers does not round to the nearest integer
from __future__ import division
# Make everything in python's symbolic math package available
from sympy import * # Make sure sympy functions are used in preference to numpy
import sympy # Make sympy. constructions available
from sympy import Rational as frac # Rename for similarity to latex
from sympy import log as ln
# Print symbolic expressions nicely
init_printing()
# We'll use the numpy `array` object for vectors
from numpy import array, cross, dot
# We'll use a custom object to keep track of variables
from Utilities.PNObjects import PNCollection
PNVariables = PNCollection()
Only the most basic variables should be defined here.
Note that we will be using (quaternion) logarithmic rotors to describe the orientations of the spins, and the orientation and velocity of the binary itself. This allows us to reduce the number of constraints in the system, and only evolve the minimal number of equations. For example, the spins are constant, so only two degrees of freedom are needed. These can be expressed without ambiguities or singularities in the form of logarithmic rotors: $\mathfrak{r}_1 = \mathfrak{r}_{\chi_1 x} \hat{x} + \mathfrak{r}_{\chi_1 y} \hat{y}$, so that $\vec{\chi}_1 = \lvert \chi_1 \rvert\, e^{\mathfrak{r}_1}\, \hat{z}\, e^{-\mathfrak{r}_1}$. This may look complicated, but it performs very well numerically.
We will still be able to write and manipulate the PN equations directly in terms of familiar quantities like $\vec{S}_1 \cdot \hat{\ell}$, etc., but the fundamental objects will be the rotors, which means that the substitutions made for code output will automatically be in terms of the rotors.
In [2]:
# Dimensionful quantities, just in case anybody uses them...
PNVariables.AddBasicConstants('G, c')
# Masses of objects 1 and 2.
PNVariables.AddBasicConstants('M1')
PNVariables.AddBasicConstants('M2')
# Angular speed of separation vector
PNVariables.AddBasicVariables('v', positive=True)
# Tidal deformabilities, in units where the total mass is 1
PNVariables.AddBasicConstants('lambda1, lambda2')
# Spin vectors (assumed to be constant)
PNVariables.AddBasicVariables('chi1_x, chi1_y, chi1_z')
PNVariables.AddBasicVariables('chi2_x, chi2_y, chi2_z')
# Orbital angular-velocity unit vector ("Newtonian" angular momentum direction)
PNVariables.AddBasicVariables('ellHat_x, ellHat_y, ellHat_z')
# Orbital separation unit vector
PNVariables.AddBasicVariables('nHat_x, nHat_y, nHat_z')
Any variable that can be derived from the variables above should be put in this section.
These variables should probably be left in arbitrary form, unless a particular simplification is desired. The substitutions dictionary should map from the general names and their definitions in terms of basic variables. In numerical codes, their values can be calculated once per time step and then stored, so that the values do not have to be re-calculated every time they appear in an expression.
Various common combinations of the two masses:
In [3]:
PNVariables.AddDerivedConstant('M', M1+M2)
PNVariables.AddDerivedConstant('delta', (M1-M2)/M)
PNVariables.AddDerivedConstant('nu', M1*M2/M**2)
PNVariables.AddDerivedConstant('nu__2', (M1*M2/M**2)**2)
PNVariables.AddDerivedConstant('nu__3', (M1*M2/M**2)**3)
PNVariables.AddDerivedConstant('q', M1/M2)
The system's vector basis is given by $(\hat{\ell}, \hat{n}, \hat{\lambda})$, and will be computed by the code in terms of the fundamental logarithmic rotors defined above. Here, we give all the substitutions that will be needed in the code.
In [4]:
PNVariables.AddDerivedVariable('ellHat', array([ellHat_x, ellHat_y, ellHat_z]), datatype='std::vector<double>')
PNVariables.AddDerivedVariable('nHat', array([nHat_x, nHat_y, nHat_z]), datatype='std::vector<double>')
PNVariables.AddDerivedVariable('lambdaHat', cross(ellHat.substitution, nHat.substitution), datatype='std::vector<double>')
# Components of lambdaHat are defined in terms of components of ellHat and nHat
for i,d in zip([0,1,2],['x','y','z']):
PNVariables.AddDerivedVariable('lambdaHat_'+d, lambdaHat.substitution[i])
Various spin components and combinations:
In [5]:
PNVariables.AddDerivedVariable('chiVec1', array([chi1_x, chi1_y, chi1_z]), datatype='std::vector<double>')
PNVariables.AddDerivedVariable('chiVec2', array([chi2_x, chi2_y, chi2_z]), datatype='std::vector<double>')
PNVariables.AddDerivedVariable('chi1Mag', sqrt(chi1_x**2 + chi1_y**2 + chi1_z**2))
PNVariables.AddDerivedVariable('chi2Mag', sqrt(chi2_x**2 + chi2_y**2 + chi2_z**2))
PNVariables.AddDerivedConstant('chi1chi1', dot(chiVec1.substitution, chiVec1.substitution))
PNVariables.AddDerivedVariable('chi1chi2', dot(chiVec1.substitution, chiVec2.substitution))
PNVariables.AddDerivedConstant('chi2chi2', dot(chiVec2.substitution, chiVec2.substitution))
PNVariables.AddDerivedVariable('chi1_ell', dot(chiVec1.substitution, ellHat.substitution))
PNVariables.AddDerivedVariable('chi1_n', dot(chiVec1.substitution, nHat.substitution))
PNVariables.AddDerivedVariable('chi1_lambda', dot(chiVec1.substitution, lambdaHat.substitution))
PNVariables.AddDerivedVariable('chi2_ell', dot(chiVec2.substitution, ellHat.substitution))
PNVariables.AddDerivedVariable('chi2_n', dot(chiVec2.substitution, nHat.substitution))
PNVariables.AddDerivedVariable('chi2_lambda', dot(chiVec2.substitution, lambdaHat.substitution))
PNVariables.AddDerivedConstant('sqrt1Mchi1chi1', sqrt(1-chi1chi1))
PNVariables.AddDerivedConstant('sqrt1Mchi2chi2', sqrt(1-chi2chi2))
PNVariables.AddDerivedVariable('S', chiVec1.substitution*M1**2 + chiVec2.substitution*M2**2, datatype=chiVec1.datatype)
PNVariables.AddDerivedVariable('S_ell', chi1_ell*M1**2 + chi2_ell*M2**2)
PNVariables.AddDerivedVariable('S_n', chi1_n*M1**2 + chi2_n*M2**2)
PNVariables.AddDerivedVariable('S_lambda', chi1_lambda*M1**2 + chi2_lambda*M2**2)
PNVariables.AddDerivedVariable('Sigma', M*(chiVec2.substitution*M2 - chiVec1.substitution*M1), datatype=chiVec1.datatype)
PNVariables.AddDerivedVariable('Sigma_ell', M*(chi2_ell*M2 - chi1_ell*M1))
PNVariables.AddDerivedVariable('Sigma_n', M*(chi2_n*M2 - chi1_n*M1))
PNVariables.AddDerivedVariable('Sigma_lambda', M*(chi2_lambda*M2 - chi1_lambda*M1))
PNVariables.AddDerivedVariable('chi_s', (chiVec1.substitution + chiVec2.substitution)/2, datatype=chiVec1.datatype)
PNVariables.AddDerivedVariable('chi_s_ell', (chi1_ell+chi2_ell)/2)
PNVariables.AddDerivedVariable('chi_s_n', (chi1_n+chi2_n)/2)
PNVariables.AddDerivedVariable('chi_s_lambda', (chi1_lambda+chi2_lambda)/2)
PNVariables.AddDerivedVariable('chi_a', (chiVec1.substitution - chiVec2.substitution)/2, datatype=chiVec1.datatype)
PNVariables.AddDerivedVariable('chi_a_ell', (chi1_ell-chi2_ell)/2)
PNVariables.AddDerivedVariable('chi_a_n', (chi1_n-chi2_n)/2)
PNVariables.AddDerivedVariable('chi_a_lambda', (chi1_lambda-chi2_lambda)/2)
Other functions of the angular velocity that find frequent use:
In [6]:
PNVariables.AddDerivedVariable('x', v**2)
PNVariables.AddDerivedVariable('Omega_orb', (v**3)/M)
PNVariables.AddDerivedVariable('logv', log(v))