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
m1PlusetaKKsq = m1PlusetaKK*m1PlusetaKK
invm1PlusetaKK = 1/m1PlusetaKK
k1sq = k1*k1
k1cu = k1sq*k1
k1ft = k1cu*k1
k2 = 0.5*k1*(k1 - 4.*m1PlusetaKK) - a2*m1PlusetaKKsq*k0
k2sq = k2*k2
k3 = -k1cu/3. + k1*k2 + k1sq*m1PlusetaKK - 2.*(k2 - m1PlusetaKK)*m1PlusetaKK - a2*k1*m1PlusetaKKsq
k4 = 1./12.*(6*a2*(k1sq - 2*k2)*m1PlusetaKKsq + 3*k1ft - 8*m1PlusetaKK*k1cu - 12*k2*k1sq + 12*(2*m1PlusetaKK*k2 + k3)*k1 + 12*(94./3. - 41./32.*sp.pi*sp.pi)*m1PlusetaKKsq + 6*(k2*k2 - 4*k3*m1PlusetaKK))
k5 = m1PlusetaKKsq*(-4237./60. + 128./5.*EMgamma + 2275./512.*sp.pi*sp.pi - a2*(k1cu - 3.*k1*k2 + 3.*k3)/3. - (k1ft*k1 - 5.*k1cu*k2 + 5.*k1*k2sq + 5.*k1sq*k3 - 5.*k2*k3 - 5.*k1*k4)/(5.*m1PlusetaKKsq) + (k1ft - 4.*k1sq*k2 + 2.*k2sq + 4.*k1*k3 - 4.*k4)/(2*m1PlusetaKK) + (256./5.)*sp.log(2))
k5l = (64./5.)*m1PlusetaKKsq
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(sp.Abs(onepluslogarg))
deltaU = sp.Abs(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 = sp.Abs(w2*w2 - a2*deltaT*xi2)
rho2xi2Lambda = rho2*xi2*Lambda
invrho2xi2Lambda = 1/rho2xi2Lambda
invrho2 = 1/rho2
invxi2 = 1/xi2
invLambda = 1/Lambda
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)
b3 = 0
bb3 = 0
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(sp.Abs(deltaT*deltaR))/w2
csi1  =  1 + (1-sp.Abs(1-tortoise)) * (csi - 1)
csi2  =  1 + (sp.Rational(1,2)-0.5*sp.sign(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
deltaSigmaStarUSCOREx = deltaSigmaStarUSCOREx2 + d1v2*sigmaKerr0*etau3
deltaSigmaStarUSCOREy = deltaSigmaStarUSCOREy2 + d1v2*sigmaKerr1*etau3
deltaSigmaStarUSCOREz = deltaSigmaStarUSCOREz2 + 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 + dheffSSv2*Hpt1
Hreal = sp.sqrt(1 + 2*eta*(sp.Abs(H) - 1))