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 [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')
SVec = Vector('S', r'\vec{S}', [S_n, S_lambda, S_ell])(t)
SigmaVec = Vector('SigmaVec', r'\vec{\Sigma}', [Sigma_n, Sigma_lambda, Sigma_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_lambda*sin(Omega*t) + S_n*cos(Omega*t)])(t)
SigmaVectimeslambdaHat = Vector('SigmaVectimeslambdaHat', r'\vec{\Sigma} \times \hat{\lambda}',
[-Sigma_ell*cos(Omega*t), -Sigma_ell*sin(Omega*t), Sigma_lambda*sin(Omega*t) + Sigma_n*cos(Omega*t)])(t)
nHattimesSVec = Vector('nHattimesSVec', r'\hat{n} \times \vec{S}',
[S_ell*sin(Omega*t), -S_ell*cos(Omega*t), S_lambda*cos(Omega*t) - S_n*sin(Omega*t)])(t)
nHattimesSigmaVec = Vector('nHattimesSigmaVec', r'\hat{n} \times \vec{\Sigma}',
[Sigma_ell*sin(Omega*t), -Sigma_ell*cos(Omega*t), Sigma_lambda*cos(Omega*t) - Sigma_n*sin(Omega*t)])(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_n*sin(Omega*t)+S_lambda*cos(Omega*t)) # Symbol(r'\vec{S} \cdot \vec{v}')
SigmaVecdotv = v*(-Sigma_n*sin(Omega*t)+Sigma_lambda*cos(Omega*t)) # Symbol(r'\vec{\Sigma} \cdot \vec{v}')
nHatdotSVec = S_n*cos(Omega*t)+S_lambda*sin(Omega*t) # Symbol(r'\hat{n} \cdot \vec{S}')
nHatdotSigmaVec = Sigma_n*cos(Omega*t)+Sigma_lambda*sin(Omega*t) # 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 [7]:
%%time
MathKernel='/Applications/Local/Mathematica.app/Contents/MacOS/MathKernel'
# This little `math` script does the conversion we need.
# It has to output to a temporary file because I know of no
# better way to get the correct format out of Mathematica.
# This whole thing
MathCommand = r"""
\[Nu] = nu;
\[Mu] = mu;
\[CapitalDelta] = delta;
\[Delta] = delta;
\[Lambda] = lambda;
\[CapitalSigma] = Sigma;
\[Gamma] = gamma;
\[Epsilon] = epsilon;
\[Gamma]\[Nu] = gamma*nu;
\[ScriptL] = ell;
(*Print[InputForm[ToExpression["{0}", TeXForm]]];*)
ToExpression["{0}", TeXForm] >> /tmp/Math.cpp;
Exit[];
"""
# # I believe all the results from this group:
# SrcDir = "source/BlanchetLivRev/"
# MultipoleNames = ["I_ij", "I_ij_A", "I_ij_B", "I_ij_C", "J_ij", "I_ijk", "J_ijk", "I_ijkl", "J_ijkl"]
# # I believe all the results from this group:
# SrcDir = "source/BoheEtAl/eqs/"
# MultipoleNames = ["JijCMS", "IijklCMS", "JijklCMS","JijkCMS",]
# # I mostly believe the results from this group, but
# # I don't have the patience to check them all:
# SrcDir = "source/BoheEtAl/eqs/"
# MultipoleNames = ["IijCMS", "IijkCMS", "JijklCMS",]
SrcDir = "source/BoheEtAl/eqs/"
MultipoleNames = ["IijklCMS", "JijklCMS","IijkCMS", "JijkCMS", "IijCMS", "JijCMS",]
multiplication = r'\,'
def TensorCreator(coefficient, vector):
"Used in the regex wars (below)"
def func(match):
indices = list(match.group(1))
if not coefficient:
return vector + '_{%s}'%match.group(1)
return vector + '_{%s}'%match.group(1) + multiplication + coefficient + '^{%s}'%len(indices)
return func
prereplacements = [
[r'\\\\',r''], # Line breask are meaningless
[r'&',r''], # Alignments are meaningless
[r'\\nonumber',r''], # I wasn't gonna anyway
[r'~',r''], # Explicit spacing is useless
[r'\\q*uad',r''], # Explicit spacing is useless
[r'.*?= *',r''],
[r' *\+ *\\mathcal{O}.*',r''],
[r'\(n,\\*(Sigma|S) *,v\)',
r'v\,\mathrm{ndot\1dotlambdaHat}'], # Catch Bohé et al.'s use of (n,S,v), (n,Sigma,v), for the triple product
[r'\(\\*(n|v|S|Sigma) *\\*(n|v|S|Sigma) *\)',
r'\mathrm{\1dot\2}'], # Catch Bohé et al.'s use of (nv), etc., for the dot product
#[r'\\,', multiplication], # This gets us a good chunk of the way, if we're lucky
[r'(\\langle|\\rangle|<|>)',r''], # The YlmTensors are already STF
[r'\\epsilon_\{ab *i\}(.*?)(\{*.\}*)(_.*?)a(.*?)(\{*.\}*)(_.*?)b',
r'r\,v\,\ell_i\, \1\2\3 \4\5\6 '], # Do the cross product
[r'\{*x\}*_(\{\}| )', ''], # Remove stubs left from above
[r'\{*v\}*_(\{\}| )', ''], # Remove stubs left from above
[r'\\left(\(|\[|\\{)',r'('],
[r'\\right(\)|\]|\\})',r')'],
[r'\\(B|b)ig+(\(|\[|\\{)',r'('],
[r'\\(B|b)ig+(\)|\]|\\})',r')'],
[r'\\left\.',r''], # Do this BEFORE replacing \left
[r'\\right\.',r''], # Do this BEFORE replacing \right
[r'\\left',r''],
[r'\\right',r''],
[r'\(\\mathbf{\\*(.*?)}\\times\\mathbf{\\*(.*?)}\)',
r'\mathrm{\1times\2}'],
[r'\^\{([ijkl])\}',r'_{\1}'],
[r'n_{ *(.*?) *}',TensorCreator('', 'n')],
[r'x_{ *(.*?) *}',TensorCreator('r', 'n')],
[r'\{x\}_{ *(.*?) *}',TensorCreator('r', 'n')],
[r'v_{ *(.*?) *}',TensorCreator('v', r'\lambda')],
[r'\{v\}_{ *(.*?) *}',TensorCreator('v', r'\lambda')],
[r' *',r' '], # Tighten up spacing
[r'\(',r'\,('], # Doesn't hurt to assume multiplication
[r'([^\\]),',r'\1'], # Remove comma at the end of an equation
[r'([^\\])\.',r'\1'], # Remove period at the end of an equation
[r'\\ln *\\*,*',r'\log'],
['\\\\',r'\\\\'], # Escape backslashes for Mathematica
[r"'",r"Prm"],
#[r"'",r"'\''"], # Pull the ugly command-line single-quote hack AFTER escaping backslashes
]
postreplacements = [
[r'\^',r'**'],
[r'lambda(?!Hat)',r'lambdaHat'],
[r'\[',r'('],
[r'\]',r')'],
[r'Sigma',r'SigmaVec'],
[r'S(?!igma|ubscript)',r'SVec'],
[r'n(?!u)',r'nHat'],
[r'ell',r'ellHat'],
[r'Subscript\(((SigmaVec|SVec|nHat)times)v', r'v*Subscript(\1lambdaHat'],
#[r'Subscript',r'TensorProductCreator'],
]
for MultipoleName in MultipoleNames:
print("\n"+MultipoleName)
with open(SrcDir+MultipoleName+".tex", "r") as myfile:
data=myfile.read().replace('\n', ' ')
#print("Raw content of {0}:\n".format(MultipoleName)); print(data)
for replacement in prereplacements:
data = re.sub(replacement[0], replacement[1], data)
#print("\nbecomes\n"); print("""'{0}'""".format(data))
#! echo '{MathCommand.format(data)}'
#continue # Skip Mathematica stuff for now
! {MathKernel} -run '{MathCommand.format(data)}' >& /dev/null
MathTmp = !cat /tmp/Math.cpp
Math[MultipoleName] = ''.join(MathTmp)
for replacement in postreplacements:
Math[MultipoleName] = re.sub(replacement[0], replacement[1], Math[MultipoleName])
#print("\nbecomes\n"); print(Math[MultipoleName]); print("\nbecomes\n")
Sym[MultipoleName] = sympify(Math[MultipoleName], ns)
#display(Sym[MultipoleName])
with open('results/raw/'+MultipoleName+'.tex', 'w') as the_file:
the_file.write(latex(Sym[MultipoleName]))
#print(" reduces to")
Sym[MultipoleName] = ReduceExpr(Sym[MultipoleName]).compress()
#display(Sym[MultipoleName])
with open('results/'+MultipoleName+'.tex', 'w') as the_file:
the_file.write(latex(Sym[MultipoleName]))
with open('results/'+MultipoleName+'.py', 'w') as the_file:
the_file.write(repr(Sym[MultipoleName]))
#print("\n\n")
In [8]:
%%time
AsymmetricModeTerms = {}
for I_L,J_L in [["IijklCMS", "JijklCMS"], ["IijkCMS", "JijkCMS"], ["IijCMS", "JijCMS"]]:
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(sympy.I*sphericalharmonictensors.Vlm(V_L, m)/sympy.sqrt(2))
else :
tmp = sympy.simplify(-sphericalharmonictensors.Ulm(U_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, 3)
AsymmetricModeTerms[ell][m] = sympy.simplify(sympy.simplify(tmp).doit().subs(t,0))
display(AsymmetricModeTerms[ell][m])
with open('results/hAsymmetric_{0}_{1}.tex'.format(ell,m), 'w') as the_file:
the_file.write(latex(AsymmetricModeTerms[ell][m]))
with open('results/hAsymmetric_{0}_{1}.py'.format(ell,m), 'w') as the_file:
the_file.write(repr(AsymmetricModeTerms[ell][m]))
with open('results/hAsymmetric_{0}_{1}.c'.format(ell,m), 'w') as the_file:
the_file.write(ccode(N(AsymmetricModeTerms[ell][m])))
In [8]:
-272*S_lambda*delta + 250*i*S_n*delta + 319*Sigma_lambda*nu - 125*Sigma_lambda -239*i*Sigma_n*nu + 19*i*Sigma_n - 293*i*Sigma_n*nu
In [9]:
LeadingOrder22 = (2*nu*v.subs(t,0)**2)*sqrt(16*pi/5)
( DropHigherRelativePowers(AsymmetricModeTerms[2][2],0) ) / ( LeadingOrder22 )
Out[9]:
In [10]:
LeadingOrder44 = (2*nu*v.subs(t,0)**2)*sqrt(16*pi/5)*v.subs(t,0)**2*8*sqrt(frac(5,7))*(-1 + 3*nu)/9
( DropHigherRelativePowers(AsymmetricModeTerms[4][4],0) ) / ( LeadingOrder44 )
Out[10]:
In [10]:
In [10]: