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 the following, the ODE system directly integrates the quaternion rotors describing the spin directions and orbital angular velocity vector. This allows us to reduce the number of variables in the system and automatically satisfy certain constraints, which should lead to better numerical behavior.
Note that in this formulation, Phi is a completely unnecessary variable. It doesn't need to be computed because it is implict in the definition of the frame. It is calculated here simply because it is easy to do, and might be handy to have around. Of course, it could be derived from knowledge of the frame alone.
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 = True
from sympy import N
from sympy import Rational as frac # Rename for similarity to latex
In [2]:
with open('PNApproximants_Q.ipp', 'w') as f :
f.write("""// File produced automatically by OrbitalEvolutionCodeGen_Q.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_Q : public TaylorTn_Q {{
private:
{Declarations}
double Phi;
const bool EvolveSpin1, EvolveSpin2;
public:
TaylorTn_{PNOrbitalEvolutionOrderString}PN_Q({InputArguments}) :
{Initializations},
Phi(0.0), EvolveSpin1(S_chi1.normsquared()>1e-12), EvolveSpin2(S_chi2.normsquared()>1e-12)
{{ }}
void Recalculate(double t, const double* y) {{
v = y[0];
rfrak_chi1_x = y[1];
rfrak_chi1_y = y[2];
rfrak_chi2_x = y[3];
rfrak_chi2_y = y[4];
rfrak_frame_x = y[5];
rfrak_frame_y = y[6];
rfrak_frame_z = y[7];
Phi = y[8];
{Evaluations}
}}
{MemberFunctions}
""".format(PNOrbitalEvolutionOrderString=PNOrbitalEvolutionOrderString,
InputArguments=CodeConstructor.CppInputArguments(11),
Declarations=CodeConstructor.CppDeclarations(2),
Initializations=CodeConstructor.CppInitializations(4),
Evaluations=CodeConstructor.CppEvaluations(4),
MemberFunctions=CodeOutput.CodeConstructor(PNVariables, PrecessionVelocities).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) {{
std::vector<double> rfrak_frame(3);
rfrak_frame[0] = y[5];
rfrak_frame[1] = y[6];
rfrak_frame[2] = y[7];
const std::vector<double> rfrakdot_frame = FrameFromAngularVelocity_Integrand(rfrak_frame, OmegaVec().vec());
dydt[0] = dvdt;
if(EvolveSpin1) {{
FrameFromAngularVelocity_2D_Integrand(y[1], y[2],
(S_chi1.inverse()*OmegaVec_chiVec_1()*S_chi1).vec(),
dydt[1], dydt[2]);
}} else {{
dydt[1] = 0.0;
dydt[2] = 0.0;
}}
if(EvolveSpin2) {{
FrameFromAngularVelocity_2D_Integrand(y[3], y[4],
(S_chi2.inverse()*OmegaVec_chiVec_2()*S_chi2).vec(),
dydt[3], dydt[4]);
}} else {{
dydt[3] = 0.0;
dydt[4] = 0.0;
}}
dydt[5] = rfrakdot_frame[0];
dydt[6] = rfrakdot_frame[1];
dydt[7] = rfrakdot_frame[2];
dydt[8] = v*v*v/M;
return GSL_SUCCESS; // GSL expects this if everything went well
}}
}}; // class TaylorTn_{PNOrbitalEvolutionOrderString}PN_Q : public TaylorTn_Q
""".format(PNOrbitalEvolutionOrderString=PNOrbitalEvolutionOrderString))
print("All done")