Macros
#define	Power(A, B) pow(A,B)

#define	Sqrt(A) sqrt(A)

#define	Log(A) log(A)

Variables
sigmaKerr	data [0] = s1Vec->data[0] + s2Vec->data[0]

const double	invmass1mass2 = 1./(mass1*mass2)

const double	mass2overmass1 = invmass1mass2((mass2)(mass2))

const double	mass1overmass2 = invmass1mass2((mass1)(mass1))

const double	s1dots1 = ((s1Vec->data[0])(s1Vec->data[0])) + ((s1Vec->data[1])(s1Vec->data[1])) + ((s1Vec->data[2])*(s1Vec->data[2]))

const double	s2dots2 = ((s2Vec->data[0])(s2Vec->data[0])) + ((s2Vec->data[1])(s2Vec->data[1])) + ((s2Vec->data[2])*(s2Vec->data[2]))

const double	r2 = ((x->data[0])(x->data[0])) + ((x->data[1])(x->data[1])) + ((x->data[2])*(x->data[2]))

const double	r = Sqrt(r2)

const double	u = 1./r

const double	u2 = ((u)*(u))

const double	u3 = u*u2

const double	u4 = ((u2)*(u2))

const double	u5 = u*u4

const double	etau3 = eta*u3

const double	etau4 = eta*u4

const double	nx = u*x->data[0]

const double	ny = u*x->data[1]

const double	nz = u*x->data[2]

const double	sKerrUSCOREx = sigmaKerr->data[0]

const double	sKerrUSCOREy = sigmaKerr->data[1]

const double	sKerrUSCOREz = sigmaKerr->data[2]

const double	sStarUSCOREx = sigmaStar->data[0]

const double	sStarUSCOREy = sigmaStar->data[1]

const double	sStarUSCOREz = sigmaStar->data[2]

const double	a2 = ((sKerrUSCOREx)(sKerrUSCOREx)) + ((sKerrUSCOREy)(sKerrUSCOREy)) + ((sKerrUSCOREz)*(sKerrUSCOREz))

const double	a4 = ((a2)*(a2))

const double	a = Sqrt(a2)

const double	m1PlusetaKK = -1. + coeffs->KK*eta

const double	invm1PlusetaKK = 1./m1PlusetaKK

coeffs	k0 = coeffs->KK*(-1. + m1PlusetaKK)

coeffs	k1 = -2.(coeffs->k0 + coeffs->KK)m1PlusetaKK

coeffs	k2 = c0k2 + a2*c1k2

coeffs	k3 = c0k3 + a2*c1k3

coeffs	k4 = c0k4 + a2c1k4 + a4c2k4

coeffs	k5 = c0k5 + a2c1k5 + a4c2k5

const double	e3USCOREx = e3_x

const double	e3USCOREy = e3_y

const double	e3USCOREz = e3_z

const double	costheta = e3USCORExnx + e3USCOREyny + e3USCOREz*nz

const double	xi2 = 1. - ((costheta)*(costheta))

const double	xiUSCOREx = -(e3USCOREzny) + e3USCOREynz

const double	xiUSCOREy = e3USCOREznx - e3USCORExnz

const double	xiUSCOREz = -(e3USCOREynx) + e3USCORExny

const double	vx = -(nzxiUSCOREy) + nyxiUSCOREz

const double	vy = nzxiUSCOREx - nxxiUSCOREz

const double	vz = -(nyxiUSCOREx) + nxxiUSCOREy

const double	w2 = a2 + r2

const double	rho2 = a2((costheta)(costheta)) + r2

const double	bulk = ((invm1PlusetaKK)(invm1PlusetaKK)) + 2.invm1PlusetaKKu + a2u2

const double	logu = Log(u)

const double	logarg = coeffs->k1u + coeffs->k2u2 + coeffs->k3u3 + coeffs->k4u4 + coeffs->k5u5 + coeffs->k5llogu*u5

const double	onepluslogarg = 1. + logarg

const double	invonepluslogarg = 1./onepluslogarg

const double	logTerms = 1. + coeffs->k0eta + etaLog(onepluslogarg)

const double	deltaU = bulk*logTerms

const double	deltaT = deltaU*r2

const double	deltaUUSCOREupt7 = coeffs->k5 + coeffs->k5l*logu

const double	deltaUUSCOREupt6 = 4.coeffs->k4 + 5.deltaUUSCOREupt7*u

