To do:
Note: At the moment, this notebook computes radiated spherical harmonic modes from source multipole tensors. It is designed only for the relatively low-order requirements of the spin parts of the waveform. Thus, it generally ignores terms like $\hat{n} \cdot \vec{v}$. (That is, the velocity is assumed to be orthogonal to the separation between the black holes.) Also, it ignores the intermediate multipoles needed to translate between source and radiated multipoles; they enter at somewhat higher order than we need.
It would be easy enough to modify this notebook to account for the unaccounted-for terms, but this has not been done yet.
We have a couple formulas at lowest order that can guide us in our choice of terms to keep:
$dv/dt = -\mathcal{L}/(dE/dv) \approx -(2v^{10}/5)/(-v/8) = v^9/20$
$d^2v/dt^2 = 9v^8/20\, dv/dt = 9 v^{17}/400$
\begin{align*} \hat{n} \cdot \vec{v} &\approx v\, \frac{dr}{dt} \\ &= v\, \frac{d(1/v^2)}{dt} \\ &= -\frac{2}{v^2}\, \frac{dv}{dt} \\ &= -\frac{v^7}{10} \end{align*}Clearing the cache is not strictly necessary, but I have found the cache can sometimes cause trouble.
In [1]:
from sympy.core.cache import clear_cache
clear_cache()
In [2]:
from __future__ import division, print_function
import sympy
from sympy import *
from sympy.printing import ccode
from sympy import Rational as frac
import re
import simpletensors
from simpletensors import Vector, xHat, yHat, zHat
from simpletensors import TensorProduct, SymmetricTensorProduct, Tensor, ReduceExpr
import sphericalharmonictensors
init_printing()
var('mu,ell,A,B,C,epsilon')
var('jk,jkl,jklm,jklmn,jklmno,jklmnop,jklmnopq') # Index sets Blanchet uses
var('i,j,k,l,m,n,o') # Indices Bohé et al. use
#var('Subscript', cls=Function); # Mathematica will return this, so let's just pretend it exists
var('vartheta, varphi')
var('nu, m, delta, c, t, v0')
# These are related scalar functions of time
var('v, gamma, r', cls=Function)
v = v(t)
x = v**2
Omega = v**3
# gamma = v**2*(1 + (1-nu/3)*v**2 + (1-65*nu/12)*v**4) # APPROXIMATELY!!! Change this later
#r = 1/gamma
gamma = gamma(t)
r = r(t)
#r = 1/v**2 # APPROXIMATELY!!! Just for testing. Change this later
# These get redefined momentarily, but have to exist first
var('nHat, lambdaHat, ellHat', cls=Function)
# And now we define them as vector functions of time
nHat = Vector('nHat', r'\hat{n}', [cos(Omega*t),sin(Omega*t),0,])(t)
lambdaHat = Vector('lambdaHat', r'\hat{\lambda}', [-sin(Omega*t),cos(Omega*t),0,])(t)
ellHat = Vector('ellHat', r'\hat{\ell}', [0,0,1,])(t)
# # These are the spin functions -- first, the individual components as regular sympy.Function objects; then the vectors themselves
# var('S_n, S_lambda, S_ell', cls=Function)
# var('Sigma_n, Sigma_lambda, Sigma_ell', cls=Function)
# SigmaVec = Vector('SigmaVec', r'\vec{\Sigma}', [Sigma_n(t), Sigma_lambda(t), Sigma_ell(t)])(t)
# SVec = Vector('S', r'\vec{S}', [S_n(t), S_lambda(t), S_ell(t)])(t)
# These are the spin vector functions, treated as constant at this order
# (note that this is w.r.t. the x,y,z basis, rather than n,lambda,ell)
var('S_n, S_lambda, S_ell')
var('Sigma_n, Sigma_lambda, Sigma_ell')
SigmaVec = Vector('SigmaVec', r'\vec{\Sigma}', [Sigma_n*cos(Omega*t) - Sigma_lambda*sin(Omega*t),
Sigma_n*sin(Omega*t) + Sigma_lambda*cos(Omega*t), Sigma_ell])(t)
SVec = Vector('S', r'\vec{S}', [S_n*cos(Omega*t) - S_lambda*sin(Omega*t),
S_n*sin(Omega*t) + S_lambda*cos(Omega*t), S_ell])(t)
# These are some cross products that appear in the spin terms, assuming $\vec{v} = v\hat{\lambda}$
SVectimeslambdaHat = Vector('SVectimeslambdaHat', r'\vec{S} \times \hat{\lambda}',
[-S_ell*cos(Omega*t), -S_ell*sin(Omega*t), S_n])(t)
SigmaVectimeslambdaHat = Vector('SigmaVectimeslambdaHat', r'\vec{\Sigma} \times \hat{\lambda}',
[-Sigma_ell*cos(Omega*t), -Sigma_ell*sin(Omega*t), Sigma_n])(t)
nHattimesSVec = Vector('nHattimesSVec', r'\hat{n} \times \vec{S}',
[S_ell*sin(Omega*t), -S_ell*cos(Omega*t), S_lambda])(t)
nHattimesSigmaVec = Vector('nHattimesSigmaVec', r'\hat{n} \times \vec{\Sigma}',
[Sigma_ell*sin(Omega*t), -Sigma_ell*cos(Omega*t), Sigma_lambda])(t)
nHattimeslambdaHat = ellHat
nHatdotSVecdotlambdaHat = -S_ell # Symbol(r'\left[ \hat{n} \cdot \left( \vec{S} \times \hat{\lambda} \right) \right]')
nHatdotSigmaVecdotlambdaHat = -Sigma_ell # Symbol(r'\left[ \hat{n} \cdot \left( \vec{\Sigma} \times \hat{\lambda} \right) \right]')
nHatdotv = diff(r,t,1) # Symbol(r'\hat{n} \cdot \hat{v}')
# These ignore the v^7 contribution from v along nHat
SVecdotv = v*S_lambda # Symbol(r'\vec{S} \cdot \vec{v}')
SigmaVecdotv = v*Sigma_lambda # Symbol(r'\vec{\Sigma} \cdot \vec{v}')
nHatdotSVec = S_n # Symbol(r'\hat{n} \cdot \vec{S}')
nHatdotSigmaVec = Sigma_n # Symbol(r'\hat{n} \cdot \vec{\Sigma}')
In [3]:
SimplifyVectorDerivatives = {diff(nHat,t):Omega*lambdaHat, diff(lambdaHat,t):-Omega*nHat}
DropvDerivatives = {diff(v,t,vDerivs)**vDerivPower:0
for vDerivs in range(1,3)
for vDerivPower in range(1,5)
if vDerivs+vDerivPower>=2} # dv/dt \approx v^9 / 20
DropvDerivatives.update({diff(gamma,t,gammaDerivs)**gammaDerivPower:0
for gammaDerivs in range(1,3)
for gammaDerivPower in range(1,5)
if gammaDerivs+gammaDerivPower>=2})
DropvDerivatives.update({diff(r,t,rDerivs)**rDerivPower:0
for rDerivs in range(1,3)
for rDerivPower in range(1,5)
if rDerivs+rDerivPower>=2})
DropvDerivatives.update({'nHatdotv':0}) # n dot v \approx - v^7 / 10
ApproxvDerivatives = {diff(v,t): v**9/20,
diff(gamma,t): v**10/10,
diff(r,t): -v**6/10,
}
SpinComponents = {'SVecdotv':v0*S_lambda, 'SigmaVecdotv':v0*Sigma_lambda,
'nHatdotSVec':S_n, 'nHatdotSigmaVec':Sigma_n, }
# powers.get(v,0) + 2*powers.get(x,0) + 3*powers.get(Omega,0) + 2*powers.get(gamma,0) - 2*powers.get(r,0)
In [4]:
def Subscript(vec, indices):
#print("Subscript({0}, {1})".format(vec, indices))
try:
indices_str = indices.name
except AttributeError:
indices_str = indices
#print("About to construct TensorProduct([{0},]*{1}, coefficient=1, symmetric=True)".format(vec, len(indices_str)))
return TensorProduct([vec,]*len(indices_str), coefficient=1, symmetric=True)
In [5]:
Math={}
Sym={}
ns = locals()
In [6]:
def CollectivePower(expr):
"""Given some product, return 2 * total PN order
This treats v as 1, x**2 as 2, r as -2, etc., not accounting
for any contribution to PN order from spins
"""
powers = expr.subs({v.subs(t,0):v,r.subs(t,0):r,}).as_powers_dict()
return powers.get(v,0) + 2*powers.get(x,0) + 3*powers.get(Omega,0) + 2*powers.get(gamma,0) - 2*powers.get(r,0) \
+ 9*powers.get(diff(v,t,1), 0) + 17*powers.get(diff(v,t,2), 0) + 25*powers.get(diff(v,t,3), 0) \
+ 10*powers.get(diff(gamma,t,1), 0) + 18*powers.get(diff(gamma,t,2), 0) + 26*powers.get(diff(gamma,t,3), 0) \
+ 6*powers.get(diff(r,t,1), 0) + 14*powers.get(diff(r,t,2), 0) + 22*powers.get(diff(r,t,3), 0)
def DropHigherRelativePowers(expression, relative_order=3):
try:
expression = sympy.expand(expression)
if(expression.func==sympy.Add):
lowest_order = min( [ CollectivePower(arg) for arg in expression.args ] )
expression = sum( [ arg for arg in expression.args if CollectivePower(arg)<=relative_order+lowest_order ] )
try:
return expression.simplify()
except:
if(expression==0.0) :
return sympy.sympify(0)
return expression
else:
return expression
except:
return expression
def DropHigherAbsolutePowers(expression, absolute_order=3):
try:
expression = sympy.expand(expression)
if(expression.func==sympy.Add):
expression = sum( [ arg for arg in expression.args if CollectivePower(arg)<=absolute_order ] )
try:
return expression.simplify()
except:
if(expression==0.0) :
return sympy.sympify(0)
return expression
else:
if(CollectivePower(expression)<=absolute_order) :
return expression
else:
return sympy.sympify(0)
except:
return expression
def DropHigherAbsolutePowersOfTensorProduct(tp, absolute_order=3):
return simpletensors.TensorProduct([c for c in tp],
coefficient=DropHigherAbsolutePowers(tp.coefficient, absolute_order), symmetric=tp.symmetric)
def DropHigherRelativePowersOfTensorProduct(tp, relative_order=3):
return simpletensors.TensorProduct([c for c in tp],
coefficient=DropHigherRelativePowers(tp.coefficient, relative_order), symmetric=tp.symmetric)
def DropHigherAbsolutePowersOfTensor(t, absolute_order=3):
return sum( [DropHigherAbsolutePowersOfTensorProduct(tp, absolute_order) for tp in t ] )
def DropHigherRelativePowersOfTensor(t, relative_order=3):
lowest_order = sys.maxint
for term in t:
expression = sympy.expand(term.coefficient)
try:
if(expression.func==sympy.Add):
for arg in expression.args:
lowest_order = min(lowest_order, CollectivePower(arg))
else:
lowest_order = min(lowest_order, CollectivePower(expression))
except:
pass
return sum( [DropHigherAbsolutePowersOfTensorProduct(tp, relative_order+lowest_order) for tp in t ] )
In [24]:
Sym['Iij'] = nu*r**2*nHat*nHat
Sym['Jij'] = -nu*delta*r**2*v*nHat*ellHat
Sym['Iijk'] = -nu*delta*r**3*nHat*nHat*nHat
Sym['Jijk'] = nu*r**3*v*nHat*nHat*ellHat
In [25]:
%%time
AsymmetricModeTerms = {}
for I_L,J_L in [["Iijk", "Jijk"], ["Iij", "Jij"]]:
U_L,V_L = [Sym[I_L],Sym[J_L]]
ell = U_L.rank
print('ell={0}'.format(ell)); sys.stdout.flush()
AsymmetricModeTerms[ell] = {}
for m in range(0,ell+1):
print('(ell,m)=({0},{1})'.format(ell,m)); sys.stdout.flush()
if((ell+m)%2==0) :
tmp = sympy.simplify(-sphericalharmonictensors.Ulm(U_L, m)/sympy.sqrt(2))
else :
tmp = sympy.simplify(sympy.I*sphericalharmonictensors.Vlm(V_L, m)/sympy.sqrt(2))
for i_order in range(ell):
tmp = diff(tmp, t)
tmp = tmp.xreplace(SimplifyVectorDerivatives).subs(dict(DropvDerivatives.items()+ApproxvDerivatives.items())).subs(SpinComponents)
tmp = DropHigherRelativePowers(tmp, 1)
AsymmetricModeTerms[ell][m] = sympy.simplify(sympy.simplify(tmp).doit().subs(t,0))
display(AsymmetricModeTerms[ell][m])
with open('results/hSymmetric_{0}_{1}.tex'.format(ell,m), 'w') as the_file:
the_file.write(latex(AsymmetricModeTerms[ell][m]))
with open('results/hSymmetric_{0}_{1}.py'.format(ell,m), 'w') as the_file:
the_file.write(repr(AsymmetricModeTerms[ell][m]))
with open('results/hSymmetric_{0}_{1}.c'.format(ell,m), 'w') as the_file:
the_file.write(ccode(N(AsymmetricModeTerms[ell][m])))
In [10]: