The Hamiltonian expression in this module has been validated against a trusted code: the Mathematica expressions for the LALSuite SEOBNRv3_opt code.
This module documents and splits up the reduced spinning effective one-body Hamiltonian derived in Barausse and Buonanno (2010) and Pan, Buonanno, Buchman, et. al (2010) for use in efficient gravitational waveform generation software. This Cartesian formulation is used to generate gravitational waveforms in LALSuite's SEOBNRv3 approximant.
# import sympy as sp
# x,y,z,s1x,s1y,s1z,s2x,s2y,s2z = sp.symbols("x y z s1x s1y s1z s2x s2y s2z")
# blah = x+y
# g = sp.sin(blah)
# print(sp.mathematica_code(g))
import sympy as sp
m1,m2,x,y,z,px,py,pz,s1x,s1y,s1z,s2x,s2y,s2z,eta = sp.symbols("m1 m2 x y z px py pz s1x s1y s1z s2x s2y s2z eta",real=True)
c0k2,c1k2,c0k3,c1k3,c0k4,c1k4,c2k4,c0k5,c1k5,c2k5 = sp.symbols("c0k2 c1k2 c0k3 c1k3 c0k4 c1k4 c2k4 c0k5 c1k5 c2k5",real=True)
KK,k5l,b3,bb3,d1,d1v2,dheffSS,dheffSSv2 = sp.symbols("KK k5l b3 bb3 d1 d1v2 dheffSS dheffSSv2",real=True)
tortoise,copysignresult = sp.symbols("tortoise copysignresult",real=True)
sigmaKerr0 = s1x + s2x
sigmaKerr1 = s1y + s2y
sigmaKerr2 = s1z + s2z
invm1m2 = 1/(m1*m2)
m2overm1 = m2*m2*invm1m2
m1overm2 = m1*m1*invm1m2
sigmaStar0 = (m2overm1)*s1x + (m1overm2)*s2x
sigmaStar1 = (m2overm1)*s1y + (m1overm2)*s2y
sigmaStar2 = (m2overm1)*s1z + (m1overm2)*s2z
s1dots1 = s1x*s1x + s1y*s1y + s1z*s1z
s2dots2 = s2x*s2x + s2y*s2y + s2z*s2z
r2 = x*x + y*y + z*z
r = sp.sqrt(r2)
u = 1/r
u2 = u*u
u3 = u2*u
u4 = u2*u2
u5 = u4*u
etau3 = eta*u3
etau4 = eta*u4
nx = x*u
ny = y*u
nz = z*u
sKerrUSCOREx = sigmaKerr0
sKerrUSCOREy = sigmaKerr1
sKerrUSCOREz = sigmaKerr2
sStarUSCOREx = sigmaStar0
sStarUSCOREy = sigmaStar1
sStarUSCOREz = sigmaStar2
a2 = sKerrUSCOREx*sKerrUSCOREx + sKerrUSCOREy*sKerrUSCOREy + sKerrUSCOREz*sKerrUSCOREz
a4 = a2*a2
a = sp.sqrt(a2)
inva = 1/a
m1PlusetaKK = -1 + eta*KK
invm1PlusetaKK = 1/m1PlusetaKK
k0 = KK*(m1PlusetaKK - 1)
k1 = -2*(k0 + KK)*m1PlusetaKK
k2 = c0k2 + c1k2*a2
k3 = c0k3 + c1k3*a2
k4 = c0k4 + c1k4*a2 + c2k4*a4
k5 = c0k5 + c1k5*a2 + c2k5*a4
e3USCOREx = sKerrUSCOREx*inva
e3USCOREy = sKerrUSCOREy*inva
e3USCOREz = sKerrUSCOREz*inva
costheta = e3USCOREx*nx + e3USCOREy*ny + e3USCOREz*nz
xi2 = 1 - costheta*costheta
xiUSCOREx = -e3USCOREz*ny + e3USCOREy*nz
xiUSCOREy = e3USCOREz*nx - e3USCOREx*nz
xiUSCOREz = -e3USCOREy*nx + e3USCOREx*ny
#vx = -nz*xiUSCOREy + ny*xiUSCOREz
vx = sp.symbols('vx')
#vy = nz*xiUSCOREx - nx*xiUSCOREz
vy = sp.symbols('vy')
#vz = -ny*xiUSCOREx + nx*xiUSCOREy
vz = sp.symbols('vz')
w2 = r2 + a2
rho2 = r2 + a2*costheta*costheta
bulk = invm1PlusetaKK*(invm1PlusetaKK + 2*u) + a2*u2
logu = sp.log( u )
logarg = k1*u + k2*u2 + k3*u3 + k4*u4 + k5*u5 + k5l*u5*logu
onepluslogarg = (1 + logarg)
invonepluslogarg = 1/onepluslogarg
logTerms = 1 + eta*k0 + eta*sp.log(onepluslogarg)
deltaU = bulk*logTerms
#deltaT = r2*deltaU
deltaT = sp.symbols('deltaT')
deltaUUSCOREupt7 = k5 + k5l*logu
deltaUUSCOREupt6 = 4*k4 + 5*deltaUUSCOREupt7*u
deltaUUSCOREupt5 = 3*k3 + u*deltaUUSCOREupt6
deltaUUSCOREupt4 = 2*k2 + u*deltaUUSCOREupt5
deltaUUSCOREupt3 = k1 + u*deltaUUSCOREupt4
deltaUUSCOREupt2 = invm1PlusetaKK + a2*u
deltaUUSCOREupt1 = bulk*eta*deltaUUSCOREupt3
deltaUUSCOREu = 2*deltaUUSCOREupt2*logTerms + deltaUUSCOREupt1*invonepluslogarg
#deltaTUSCOREr = 2*r*deltaU - deltaUUSCOREu
deltaTUSCOREr = sp.symbols('deltaTUSCOREr')
#Lambda = w2*w2 - a2*deltaT*xi2
Lambda = sp.symbols('Lambda')
rho2xi2Lambda = rho2*xi2*Lambda
invrho2xi2Lambda = 1/rho2xi2Lambda
invrho2 = xi2*Lambda*invrho2xi2Lambda
invxi2 = rho2*Lambda*invrho2xi2Lambda
invLambda = xi2*rho2*invrho2xi2Lambda
invLambdasq = invLambda*invLambda
rho2invLambda = rho2*invLambda
#expnu = sp.sqrt(deltaT*rho2invLambda)
expnu = sp.symbols('expnu')
expMU = sp.sqrt(rho2)
#expMUexpnu = expMU*expnu
expMUexpnu = sp.symbols('expMUexpnu')
expMUsq = expMU*expMU
expnusq = expnu*expnu
expMUsqexpnusq = expMUsq*expnusq
invexpnuexpMU = 1/expMUexpnu
#invexpMU = expnu*invexpnuexpMU
invexpMU = sp.symbols('invexpMU')
invexpMUsq = invexpMU*invexpMU
expnuinvexpMU2 = expnu*invexpMUsq
invexpMUcubinvexpnu = invexpMUsq*invexpnuexpMU
DD = 1 + sp.log(1 + 6*eta*u2 + 2*(26 - 3*eta)*etau3)
deltaR = deltaT*DD
qq = 2*eta*(4 - 3*eta)
ww = 2*a*r + b3*eta*a2*a*u + bb3*eta*a*u
B = sp.sqrt(deltaT)
sqrtdeltaT = B
sqrtdeltaR = sp.sqrt(deltaR)
deltaTsqrtdeltaR = deltaT*sqrtdeltaR
sqrtdeltaTdeltaTsqrtdeltaR = sqrtdeltaT*deltaTsqrtdeltaR
invdeltaTsqrtdeltaTsqrtdeltaR = 1./sqrtdeltaTdeltaTsqrtdeltaR
#invdeltaT = sqrtdeltaT*(sqrtdeltaR*invdeltaTsqrtdeltaTsqrtdeltaR)
invdeltaT = sp.symbols('invdeltaT')
invsqrtdeltaT = deltaTsqrtdeltaR*invdeltaTsqrtdeltaTsqrtdeltaR
invsqrtdeltaR = deltaT*sqrtdeltaT*invdeltaTsqrtdeltaTsqrtdeltaR
w = ww*invLambda
LambdaUSCOREr = 4*r*w2 - a2*deltaTUSCOREr*xi2
wwUSCOREr = 2*a - (a2*a*b3*eta)*u2 - bb3*eta*a*u2
#BR = (-deltaT*invsqrtdeltaR + deltaTUSCOREr*sp.Rational(1,2))*invsqrtdeltaT
BR = sp.symbols('BR')
#wr = (-LambdaUSCOREr*ww + Lambda*wwUSCOREr)*(invLambdasq)
wr = sp.symbols('wr')
nurpt2 = w2*(-4.*r*deltaT + w2*deltaTUSCOREr)
nurpt1 = nurpt2*invdeltaT
nur = r*invrho2 + sp.Rational(1,2)*invLambda*nurpt1
mur = r*invrho2 - invsqrtdeltaR
a2costheta = a2*costheta
wcospt2 = deltaT*ww
wcospt1 = invLambdasq*wcospt2
#wcos = -2*(a2costheta)*wcospt1
wcos = sp.symbols('wcos')
nucospt3 = invrho2*invLambda
nucospt2 = w2*nucospt3
nucospt1 = a2costheta*nucospt2
nucos = (w2 - deltaT)*nucospt1
mucos = a2costheta*invrho2
etaover12r = eta*sp.Rational(1,12)*u
csi = sp.sqrt(deltaT*deltaR)/w2
csi1 = 1 + (1-sp.Abs(1-tortoise)) * (csi - 1)
#csi2 = 1 + (sp.Rational(1,2)-copysign(0.5, 1.5-tortoise)) * (csi - 1)
csi2 = 1 + (sp.Rational(1,2)-copysignresult) * (csi - 1)
#prT = (px*nx + py*ny + pz*nz)*csi2
prT = sp.symbols('prT')
#prTtimesoneminuscsi1inv = prT*(1 - 1/csi1)
prTtimesoneminuscsi1inv = sp.symbols('prTtimesoneminuscsi1inv')
tmppx = px - nx*prTtimesoneminuscsi1inv
tmppy = py - ny*prTtimesoneminuscsi1inv
tmppz = pz - nz*prTtimesoneminuscsi1inv
#pxir = (tmppx*xiUSCOREx + tmppy*xiUSCOREy + tmppz*xiUSCOREz)*r
pxir = sp.symbols('pxir')
pvr = (tmppx*vx + tmppy*vy + tmppz*vz)*r
pvrsq = pvr*pvr
pn = tmppx*nx + tmppy*ny + tmppz*nz
pnsq = pn*pn
pr = pn
prsq = pr*pr
#pf = pxir
pf = sp.symbols('pf')
pxirsq = pxir*pxir
ptheta2 = pvrsq*invxi2
prT4=prT*prT*prT*prT
Hnspt7 = deltaR*invrho2
Hnspt6 = rho2invLambda*invxi2
Hnspt5 = qq*u2
#Hnspt4 = (1 + prT4*Hnspt5 + ptheta2*invrho2 + pf*pf*Hnspt6 + prsq*Hnspt7)
Hnspt4 = sp.symbols('Hnspt4')
Hnspt3 = deltaT*Hnspt4
Hnspt2 = rho2*Hnspt3
Hnspt1 = pf*ww
#Hns = sp.sqrt(Hnspt2*invLambda) + invLambda*Hnspt1
Hns = sp.symbols('Hns')
Qpt3 = deltaR*invrho2
Qpt2 = rho2invLambda*invxi2
Qpt1 = invrho2*invxi2
#Q = 1 + pvrsq*Qpt1 + pxirsq*Qpt2 + pnsq*Qpt3
Q = sp.symbols('Q')
#pn2 = prsq*deltaR*invrho2
pn2 = sp.symbols('pn2')
pp = Q - 1
sKerrmultfact = (-8 - 3*r*(12*pn2 - pp))
sStarmultfact = (14 + (- 30*pn2 + 4*pp)*r)
deltaSigmaStarUSCOREx1=etaover12r*(sKerrmultfact*sKerrUSCOREx + sStarmultfact*sStarUSCOREx)
deltaSigmaStarUSCOREy1=etaover12r*(sKerrmultfact*sKerrUSCOREy + sStarmultfact*sStarUSCOREy)
deltaSigmaStarUSCOREz1=etaover12r*(sKerrmultfact*sKerrUSCOREz + sStarmultfact*sStarUSCOREz)
pn2pp = pn2*pp
pp2 = pp*pp
pn2u2 = pn2*u2
ppu2 = pp*u2
pn2ppu2 = pn2pp*u2
sMultiplier1pt6 = -360*pn2*pn2 + 126*pn2pp + 3*pp2
sMultiplier1pt5 = -96*pn2pp + 23*pp2
sMultiplier1pt4 = -120*pp + 324*pn2 + sMultiplier1pt6*r
sMultiplier1pt3 = 206*pp - 282*pn2 + sMultiplier1pt5*r
sMultiplier1pt2 = 54 + sMultiplier1pt4*r
sMultiplier1pt1 = -706 + sMultiplier1pt3*r + sMultiplier1pt2*eta
#sMultiplier1 = sMultiplier1pt1*eta*u2*sp.Rational(-1,72)
sMultiplier1 = sp.symbols('sMultiplier1')
sMultiplier2pt6 = sp.Rational(45,8)*pn2*pn2u2 - sp.Rational(13,8)*pn2ppu2
sMultiplier2pt5 = pn2ppu2/4 - sp.Rational(5,16)*pp2*u2
#sMultiplier2pt4 = sp.Rational(-49,8)*pn2u2 + sp.Rational(17,12)*ppu2 + sMultiplier2pt6*r
sMultiplier2pt4 = sp.symbols('sMultiplier2pt4')
sMultiplier2pt3 = sp.Rational(-2,3)*pn2u2 - sp.Rational(109,36)*ppu2 + sMultiplier2pt5*r
sMultiplier2pt2 = sp.Rational(-7,3)*u2 + sMultiplier2pt4*r
sMultiplier2pt1 = sp.Rational(-56,9)*u2 + sMultiplier2pt3*r + sMultiplier2pt2*eta
#sMultiplier2 = sMultiplier2pt1*eta
sMultiplier2 = sp.symbols('sMultiplier2')
deltaSigmaStarUSCOREx2 = deltaSigmaStarUSCOREx1 + sMultiplier1*sigmaStar0 + sMultiplier2*sigmaKerr0
deltaSigmaStarUSCOREy2 = deltaSigmaStarUSCOREy1 + sMultiplier1*sigmaStar1 + sMultiplier2*sigmaKerr1
deltaSigmaStarUSCOREz2 = deltaSigmaStarUSCOREz1 + sMultiplier1*sigmaStar2 + sMultiplier2*sigmaKerr2
deltaSigmaStarUSCOREx3 = deltaSigmaStarUSCOREx2 + d1*sigmaStar0*etau3
deltaSigmaStarUSCOREy3 = deltaSigmaStarUSCOREy2 + d1*sigmaStar1*etau3
deltaSigmaStarUSCOREz3 = deltaSigmaStarUSCOREz2 + d1*sigmaStar2*etau3
deltaSigmaStarUSCOREx = deltaSigmaStarUSCOREx3 + d1v2*sigmaKerr0*etau3
deltaSigmaStarUSCOREy = deltaSigmaStarUSCOREy3 + d1v2*sigmaKerr1*etau3
deltaSigmaStarUSCOREz = deltaSigmaStarUSCOREz3 + d1v2*sigmaKerr2*etau3
#sx = sStarUSCOREx + deltaSigmaStarUSCOREx
sx = sp.symbols('sx')
#sy = sStarUSCOREy + deltaSigmaStarUSCOREy
sy = sp.symbols('sy')
#sz = sStarUSCOREz + deltaSigmaStarUSCOREz
sz = sp.symbols('sz')
#sxi = sx*xiUSCOREx + sy*xiUSCOREy + sz*xiUSCOREz
sxi = sp.symbols('sxi')
#sv = sx*vx + sy*vy + sz*vz
sv = sp.symbols('sv')
sn = sx*nx + sy*ny + sz*nz
s3 = sx*e3USCOREx + sy*e3USCOREy + sz*e3USCOREz
sqrtQ = sp.sqrt(Q)
oneplus2sqrtQ = 1 + 2*sqrtQ
oneplus1sqrtQ = oneplus2sqrtQ - sqrtQ
twoB1psqrtQsqrtQ = (2*B*oneplus1sqrtQ*sqrtQ)
invtwoB1psqrtQsqrtQ = 1/twoB1psqrtQsqrtQ
expMUsqsqrtQplusQ = (expMUsq)*(sqrtQ + Q)
Hwrpt4a = pxirsq*sv
#Hwrpt4 = expMUsqexpnusq*Hwrpt4a
Hwrpt4 = sp.symbols('Hwrpt4')
Hwrpt3c = pxir*sxi
#Hwrpt3b = pvr*Hwrpt3c
Hwrpt3b = sp.symbols('Hwrpt3b')
Hwrpt3a = expMUexpnu*Hwrpt3b
Hwrpt3 = B*Hwrpt3a
Hwrpt2g = sv*deltaR
Hwrpt2f = sn*sqrtdeltaR
Hwrpt2e = pvr*Hwrpt2f
Hwrpt2d = pnsq*Hwrpt2g
Hwrpt2c = pn*Hwrpt2e
Hwrpt2b = expMUsqsqrtQplusQ*sv
Hwrpt2a = xi2*(Hwrpt2b + Hwrpt2c - Hwrpt2d)
#Hwrpt2 = deltaT*Hwrpt2a
Hwrpt2 = sp.symbols('Hwrpt2')
Hwrpt1b = invtwoB1psqrtQsqrtQ*invxi2
Hwrpt1a = sqrtdeltaR*Hwrpt1b
#Hwrpt1 = invexpMUcubinvexpnu*Hwrpt1a
Hwrpt1 = sp.symbols('Hwrpt1')
#Hwr = (Hwrpt4 - Hwrpt3 + Hwrpt2)*Hwrpt1
Hwr = sp.symbols('Hwr')
Hwcospt9 = pxir*sxi
Hwcospt8 = pvr*sv
Hwcospt7 = (B*Hwcospt8 - (expMUexpnu)*Hwcospt9)
Hwcospt6 = sqrtdeltaR*Hwcospt7
#Hwcospt5 = (pvrsq - expMUsqsqrtQplusQ*xi2)
Hwcospt5 = sp.symbols('Hwcospt5')
Hwcospt4 = pn*Hwcospt6
#Hwcospt3 = -(expMUsqexpnusq*pxirsq) + deltaT*Hwcospt5
Hwcospt3 = sp.symbols('Hwcospt3')
Hwcospt2 = sn*Hwcospt3 - B*Hwcospt4
Hwcospt1 = invexpMUcubinvexpnu*Hwcospt2
#Hwcos = invtwoB1psqrtQsqrtQ*Hwcospt1
Hwcos = sp.symbols('Hwcos')
deltaTsqrtQ = deltaT*sqrtQ
invdeltatTsqrtQ = 1/deltaTsqrtQ
#HSOLpt5 = (-B + (expMUexpnu))*pxir
HSOLpt5 = sp.symbols('HSOLpt5')
HSOLpt4 = invexpMU*HSOLpt5
HSOLpt3 = expnusq*HSOLpt4
HSOLpt2 = (HSOLpt3*s3)
HSOLpt1 = HSOLpt2*invxi2
#HSOL = HSOLpt1*invdeltatTsqrtQ
HSOL = sp.symbols('HSOL')
deltaTsqrtQplusQ = (deltaT*(sqrtQ + Q))
invdeltaTsqrtQplusQ = 1/deltaTsqrtQplusQ
#HSONLmult2 = invxi2*invdeltaTsqrtQplusQ
HSONLmult2 = sp.symbols('HSONLmult2')
#HSONLmult = expnuinvexpMU2*HSONLmult2
HSONLmult = sp.symbols('HSONLmult')
HSONLpt1b = pn*xi2
HSONLpt1a = (mur*pvr - nur*pvr + (-mucos + nucos)*HSONLpt1b)
#HSONLpt1 = mur*pvr - (mucos*HSONLpt1b) + sqrtQ*HSONLpt1a
HSONLpt1 = sp.symbols('HSONLpt1')
HSONLpt2d = nur*pxir
HSONLpt2c = oneplus2sqrtQ*HSONLpt2d
HSONLpt2b = B*sxi
HSONLpt2a = expMUexpnu*HSONLpt2c
#HSONLpt2 = (sv*HSONLpt2a + HSONLpt1*HSONLpt2b)
HSONLpt2 = sp.symbols('HSONLpt2')
HSONLpt3c = sv*pxir
HSONLpt3b = oneplus1sqrtQ*HSONLpt3c
HSONLpt3a = expMUexpnu*HSONLpt3b
HSONLpt3 = -BR*HSONLpt3a + B*HSONLpt2
HSONLpt4e = sn*xi2
HSONLpt4d = oneplus2sqrtQ*HSONLpt4e
HSONLpt4c = pxir*HSONLpt4d
HSONLpt4b = nucos*HSONLpt4c
HSONLpt4a = expMUexpnu*HSONLpt4b
#HSONLpt4 = (-(B*HSONLpt4a) + HSONLpt3*sqrtdeltaR)
HSONLpt4 = sp.symbols('HSONLpt4')
HSONL = HSONLmult*HSONLpt4
#Hs = w*s3 + Hwr*wr + Hwcos*wcos + HSOL + HSONL
Hs = sp.symbols('Hs')
Hsspt1 = sp.Rational(-1,2)*(sx*sx + sy*sy + sz*sz - 3*sn*sn)
#Hss = u3*Hsspt1
Hss = sp.symbols('Hss')
sKerrdotsStar = (sKerrUSCOREx*sStarUSCOREx + sKerrUSCOREy*sStarUSCOREy + sKerrUSCOREz*sStarUSCOREz)
Hpt1 = etau4*(s1dots1 + s2dots2)
H = Hns + Hs + Hss + (dheffSS*sKerrdotsStar + dheffSSv2)*Hpt1
Hreal = sp.sqrt(1 + 2*eta*(H - 1))
#print(sp.mathematica_code(expnu))
from outputC import *
import time
start = time.time()
#outstring = "test"
outstring = outputC(Hreal,"Hrealtest","returnstring","outCverbose=False,enable_SIMD=True,enable_TYPE=False,includebraces=False")
end = time.time()
print("Finished in " + str(end - start) + " seconds.")
Finished in 0.61027598381 seconds.
print("""(* Part 1: SIMD function definitions *)
MulSIMD[x_, y_] = x*y;
DivSIMD[x_, y_] = x/y;
AddSIMD[x_, y_] = x + y;
FusedMulAddSIMD[x_, y_, z_] = x*y + z;
FusedMulSubSIMD[x_, y_, z_] = x*y - z;
SqrtSIMD[x_] = Sqrt[x];
LogSIMD[x_] = Log[x];
SubSIMD[x_, y_] = x - y;
AbsSIMD[x_] = Abs[x];
PowSIMD[x_, y_] = x^y;
""")
print("(* Part 2: SEOBNR Hamiltonian *)")
print(outstring.replace("_Rational","Rational").replace("_Integer","Integer").replace("_","USCR")).replace("(","[").replace(")","]").replace(".000000000000000000000000000000000","").replace(".00000000000000000000000000000000","").replace("000000000000000000000000000000000","").replace("00000000000000000000000000000000","").replace(".0000000000000000000000000000000","").replace("000000000000000000000000000000","")
# replace(";","").
(* Part 1: SIMD function definitions *) MulSIMD[x_, y_] = x*y; DivSIMD[x_, y_] = x/y; AddSIMD[x_, y_] = x + y; FusedMulAddSIMD[x_, y_, z_] = x*y + z; FusedMulSubSIMD[x_, y_, z_] = x*y - z; SqrtSIMD[x_] = Sqrt[x]; LogSIMD[x_] = Log[x]; SubSIMD[x_, y_] = x - y; AbsSIMD[x_] = Abs[x]; PowSIMD[x_, y_] = x^y; (* Part 2: SEOBNR Hamiltonian *) IntegerUSCR1 = 1; IntegerUSCR2 = 2; IntegerUSCRm1 = -1; tmp0 = DivSIMD[m2, m1]; tmp1 = DivSIMD[m1, m2]; Hrealtest = SqrtSIMD[FusedMulAddSIMD[AddSIMD[IntegerUSCRm1, AddSIMD[Hss, AddSIMD[Hs, FusedMulAddSIMD[FusedMulAddSIMD[s2y, s2y, FusedMulAddSIMD[s2x, s2x, FusedMulAddSIMD[s1z, s1z, FusedMulAddSIMD[s1y, s1y, FusedMulAddSIMD[s2z, s2z, MulSIMD[s1x, s1x]]]]]], MulSIMD[FusedMulAddSIMD[dheffSS, FusedMulAddSIMD[AddSIMD[s1y, s2y], FusedMulAddSIMD[tmp1, s2y, MulSIMD[tmp0, s1y]], FusedMulAddSIMD[AddSIMD[s1z, s2z], FusedMulAddSIMD[tmp1, s2z, MulSIMD[tmp0, s1z]], MulSIMD[AddSIMD[s1x, s2x], FusedMulAddSIMD[tmp1, s2x, MulSIMD[tmp0, s1x]]]]], dheffSSv2], DivSIMD[eta, MulSIMD[FusedMulAddSIMD[y, y, FusedMulAddSIMD[z, z, MulSIMD[x, x]]], FusedMulAddSIMD[y, y, FusedMulAddSIMD[z, z, MulSIMD[x, x]]]]]], Hns]]]], MulSIMD[eta, IntegerUSCR2], IntegerUSCR1]];
bigstring = """
sigmaKerr0 = s1x + s2x
sigmaKerr1 = s1y + s2y
sigmaKerr2 = s1z + s2z
invm1m2 = 1/(m1*m2)
m2overm1 = m2*m2*invm1m2
m1overm2 = m1*m1*invm1m2
sigmaStar0 = (m2overm1)*s1x + (m1overm2)*s2x
sigmaStar1 = (m2overm1)*s1y + (m1overm2)*s2y
sigmaStar2 = (m2overm1)*s1z + (m1overm2)*s2z
s1dots1 = s1x*s1x + s1y*s1y + s1z*s1z
s2dots2 = s2x*s2x + s2y*s2y + s2z*s2z
r2 = x*x + y*y + z*z
r = sp.sqrt(r2)
u = 1/r
u2 = u*u
u3 = u2*u
u4 = u2*u2
u5 = u4*u
etau3 = eta*u3
etau4 = eta*u4
nx = x*u
ny = y*u
nz = z*u
sKerrUSCOREx = sigmaKerr0
sKerrUSCOREy = sigmaKerr1
sKerrUSCOREz = sigmaKerr2
sStarUSCOREx = sigmaStar0
sStarUSCOREy = sigmaStar1
sStarUSCOREz = sigmaStar2
a2 = sKerrUSCOREx*sKerrUSCOREx + sKerrUSCOREy*sKerrUSCOREy + sKerrUSCOREz*sKerrUSCOREz
a4 = a2*a2
a = sp.sqrt(a2)
inva = 1/a
m1PlusetaKK = -1 + eta*KK
invm1PlusetaKK = 1/m1PlusetaKK
k0 = KK*(m1PlusetaKK - 1)
k1 = -2*(k0 + KK)*m1PlusetaKK
k2 = c0k2 + c1k2*a2
k3 = c0k3 + c1k3*a2
k4 = c0k4 + c1k4*a2 + c2k4*a4
k5 = c0k5 + c1k5*a2 + c2k5*a4
e3USCOREx = sKerrUSCOREx*inva
e3USCOREy = sKerrUSCOREy*inva
e3USCOREz = sKerrUSCOREz*inva
costheta = e3USCOREx*nx + e3USCOREy*ny + e3USCOREz*nz
xi2 = 1 - costheta*costheta
xiUSCOREx = -e3USCOREz*ny + e3USCOREy*nz
xiUSCOREy = e3USCOREz*nx - e3USCOREx*nz
xiUSCOREz = -e3USCOREy*nx + e3USCOREx*ny
vx = -nz*xiUSCOREy + ny*xiUSCOREz
vy = nz*xiUSCOREx - nx*xiUSCOREz
vz = -ny*xiUSCOREx + nx*xiUSCOREy
w2 = r2 + a2
rho2 = r2 + a2*costheta*costheta
bulk = invm1PlusetaKK*(invm1PlusetaKK + 2*u) + a2*u2
logu = sp.log( u )
logarg = k1*u + k2*u2 + k3*u3 + k4*u4 + k5*u5 + k5l*u5*logu
onepluslogarg = (1 + logarg)
invonepluslogarg = 1/onepluslogarg
logTerms = 1 + eta*k0 + eta*sp.log(onepluslogarg)
deltaU = bulk*logTerms
deltaT = r2*deltaU
deltaUUSCOREupt7 = k5 + k5l*logu
deltaUUSCOREupt6 = 4*k4 + 5*deltaUUSCOREupt7*u
deltaUUSCOREupt5 = 3*k3 + u*deltaUUSCOREupt6
deltaUUSCOREupt4 = 2*k2 + u*deltaUUSCOREupt5
deltaUUSCOREupt3 = k1 + u*deltaUUSCOREupt4
deltaUUSCOREupt2 = invm1PlusetaKK + a2*u
deltaUUSCOREupt1 = bulk*eta*deltaUUSCOREupt3
deltaUUSCOREu = 2*deltaUUSCOREupt2*logTerms + deltaUUSCOREupt1*invonepluslogarg
deltaTUSCOREr = 2*r*deltaU - deltaUUSCOREu
Lambda = w2*w2 - a2*deltaT*xi2
rho2xi2Lambda = rho2*xi2*Lambda
invrho2xi2Lambda = 1/rho2xi2Lambda
invrho2 = xi2*Lambda*invrho2xi2Lambda
invxi2 = rho2*Lambda*invrho2xi2Lambda
invLambda = xi2*rho2*invrho2xi2Lambda
invLambdasq = invLambda*invLambda
rho2invLambda = rho2*invLambda
expnu = sp.sqrt(deltaT*rho2invLambda)
expMU = sp.sqrt(rho2)
expMUexpnu = expMU*expnu
expMUsq = expMU*expMU
expnusq = expnu*expnu
expMUsqexpnusq = expMUsq*expnusq
invexpnuexpMU = 1/expMUexpnu
invexpMU = expnu*invexpnuexpMU
invexpMUsq = invexpMU*invexpMU
expnuinvexpMU2 = expnu*invexpMUsq
invexpMUcubinvexpnu = invexpMUsq*invexpnuexpMU
DD = 1 + sp.log(1 + 6*eta*u2 + 2*(26 - 3*eta)*etau3)
deltaR = deltaT*DD
qq = 2*eta*(4 - 3*eta)
ww = 2*a*r + b3*eta*a2*a*u + bb3*eta*a*u
B = sp.sqrt(deltaT)
sqrtdeltaT = B
sqrtdeltaR = sp.sqrt(deltaR)
deltaTsqrtdeltaR = deltaT*sqrtdeltaR
sqrtdeltaTdeltaTsqrtdeltaR = sqrtdeltaT*deltaTsqrtdeltaR
invdeltaTsqrtdeltaTsqrtdeltaR = 1./sqrtdeltaTdeltaTsqrtdeltaR
invdeltaT = sqrtdeltaT*(sqrtdeltaR*invdeltaTsqrtdeltaTsqrtdeltaR)
invsqrtdeltaT = deltaTsqrtdeltaR*invdeltaTsqrtdeltaTsqrtdeltaR
invsqrtdeltaR = deltaT*sqrtdeltaT*invdeltaTsqrtdeltaTsqrtdeltaR
w = ww*invLambda
LambdaUSCOREr = 4*r*w2 - a2*deltaTUSCOREr*xi2
wwUSCOREr = 2*a - (a2*a*b3*eta)*u2 - bb3*eta*a*u2
BR = (-deltaT*invsqrtdeltaR + deltaTUSCOREr*sp.Rational(1,2))*invsqrtdeltaT
wr = (-LambdaUSCOREr*ww + Lambda*wwUSCOREr)*(invLambdasq)
nurpt2 = w2*(-4.*r*deltaT + w2*deltaTUSCOREr)
nurpt1 = nurpt2*invdeltaT
nur = r*invrho2 + sp.Rational(1,2)*invLambda*nurpt1
mur = r*invrho2 - invsqrtdeltaR
a2costheta = a2*costheta
wcospt2 = deltaT*ww
wcospt1 = invLambdasq*wcospt2
wcos = -2*(a2costheta)*wcospt1
nucospt3 = invrho2*invLambda
nucospt2 = w2*nucospt3
nucospt1 = a2costheta*nucospt2
nucos = (w2 - deltaT)*nucospt1
mucos = a2costheta*invrho2
etaover12r = eta*sp.Rational(1,12)*u
csi = sp.sqrt(deltaT*deltaR)/w2
csi1 = 1 + (1-sp.Abs(1-tortoise)) * (csi - 1)
#csi2 = 1 + (sp.Rational(1,2)-copysign(0.5, 1.5-tortoise)) * (csi - 1)
csi2 = 1 + (sp.Rational(1,2)-copysignresult) * (csi - 1)
prT = (px*nx + py*ny + pz*nz)*csi2
prTtimesoneminuscsi1inv = prT*(1 - 1/csi1)
tmppx = px - nx*prTtimesoneminuscsi1inv
tmppy = py - ny*prTtimesoneminuscsi1inv
tmppz = pz - nz*prTtimesoneminuscsi1inv
pxir = (tmppx*xiUSCOREx + tmppy*xiUSCOREy + tmppz*xiUSCOREz)*r
pvr = (tmppx*vx + tmppy*vy + tmppz*vz)*r
pvrsq = pvr*pvr
pn = tmppx*nx + tmppy*ny + tmppz*nz
pnsq = pn*pn
pr = pn
prsq = pr*pr
pf = pxir
pxirsq = pxir*pxir
ptheta2 = pvrsq*invxi2
prT4=prT*prT*prT*prT
Hnspt7 = deltaR*invrho2
Hnspt6 = rho2invLambda*invxi2
Hnspt5 = qq*u2
Hnspt4 = (1 + prT4*Hnspt5 + ptheta2*invrho2 + pf*pf*Hnspt6 + prsq*Hnspt7)
Hnspt3 = deltaT*Hnspt4
Hnspt2 = rho2*Hnspt3
Hnspt1 = pf*ww
Hns = sp.sqrt(Hnspt2*invLambda) + invLambda*Hnspt1
Qpt3 = deltaR*invrho2
Qpt2 = rho2invLambda*invxi2
Qpt1 = invrho2*invxi2
Q = 1 + pvrsq*Qpt1 + pxirsq*Qpt2 + pnsq*Qpt3
pn2 = prsq*deltaR*invrho2
pp = Q - 1
sKerrmultfact = (-8 - 3*r*(12*pn2 - pp))
sStarmultfact = (14 + (- 30*pn2 + 4*pp)*r)
deltaSigmaStarUSCOREx1=etaover12r*(sKerrmultfact*sKerrUSCOREx + sStarmultfact*sStarUSCOREx)
deltaSigmaStarUSCOREy1=etaover12r*(sKerrmultfact*sKerrUSCOREy + sStarmultfact*sStarUSCOREy)
deltaSigmaStarUSCOREz1=etaover12r*(sKerrmultfact*sKerrUSCOREz + sStarmultfact*sStarUSCOREz)
pn2pp = pn2*pp
pp2 = pp*pp
pn2u2 = pn2*u2
ppu2 = pp*u2
pn2ppu2 = pn2pp*u2
sMultiplier1pt6 = -360*pn2*pn2 + 126*pn2pp + 3*pp2
sMultiplier1pt5 = -96*pn2pp + 23*pp2
sMultiplier1pt4 = -120*pp + 324*pn2 + sMultiplier1pt6*r
sMultiplier1pt3 = 206*pp - 282*pn2 + sMultiplier1pt5*r
sMultiplier1pt2 = 54 + sMultiplier1pt4*r
sMultiplier1pt1 = -706 + sMultiplier1pt3*r + sMultiplier1pt2*eta
sMultiplier1 = sMultiplier1pt1*eta*u2*sp.Rational(-1,72)
sMultiplier2pt6 = sp.Rational(45,8)*pn2*pn2u2 - sp.Rational(13,8)*pn2ppu2
sMultiplier2pt5 = pn2ppu2/4 - sp.Rational(5,16)*pp2*u2
sMultiplier2pt4 = sp.Rational(-49,8)*pn2u2 + sp.Rational(17,12)*ppu2 + sMultiplier2pt6*r
sMultiplier2pt3 = sp.Rational(-2,3)*pn2u2 - sp.Rational(109,36)*ppu2 + sMultiplier2pt5*r
sMultiplier2pt2 = sp.Rational(-7,3)*u2 + sMultiplier2pt4*r
sMultiplier2pt1 = sp.Rational(-56,9)*u2 + sMultiplier2pt3*r + sMultiplier2pt2*eta
sMultiplier2 = sMultiplier2pt1*eta
deltaSigmaStarUSCOREx2 = deltaSigmaStarUSCOREx1 + sMultiplier1*sigmaStar0 + sMultiplier2*sigmaKerr0
deltaSigmaStarUSCOREy2 = deltaSigmaStarUSCOREy1 + sMultiplier1*sigmaStar1 + sMultiplier2*sigmaKerr1
deltaSigmaStarUSCOREz2 = deltaSigmaStarUSCOREz1 + sMultiplier1*sigmaStar2 + sMultiplier2*sigmaKerr2
deltaSigmaStarUSCOREx3 = deltaSigmaStarUSCOREx2 + d1*sigmaStar0*etau3
deltaSigmaStarUSCOREy3 = deltaSigmaStarUSCOREy2 + d1*sigmaStar1*etau3
deltaSigmaStarUSCOREz3 = deltaSigmaStarUSCOREz2 + d1*sigmaStar2*etau3
deltaSigmaStarUSCOREx = deltaSigmaStarUSCOREx3 + d1v2*sigmaKerr0*etau3
deltaSigmaStarUSCOREy = deltaSigmaStarUSCOREy3 + d1v2*sigmaKerr1*etau3
deltaSigmaStarUSCOREz = deltaSigmaStarUSCOREz3 + d1v2*sigmaKerr2*etau3
sx = sStarUSCOREx + deltaSigmaStarUSCOREx
sy = sStarUSCOREy + deltaSigmaStarUSCOREy
sz = sStarUSCOREz + deltaSigmaStarUSCOREz
sxi = sx*xiUSCOREx + sy*xiUSCOREy + sz*xiUSCOREz
sv = sx*vx + sy*vy + sz*vz
sn = sx*nx + sy*ny + sz*nz
s3 = sx*e3USCOREx + sy*e3USCOREy + sz*e3USCOREz
sqrtQ = sp.sqrt(Q)
oneplus2sqrtQ = 1 + 2*sqrtQ
oneplus1sqrtQ = oneplus2sqrtQ - sqrtQ
twoB1psqrtQsqrtQ = (2*B*oneplus1sqrtQ*sqrtQ)
invtwoB1psqrtQsqrtQ = 1/twoB1psqrtQsqrtQ
expMUsqsqrtQplusQ = (expMUsq)*(sqrtQ + Q)
Hwrpt4a = pxirsq*sv
Hwrpt4 = expMUsqexpnusq*Hwrpt4a
Hwrpt3c = pxir*sxi
Hwrpt3b = pvr*Hwrpt3c
Hwrpt3a = expMUexpnu*Hwrpt3b
Hwrpt3 = B*Hwrpt3a
Hwrpt2g = sv*deltaR
Hwrpt2f = sn*sqrtdeltaR
Hwrpt2e = pvr*Hwrpt2f
Hwrpt2d = pnsq*Hwrpt2g
Hwrpt2c = pn*Hwrpt2e
Hwrpt2b = expMUsqsqrtQplusQ*sv
Hwrpt2a = xi2*(Hwrpt2b + Hwrpt2c - Hwrpt2d)
Hwrpt2 = deltaT*Hwrpt2a
Hwrpt1b = invtwoB1psqrtQsqrtQ*invxi2
Hwrpt1a = sqrtdeltaR*Hwrpt1b
Hwrpt1 = invexpMUcubinvexpnu*Hwrpt1a
Hwr = (Hwrpt4 - Hwrpt3 + Hwrpt2)*Hwrpt1
Hwcospt9 = pxir*sxi
Hwcospt8 = pvr*sv
Hwcospt7 = (B*Hwcospt8 - (expMUexpnu)*Hwcospt9)
Hwcospt6 = sqrtdeltaR*Hwcospt7
Hwcospt5 = (pvrsq - expMUsqsqrtQplusQ*xi2)
Hwcospt4 = pn*Hwcospt6
Hwcospt3 = -(expMUsqexpnusq*pxirsq) + deltaT*Hwcospt5
Hwcospt2 = sn*Hwcospt3 - B*Hwcospt4
Hwcospt1 = invexpMUcubinvexpnu*Hwcospt2
Hwcos = invtwoB1psqrtQsqrtQ*Hwcospt1
deltaTsqrtQ = deltaT*sqrtQ
invdeltatTsqrtQ = 1/deltaTsqrtQ
HSOLpt5 = (-B + (expMUexpnu))*pxir
HSOLpt4 = invexpMU*HSOLpt5
HSOLpt3 = expnusq*HSOLpt4
HSOLpt2 = (HSOLpt3*s3)
HSOLpt1 = HSOLpt2*invxi2
HSOL = HSOLpt1*invdeltatTsqrtQ
deltaTsqrtQplusQ = (deltaT*(sqrtQ + Q))
invdeltaTsqrtQplusQ = 1/deltaTsqrtQplusQ
HSONLmult2 = invxi2*invdeltaTsqrtQplusQ
HSONLmult = expnuinvexpMU2*HSONLmult2
HSONLpt1b = pn*xi2
HSONLpt1a = (mur*pvr - nur*pvr + (-mucos + nucos)*HSONLpt1b)
HSONLpt1 = mur*pvr - (mucos*HSONLpt1b) + sqrtQ*HSONLpt1a
HSONLpt2d = nur*pxir
HSONLpt2c = oneplus2sqrtQ*HSONLpt2d
HSONLpt2b = B*sxi
HSONLpt2a = expMUexpnu*HSONLpt2c
HSONLpt2 = (sv*HSONLpt2a + HSONLpt1*HSONLpt2b)
HSONLpt3c = sv*pxir
HSONLpt3b = oneplus1sqrtQ*HSONLpt3c
HSONLpt3a = expMUexpnu*HSONLpt3b
HSONLpt3 = -BR*HSONLpt3a + B*HSONLpt2
HSONLpt4e = sn*xi2
HSONLpt4d = oneplus2sqrtQ*HSONLpt4e
HSONLpt4c = pxir*HSONLpt4d
HSONLpt4b = nucos*HSONLpt4c
HSONLpt4a = expMUexpnu*HSONLpt4b
HSONLpt4 = (-(B*HSONLpt4a) + HSONLpt3*sqrtdeltaR)
HSONL = HSONLmult*HSONLpt4
Hs = w*s3 + Hwr*wr + Hwcos*wcos + HSOL + HSONL
Hsspt1 = sp.Rational(-1,2)*(sx*sx + sy*sy + sz*sz - 3*sn*sn)
Hss = u3*Hsspt1
sKerrdotsStar = (sKerrUSCOREx*sStarUSCOREx + sKerrUSCOREy*sStarUSCOREy + sKerrUSCOREz*sStarUSCOREz)
Hpt1 = etau4*(s1dots1 + s2dots2)
H = Hns + Hs + Hss + (dheffSS*sKerrdotsStar + dheffSSv2)*Hpt1
Hreal = sp.sqrt(1 + 2*eta*(H - 1))
"""
from outputC import *
# Split bigstring by carriage returns:
blah = bigstring.splitlines()
# Create "lr" array, which will store each left-hand side and right-hand side as strings.
lr = []
# Loop over each line in bigstring
for i in range(len(blah)):
# Ignore lines with 2 or fewer characters and those starting with #
if len(blah[i]) > 2 and blah[i][0] != "#":
# Split each line by its equals sign.
splitblah = blah[i].split("=")
# Append to the "lr" array, removing spaces, "sp." prefixes, and replacing Lambda->Lamb
# (Lambda is a protected keyword):
#lhs =
lr.append(lhrh(lhs=splitblah[0].replace(" ","").replace("Lambda","Lamb"),
rhs=splitblah[1].replace(" ","").replace("sp.","").replace("Lambda","Lamb")))
# print(lr[1].lhs)
xx = sp.Symbol('xx')
func = []
lhss = []
rhss = []
for i in range(len(lr)):
func.append(sp.sympify(sp.Function(lr[i].lhs)(xx)))
# print(i,lr[i].rhs)
lhss.append(sp.sympify(lr[i].lhs))
rhss.append(sp.sympify(lr[i].rhs))
# Next get a list of all the "free symbols" in the RHS expressions.
full_symbol_list_with_dups = []
for i in range(len(lr)):
for var in rhss[i].free_symbols:
full_symbol_list_with_dups.append(var)
full_symbol_list = superfast_uniq(full_symbol_list_with_dups)
# Next declare independent variables:
# x,y,z,px,py,pz,s1x,s1y,s1z,s2x,s2y,s2z = sp.symbols("x y z px py pz s1x s1y s1z s2x s2y s2z",real=True)
# independent_variables = [x,y,z,px,py,pz,s1x,s1y,s1z,s2x,s2y,s2z]
# Next declare input constants:
m1,m2,eta = sp.symbols("m1 m2 eta",real=True)
c0k2,c1k2,c0k3,c1k3,c0k4,c1k4,c2k4,c0k5,c1k5,c2k5 = sp.symbols("c0k2 c1k2 c0k3 c1k3 c0k4 c1k4 c2k4 c0k5 c1k5 c2k5",real=True)
KK,k5l,b3,bb3,d1,d1v2,dheffSS,dheffSSv2 = sp.symbols("KK k5l b3 bb3 d1 d1v2 dheffSS dheffSSv2",real=True)
tortoise,copysignresult = sp.symbols("tortoise copysignresult",real=True)
input_constants = [m1,m2,eta,
c0k2,c1k2,c0k3,c1k3,c0k4,c1k4,c2k4,c0k5,c1k5,c2k5,
KK,k5l,b3,bb3,d1,d1v2,dheffSS,dheffSSv2,
tortoise,copysignresult]
# Derivatives of input constants will always be zero, so
# print(full_symbol_list)
for inputconst in input_constants:
for symbol in full_symbol_list:
if str(symbol) == str(inputconst):
full_symbol_list.remove(symbol)
full_function_list = []
for symb in full_symbol_list:
func = sp.sympify(sp.Function(str(symb))(xx))
full_function_list.append(func)
for i in range(len(rhss)):
for var in rhss[i].free_symbols:
if str(var) == str(symb):
rhss[i] = rhss[i].subs(var,func)
lhss_deriv = []
rhss_deriv = []
for i in range(len(rhss)):
lhss_deriv.append(sp.sympify(str(lhss[i])+"prm"))
newrhs = sp.sympify(str(sp.diff(rhss[i],xx)).replace("(xx)","").replace(", xx","prm").replace("Derivative",""))
rhss_deriv.append(newrhs)
# rhss_deriv.append(sp.diff(rhss[i],xx))
def simplify_deriv(lhss_deriv,rhss_deriv):
lhss_deriv_simp = []
rhss_deriv_simp = []
for i in range(len(rhss_deriv)):
lhss_deriv_simp.append(lhss_deriv[i])
rhss_deriv_simp.append(rhss_deriv[i])
for i in range(len(rhss_deriv_simp)):
if rhss_deriv_simp[i] == 0:
for j in range(i+1,len(rhss_deriv_simp)):
for var in rhss_deriv_simp[j].free_symbols:
if str(var) == str(lhss_deriv_simp[i]):
rhss_deriv_simp[j] = rhss_deriv_simp[j].subs(var,0)
zero_elements_to_remove = []
for i in range(len(rhss_deriv_simp)):
if rhss_deriv_simp[i] == sp.sympify(0):
zero_elements_to_remove.append(i)
count = 0
for i in range(len(zero_elements_to_remove)):
del lhss_deriv_simp[zero_elements_to_remove[i]+count]
del rhss_deriv_simp[zero_elements_to_remove[i]+count]
count -= 1
return lhss_deriv_simp,rhss_deriv_simp
lhss_deriv_simp,rhss_deriv_simp = simplify_deriv(lhss_deriv,rhss_deriv)
lhss_deriv = lhss_deriv_simp
rhss_deriv = rhss_deriv_simp
def deriv_onevar(lhss_deriv,rhss_deriv,
xprm=0,yprm=0,zprm=0,pxprm=0,pyprm=0,pzprm=0,s1xprm=0,s1yprm=0,s1zprm=0,s2xprm=0,s2yprm=0,s2zprm=0):
lhss_deriv_new = []
rhss_deriv_new = []
for i in range(len(rhss_deriv)):
lhss_deriv_new.append(lhss_deriv[i])
rhss_deriv_new.append(rhss_deriv[i])
for i in range(len(rhss_deriv_new)):
for var in rhss_deriv_new[i].free_symbols:
if str(var)=="xprm":
rhss_deriv_new[i] = rhss_deriv_new[i].subs(var,xprm)
elif str(var)=="yprm":
rhss_deriv_new[i] = rhss_deriv_new[i].subs(var,yprm)
elif str(var)=="zprm":
rhss_deriv_new[i] = rhss_deriv_new[i].subs(var,zprm)
elif str(var)=="pxprm":
rhss_deriv_new[i] = rhss_deriv_new[i].subs(var,pxprm)
elif str(var)=="pyprm":
rhss_deriv_new[i] = rhss_deriv_new[i].subs(var,pyprm)
elif str(var)=="pzprm":
rhss_deriv_new[i] = rhss_deriv_new[i].subs(var,pzprm)
elif str(var)=="s1xprm":
rhss_deriv_new[i] = rhss_deriv_new[i].subs(var,s1xprm)
elif str(var)=="s1yprm":
rhss_deriv_new[i] = rhss_deriv_new[i].subs(var,s1yprm)
elif str(var)=="s1zprm":
rhss_deriv_new[i] = rhss_deriv_new[i].subs(var,s1zprm)
elif str(var)=="s2xprm":
rhss_deriv_new[i] = rhss_deriv_new[i].subs(var,s2xprm)
elif str(var)=="s2yprm":
rhss_deriv_new[i] = rhss_deriv_new[i].subs(var,s2yprm)
elif str(var)=="s2zprm":
rhss_deriv_new[i] = rhss_deriv_new[i].subs(var,s2zprm)
lhss_deriv_simp,rhss_deriv_simp = simplify_deriv(lhss_deriv_new,rhss_deriv_new)
return lhss_deriv_simp,rhss_deriv_simp
lhss_deriv_x,rhss_deriv_x = deriv_onevar(lhss_deriv,rhss_deriv, \
xprm=1,yprm=0,zprm=0,pxprm=0,pyprm=0,pzprm=0,s1xprm=0,s1yprm=0,s1zprm=0,s2xprm=0,s2yprm=0,s2zprm=0)
lhss_deriv_y,rhss_deriv_y = deriv_onevar(lhss_deriv,rhss_deriv,
xprm=0,yprm=1,zprm=0,pxprm=0,pyprm=0,pzprm=0,s1xprm=0,s1yprm=0,s1zprm=0,s2xprm=0,s2yprm=0,s2zprm=0)
lhss_deriv_z,rhss_deriv_z = deriv_onevar(lhss_deriv,rhss_deriv,
xprm=0,yprm=0,zprm=1,pxprm=0,pyprm=0,pzprm=0,s1xprm=0,s1yprm=0,s1zprm=0,s2xprm=0,s2yprm=0,s2zprm=0)
lhss_deriv_px,rhss_deriv_px = deriv_onevar(lhss_deriv,rhss_deriv,
xprm=0,yprm=0,zprm=0,pxprm=1,pyprm=0,pzprm=0,s1xprm=0,s1yprm=0,s1zprm=0,s2xprm=0,s2yprm=0,s2zprm=0)
lhss_deriv_py,rhss_deriv_py = deriv_onevar(lhss_deriv,rhss_deriv,
xprm=0,yprm=0,zprm=0,pxprm=0,pyprm=1,pzprm=0,s1xprm=0,s1yprm=0,s1zprm=0,s2xprm=0,s2yprm=0,s2zprm=0)
lhss_deriv_pz,rhss_deriv_pz = deriv_onevar(lhss_deriv,rhss_deriv,
xprm=0,yprm=0,zprm=0,pxprm=0,pyprm=1,pzprm=1,s1xprm=0,s1yprm=0,s1zprm=0,s2xprm=0,s2yprm=0,s2zprm=0)
lhss_deriv_s1x,rhss_deriv_s1x = deriv_onevar(lhss_deriv,rhss_deriv,
xprm=0,yprm=0,zprm=0,pxprm=0,pyprm=0,pzprm=0,s1xprm=1,s1yprm=0,s1zprm=0,s2xprm=0,s2yprm=0,s2zprm=0)
lhss_deriv_s1y,rhss_deriv_s1y = deriv_onevar(lhss_deriv,rhss_deriv,
xprm=0,yprm=0,zprm=0,pxprm=0,pyprm=0,pzprm=0,s1xprm=0,s1yprm=1,s1zprm=0,s2xprm=0,s2yprm=0,s2zprm=0)
lhss_deriv_s1z,rhss_deriv_s1z = deriv_onevar(lhss_deriv,rhss_deriv,
xprm=0,yprm=0,zprm=0,pxprm=0,pyprm=0,pzprm=0,s1xprm=0,s1yprm=0,s1zprm=1,s2xprm=0,s2yprm=0,s2zprm=0)
lhss_deriv_s2x,rhss_deriv_s2x = deriv_onevar(lhss_deriv,rhss_deriv,
xprm=0,yprm=0,zprm=0,pxprm=0,pyprm=0,pzprm=0,s1xprm=0,s1yprm=0,s1zprm=0,s2xprm=1,s2yprm=0,s2zprm=0)
lhss_deriv_s2y,rhss_deriv_s2y = deriv_onevar(lhss_deriv,rhss_deriv,
xprm=0,yprm=0,zprm=0,pxprm=0,pyprm=0,pzprm=0,s1xprm=0,s1yprm=0,s1zprm=0,s2xprm=0,s2yprm=1,s2zprm=0)
lhss_deriv_s2z,rhss_deriv_s2z = deriv_onevar(lhss_deriv,rhss_deriv,
xprm=0,yprm=0,zprm=0,pxprm=0,pyprm=0,pzprm=0,s1xprm=0,s1yprm=0,s1zprm=0,s2xprm=0,s2yprm=0,s2zprm=1)
# CSE_results = sp.cse(rhss_deriv_px, sp.numbered_symbols("tmp"), order='canonical')
# for commonsubexpression in CSE_results[0]:
# print(" "+str(commonsubexpression[0])+" = "+str(commonsubexpression[1]))
# for i,result in enumerate(CSE_results[1]):
# print("rhss_deriv_px = "+str(result))
#for commonsubexpression in CSE_results[0]:
outstring = "/* SEOBNR Hamiltonian expression: */\n"
outstringsp = ""
outsplhs = []
outsprhs = []
for i in range(len(lr)):
outstring += outputC(sp.sympify(lr[i].rhs),lr[i].lhs,"returnstring","outCverbose=False,includebraces=False,CSE_enable=False")
outstringsp += lr[i].lhs+" = "+lr[i].rhs+"\n"
outsplhs.append(sp.sympify(lr[i].lhs))
outsprhs.append(sp.sympify(lr[i].rhs))
outstring += "\n\n\n/* SEOBNR \partial_x H expression: */\n"
for i in range(len(lhss_deriv_x)):
outstring += outputC(rhss_deriv_x[i],str(lhss_deriv_x[i]),"returnstring","outCverbose=False,includebraces=False,CSE_enable=False")
outstringsp += str(lhss_deriv_x[i])+" = "+str(rhss_deriv_x[i])+"\n"
outsplhs.append(lhss_deriv_x[i])
outsprhs.append(rhss_deriv_x[i])
# for i in range(len(outsplhs)):
# for j in range(i+1,len(outsplhs)):
# outsprhs[j] = outsprhs[j].subs(outsplhs[i],outsprhs[i])
with open("/tmp/sympy_expression.py","w") as file:
file.write("""
import sympy as sp
from outputC import *
m1,m2,x,y,z,px,py,pz,s1x,s1y,s1z,s2x,s2y,s2z,eta = sp.symbols("m1 m2 x y z px py pz s1x s1y s1z s2x s2y s2z eta",real=True)
c0k2,c1k2,c0k3,c1k3,c0k4,c1k4,c2k4,c0k5,c1k5,c2k5 = sp.symbols("c0k2 c1k2 c0k3 c1k3 c0k4 c1k4 c2k4 c0k5 c1k5 c2k5",real=True)
KK,k5l,b3,bb3,d1,d1v2,dheffSS,dheffSSv2 = sp.symbols("KK k5l b3 bb3 d1 d1v2 dheffSS dheffSSv2",real=True)
tortoise,copysignresult = sp.symbols("tortoise copysignresult",real=True)
""")
for i in range(len(lr)):
file.write(lr[i].lhs+" = "+"sp.symbols(\""+lr[i].lhs+"\")\n")
file.write("\n")
for i in range(len(lhss_deriv_x)):
file.write(str(lhss_deriv_x[i]).replace("prm","prm_x")+" = "+str(rhss_deriv_x[i]).replace("sqrt(","sp.sqrt(").replace("log(","sp.log(").replace("Abs(","sp.Abs(").replace("prm","prm_x")+"\n")
for i in range(len(lhss_deriv_y)):
file.write(str(lhss_deriv_y[i]).replace("prm","prm_y")+" = "+str(rhss_deriv_y[i]).replace("sqrt(","sp.sqrt(").replace("log(","sp.log(").replace("Abs(","sp.Abs(").replace("prm","prm_y")+"\n")
for i in range(len(lhss_deriv_z)):
file.write(str(lhss_deriv_z[i]).replace("prm","prm_z")+" = "+str(rhss_deriv_z[i]).replace("sqrt(","sp.sqrt(").replace("log(","sp.log(").replace("Abs(","sp.Abs(").replace("prm","prm_z")+"\n")
for i in range(len(lhss_deriv_px)):
file.write(str(lhss_deriv_px[i]).replace("prm","prm_px")+" = "+str(rhss_deriv_px[i]).replace("sqrt(","sp.sqrt(").replace("log(","sp.log(").replace("Abs(","sp.Abs(").replace("prm","prm_px")+"\n")
for i in range(len(lhss_deriv_py)):
file.write(str(lhss_deriv_py[i]).replace("prm","prm_py")+" = "+str(rhss_deriv_py[i]).replace("sqrt(","sp.sqrt(").replace("log(","sp.log(").replace("Abs(","sp.Abs(").replace("prm","prm_py")+"\n")
for i in range(len(lhss_deriv_pz)):
file.write(str(lhss_deriv_pz[i]).replace("prm","prm_pz")+" = "+str(rhss_deriv_pz[i]).replace("sqrt(","sp.sqrt(").replace("log(","sp.log(").replace("Abs(","sp.Abs(").replace("prm","prm_pz")+"\n")
for i in range(len(lhss_deriv_s1x)):
file.write(str(lhss_deriv_s1x[i]).replace("prm","prm_s1x")+" = "+str(rhss_deriv_s1x[i]).replace("sqrt(","sp.sqrt(").replace("log(","sp.log(").replace("Abs(","sp.Abs(").replace("prm","prm_s1x")+"\n")
for i in range(len(lhss_deriv_s1y)):
file.write(str(lhss_deriv_s1y[i]).replace("prm","prm_s1y")+" = "+str(rhss_deriv_s1y[i]).replace("sqrt(","sp.sqrt(").replace("log(","sp.log(").replace("Abs(","sp.Abs(").replace("prm","prm_s1y")+"\n")
for i in range(len(lhss_deriv_s1z)):
file.write(str(lhss_deriv_s1z[i]).replace("prm","prm_s1z")+" = "+str(rhss_deriv_s1z[i]).replace("sqrt(","sp.sqrt(").replace("log(","sp.log(").replace("Abs(","sp.Abs(").replace("prm","prm_s1z")+"\n")
for i in range(len(lhss_deriv_s2x)):
file.write(str(lhss_deriv_s2x[i]).replace("prm","prm_s2x")+" = "+str(rhss_deriv_s2x[i]).replace("sqrt(","sp.sqrt(").replace("log(","sp.log(").replace("Abs(","sp.Abs(").replace("prm","prm_s2x")+"\n")
for i in range(len(lhss_deriv_s2y)):
file.write(str(lhss_deriv_s2y[i]).replace("prm","prm_s2y")+" = "+str(rhss_deriv_s2y[i]).replace("sqrt(","sp.sqrt(").replace("log(","sp.log(").replace("Abs(","sp.Abs(").replace("prm","prm_s2y")+"\n")
for i in range(len(lhss_deriv_s2z)):
file.write(str(lhss_deriv_s2z[i]).replace("prm","prm_s2z")+" = "+str(rhss_deriv_s2z[i]).replace("sqrt(","sp.sqrt(").replace("log(","sp.log(").replace("Abs(","sp.Abs(").replace("prm","prm_s2z")+"\n")
file.write("""
sp.cse(Hrealprm_x)
# outputC([Hrealprm_x,Hrealprm_y,Hrealprm_z,Hrealprm_px,Hrealprm_py,Hrealprm_pz,
# Hrealprm_s1x,Hrealprm_s1y,Hrealprm_s1z,Hrealprm_s2x,Hrealprm_s2y,Hrealprm_s2z],
# ["Hrealprm_x","Hrealprm_y","Hrealprm_z","Hrealprm_px","Hrealprm_py","Hrealprm_pz",
# "Hrealprm_s1x","Hrealprm_s1y","Hrealprm_s1z","Hrealprm_s2x","Hrealprm_s2y","Hrealprm_s2z"],
# "/tmp/outC.h","outCverbose=False,includebraces=False,enable_SIMD=True")
""")
# !cp /tmp/sympy_expression.py .
# !pypy sympy_expression.py
!cp /tmp/sympy_expression.py .
!pypy sympy_expression.py
# with open("/tmp/outC.h","r") as file:
# file.read()