const double	deltaUUSCOREupt5 = 3.coeffs->k3 + deltaUUSCOREupt6u

const double	deltaUUSCOREupt4 = 2.coeffs->k2 + deltaUUSCOREupt5u

const double	deltaUUSCOREupt3 = coeffs->k1 + deltaUUSCOREupt4*u

const double	deltaUUSCOREupt2 = invm1PlusetaKK + a2*u

const double	deltaUUSCOREupt1 = bulkdeltaUUSCOREupt3eta

const double	deltaUUSCOREu = deltaUUSCOREupt1invonepluslogarg + 2.deltaUUSCOREupt2*logTerms

const double	deltaTUSCOREr = -deltaUUSCOREu + 2.deltaUr

const double	Lambda = ((w2)(w2)) - a2deltaT*xi2

const double	rho2xi2Lambda = Lambdarho2xi2

const double	invrho2xi2Lambda = 1./rho2xi2Lambda

const double	invrho2 = invrho2xi2LambdaLambdaxi2

const double	invxi2 = invrho2xi2LambdaLambdarho2

const double	invLambda = invrho2xi2Lambdarho2xi2

const double	invLambdasq = ((invLambda)*(invLambda))

const double	rho2invLambda = invLambda*rho2

const double	expnu = Sqrt(deltaT*rho2invLambda)

const double	expMU = Sqrt(rho2)

const double	expMUexpnu = expMU*expnu

const double	expMUsq = ((expMU)*(expMU))

const double	expnusq = ((expnu)*(expnu))

const double	expMUsqexpnusq = expMUsq*expnusq

const double	invexpnuexpMU = 1./expMUexpnu

const double	invexpMU = expnu*invexpnuexpMU

const double	invexpMUsq = ((invexpMU)*(invexpMU))

const double	expnuinvexpMU2 = expnu*invexpMUsq

const double	invexpMUcubinvexpnu = invexpMUsq*invexpnuexpMU

const double	DD = 1. + Log(1 + (52. - 6.eta)etau3 + 6.etau2)

const double	deltaR = DD*deltaT

const double	qq = -6.eta(-1.3333333333333333 + 1.*eta)

const double	ww = a(2.r + (a2coeffs->b3 + coeffs->bb3)eta*u)

const double	B = Sqrt(deltaT)

const double	sqrtdeltaT = B

const double	sqrtdeltaR = Sqrt(deltaR)

const double	deltaTsqrtdeltaR = deltaT*sqrtdeltaR

const double	sqrtdeltaTdeltaTsqrtdeltaR = deltaTsqrtdeltaR*sqrtdeltaT

const double	invdeltaTsqrtdeltaTsqrtdeltaR = 1./sqrtdeltaTdeltaTsqrtdeltaR

const double	invdeltaT = invdeltaTsqrtdeltaTsqrtdeltaRsqrtdeltaRsqrtdeltaT

const double	invsqrtdeltaT = deltaTsqrtdeltaR*invdeltaTsqrtdeltaTsqrtdeltaR

const double	invsqrtdeltaR = deltaTinvdeltaTsqrtdeltaTsqrtdeltaRsqrtdeltaT

const double	w = invLambda*ww

const double	LambdaUSCOREr = 4.rw2 - a2deltaTUSCORErxi2

const double	wwUSCOREr = -1.a(-2. + a2coeffs->b3etau2 + coeffs->bb3eta*u2)

const double	BR = (0.5deltaTUSCOREr - deltaTinvsqrtdeltaR)*invsqrtdeltaT

const double	wr = invLambdasq(-(LambdaUSCORErww) + Lambda*wwUSCOREr)

const double	nurpt2 = w2(-4.deltaTr + deltaTUSCORErw2)

const double	nurpt1 = invdeltaT*nurpt2

const double	nur = 0.5invLambdanurpt1 + invrho2*r

const double	mur = -invsqrtdeltaR + invrho2*r

const double	a2costheta = a2*costheta

const double	wcospt2 = deltaT*ww

const double	wcospt1 = invLambdasq*wcospt2

const double	wcos = -2.a2costhetawcospt1

const double	nucospt3 = invLambda*invrho2

const double	nucospt2 = nucospt3*w2

const double	nucospt1 = a2costheta*nucospt2

const double	nucos = nucospt1*(-deltaT + w2)

const double	mucos = a2costheta*invrho2

const double	etaover12r = 0.08333333333333333etau

const double	csi = Sqrt(deltaR*deltaT)/w2