Finally, we write the code file itself. Most of the text below will appear literally in the file. Statements in single curly braces will be substituted by the format statement at the end of the string; double curly braces will just appear as literal braces in the output.
In [1]:
# Always run this first
# NOTE: Do not define new basic variables in this notebook;
# define them in Variables_Q.ipynb. Use this notebook
# to define new expressions built from those variables.
from __future__ import division # This needs to be here, even though it's in Variables.py
import sys
sys.path.insert(0, '..') # Look for modules in directory above this one
execfile('../Utilities/ExecNotebook.ipy')
UseQuaternions = False
from sympy import N
from sympy import Rational as frac # Rename for similarity to latex
In [2]:
with open('PNApproximants.ipp', 'w') as f :
f.write("""// File produced automatically by OrbitalEvolutionCodeGen.ipynb""")
for PNOrbitalEvolutionOrder in [frac(n,2) for n in range(0,13)]:
print("Working on {0} PN...".format(PNOrbitalEvolutionOrder))
PNOrbitalEvolutionOrderString = str(N(PNOrbitalEvolutionOrder,2)).replace('.','p')
%cd -q ../PNTerms
execnotebook('OrbitalEvolution.ipynb')
ExpressionsForMemberFunctions = CodeOutput.CodeConstructor(PNVariables, PrecessionVelocities)
execnotebook('AngularMomentum.ipynb')
%cd -q -
PrecessionVelocities.AddDerivedVariable('OrbitalAngularMomentum',
AngularMomentumExpression(PNOrder=PNOrbitalEvolutionOrder),
datatype=ellHat.datatype)
CodeConstructor.AddDependencies(PrecessionVelocities)
# Start the class, write the declarations, initializations, etc.
f.write("""
class TaylorTn_{PNOrbitalEvolutionOrderString}PN : public TaylorTn {{
private:
{Declarations}
double Phi, gamma;
public:
TaylorTn_{PNOrbitalEvolutionOrderString}PN({InputArguments}) :
{Initializations}, Phi(0.0), gamma(0.0)
{{
ellHat[0] = ellHat_x;
ellHat[1] = ellHat_y;
ellHat[2] = ellHat_z;
nHat[0] = nHat_x;
nHat[1] = nHat_y;
nHat[2] = nHat_z;
}}
void Recalculate(double t, const double* y) {{
v = y[0];
chi1_x = y[1];
chi1_y = y[2];
chi1_z = y[3];
chi2_x = y[4];
chi2_y = y[5];
chi2_z = y[6];
ellHat_x = y[7];
ellHat_y = y[8];
ellHat_z = y[9];
Phi = y[10];
gamma = y[11];
const Quaternion Rax =
sqrtOfRotor(-normalized(Quaternion(0., ellHat_x, ellHat_y, ellHat_z))*zHat);
const Quaternion R = Rax * exp(((gamma+Phi)/2.)*zHat);
const Quaternion nHatQ = R*xHat*R.conjugate();
nHat_x = nHatQ[1];
nHat_y = nHatQ[2];
nHat_z = nHatQ[3];
{Evaluations}
}}
{MemberFunctions}
""".format(PNOrbitalEvolutionOrderString=PNOrbitalEvolutionOrderString,
InputArguments=CodeConstructor.CppInputArguments(17),
Declarations=CodeConstructor.CppDeclarations(2),
Initializations=CodeConstructor.CppInitializations(4),
Evaluations=CodeConstructor.CppEvaluations(4),
MemberFunctions=ExpressionsForMemberFunctions.CppExpressionsAsFunctions(2)))
# Write the external interfaces to the ODE stepper
for n,Expressions in zip([1,4,5], [T1Expressions,T4Expressions,T5Expressions]) :
f.write("""
int TaylorT{n}(double t, const double* y, double* dydt) {{
Recalculate(t, y);
if(v>=1.0) {{ return GSL_EDOM; }} // Beyond domain of PN validity
{Computations}
if(dvdt_T{n}<1.0e-12) {{ return GSL_EDIVERGE; }} // v is decreasing
return CommonRHS(dvdt_T{n}, y, dydt);
}}
""".format(n=n,
Computations=CodeOutput.CodeConstructor(PNVariables, Expressions).CppEvaluateExpressions() ))
# Finish up, writing the common RHS function and closing the class
f.write("""
int CommonRHS(const double dvdt, const double* y, double* dydt) {{
const std::vector<double> Omega1 = OmegaVec_chiVec_1();
const std::vector<double> Omega2 = OmegaVec_chiVec_2();
const std::vector<double> Omega = OmegaVec_ellHat();
dydt[0] = dvdt;
cross(&Omega1[0], &y[1], &dydt[1]);
cross(&Omega2[0], &y[4], &dydt[4]);
cross(&Omega[0], &y[7], &dydt[7]);
dydt[10] = v*v*v/m;
const Quaternion adot(0., dydt[7], dydt[8], dydt[9]);
const Quaternion Rax =
sqrtOfRotor(-Quaternion(0., ellHat_x, ellHat_y, ellHat_z).normalized()*zHat);
const Quaternion Raxdot = ( (-1.0/std::sqrt(2+2*y[9]))*adot*zHat - (dydt[9]/(2+2*y[9]))*Rax );
dydt[11] = -2*(Rax.conjugate() * Raxdot)[3];
return GSL_SUCCESS; // GSL expects this if everything went well
}}
}}; // class TaylorTn_{PNOrbitalEvolutionOrderString}PN : public TaylorTn
""".format(PNOrbitalEvolutionOrderString=PNOrbitalEvolutionOrderString))
print("All done")