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')
from sympy import N
from sympy import Rational as frac # Rename for similarity to latex
execnotebook('../PNTerms/Variables_Q.ipynb')
from Utilities import CodeOutput
execnotebook('../PNTerms/WaveformModes.ipynb')
In [2]:
WaveformModeTerms = [WaveformModes_NoSpin, WaveformModes_Spin_Symm, WaveformModes_Spin_Asymm]
for Term in WaveformModeTerms:
PNVariables.update(Term)
We use the co-orbital frame to compute waveforms, rather than the $(x,y,z)$ frame used to computer the orbital dynamics. This means that we need to change a few of the definitions of the variables that go into our calculations. The following commands overwrite the variables defined in Variables_Q.ipynb.
In [3]:
# For some reason I have to both overwrite these variables and pop them...
PNVariables.AddVariable('nHat', constant=True, fundamental=True,
substitution_atoms=[], substitution='Quaternions::xHat', datatype='Quaternions::Quaternion');
PNVariables.AddVariable('lambdaHat', constant=True, fundamental=True,
substitution_atoms=[], substitution='Quaternions::yHat', datatype='Quaternions::Quaternion');
PNVariables.AddVariable('ellHat', constant=True, fundamental=True,
substitution_atoms=[], substitution='Quaternions::zHat', datatype='Quaternions::Quaternion');
PNVariables.pop(nHat);
PNVariables.pop(lambdaHat);
PNVariables.pop(ellHat);
PNVariables.AddBasicVariables('chiVec1,chiVec2', datatype='Quaternions::Quaternion')
PNVariables.AddDerivedVariable('chi1_n', substitution_atoms=[chiVec1,nHat], substitution='chiVec1[1]');
PNVariables.AddDerivedVariable('chi1_lambda', substitution_atoms=[chiVec1,lambdaHat], substitution='chiVec1[2]');
PNVariables.AddDerivedVariable('chi1_ell', substitution_atoms=[chiVec1,ellHat], substitution='chiVec1[3]');
PNVariables.AddDerivedVariable('chi2_n', substitution_atoms=[chiVec2,nHat], substitution='chiVec2[1]');
PNVariables.AddDerivedVariable('chi2_lambda', substitution_atoms=[chiVec2,lambdaHat], substitution='chiVec2[2]');
PNVariables.AddDerivedVariable('chi2_ell', substitution_atoms=[chiVec2,ellHat], substitution='chiVec2[3]');
In [4]:
%%time
with open('PNWaveformModes.ipp', 'w') as f :
f.write("""// File produced automatically by WaveformModeCodeGen.ipynb
class WaveformModes_Base {{
public:
virtual std::vector<std::complex<double> > operator()(
const double v_k,
const double chi1_n_k, const double chi1_lambda_k, const double chi1_ell_k,
const double chi2_n_k, const double chi2_lambda_k, const double chi2_ell_k) = 0;
virtual std::vector<std::complex<double> > operator()(
const double v_k, const std::vector<double>& chi1, const std::vector<double>& chi2)
{{
return this->operator()(v_k, chi1[0], chi1[1], chi1[2], chi2[0], chi2[1], chi2[2]);
}}
}};
const unsigned int ellMax = {ellMax};
const std::complex<double> I(0,1.0);
inline std::complex<double> conjugate(const std::complex<double>& a) {{ return std::conj(a); }}""".format(ellMax=ellMax))
for PNWaveformModeOrder in [frac(n,2) for n in range(0,8)]:
print("Working on {0} PN...".format(PNWaveformModeOrder))
PNWaveformModeOrderString = str(N(PNWaveformModeOrder,2)).replace('.','p')
ModeExpressions = PNCollection()
from textwrap import TextWrapper
wrapper = TextWrapper(width=120)
wrapper.initial_indent = ' '*4
wrapper.subsequent_indent = wrapper.initial_indent+' '
LM = [[ell,m] for ell in range(2,ellMax+1) for m in range(-ell,ell+1)]
Evaluations = [wrapper.fill('std::vector<std::complex<double> > Modes({0});'.format(len(LM))),
wrapper.fill('std::complex<double> Symm, Asymm;')]
for ell in range(2, ellMax+1):
for m in range(0, ell+1):
Symm = (SymmetricWaveformModes([ell,m],PNOrder=PNWaveformModeOrder)).subs(log(v), logv)
Asymm = (AsymmetricWaveformModes([ell,m],PNOrder=PNWaveformModeOrder)).subs(log(v), logv)
code1 = "Symm = {0};".format(sympy.ccode(Symm))
code2 = "Asymm = {0};".format(sympy.ccode(Asymm))
code3 = "Modes[{2}] = Symm + Asymm;".format(ell, m, LM.index([ell,m]))
code4 = "Modes[{2}] = {3}std::conj(Symm - Asymm);".format(ell, -m, LM.index([ell,-m]),
('' if ((ell%2)==0) else '-'))
ModeExpressions.AddDerivedVariable('rhOverM_{0}_{1}_Symm'.format(ell,m),
Symm, datatype='std::complex<double>')
ModeExpressions.AddDerivedVariable('rhOverM_{0}_{1}_Asymm'.format(ell,m),
Asymm, datatype='std::complex<double>')
Evaluations.append(wrapper.fill("// (ell, m) = ({0}, +/- {1})".format(ell, m)))
Evaluations.append(wrapper.fill(code1))
Evaluations.append(wrapper.fill(code2))
Evaluations.append(wrapper.fill(code3))
if m!=0:
Evaluations.append(wrapper.fill(code4))
CodeConstructor = CodeOutput.CodeConstructor(PNVariables, ModeExpressions)
if(PNWaveformModeOrder>0.9):
Updates = """ chiVec1 = Quaternions::Quaternion(0., chi1_n_k, chi1_lambda_k, chi1_ell_k);
chiVec2 = Quaternions::Quaternion(0., chi2_n_k, chi2_lambda_k, chi2_ell_k);
"""
else:
Updates = ""
# Start the class, write the declarations, initializations, etc.
f.write("""
class WaveformModes_{PNWaveformModeOrderString}PN : public WaveformModes_Base {{
private:
{Declarations}
public:
WaveformModes_{PNWaveformModeOrderString}PN({InputArguments}) :
{Initializations}
{{ }}
using WaveformModes_Base::operator();
std::vector<std::complex<double> > operator()(
const double v_k,
const double chi1_n_k, const double chi1_lambda_k, const double chi1_ell_k,
const double chi2_n_k, const double chi2_lambda_k, const double chi2_ell_k)
{{
v = v_k;
{Updates}
{Evaluations}
{Computations}
return Modes;
}}
}}; // class WaveformModes_{PNWaveformModeOrderString}PN : public WaveformModes_Base
""".format(PNWaveformModeOrderString=PNWaveformModeOrderString,
InputArguments=CodeConstructor.CppInputArguments(22),
Declarations=CodeConstructor.CppDeclarations(2),
Initializations=CodeConstructor.CppInitializations(4),
Updates=Updates,
Evaluations=CodeConstructor.CppEvaluations(4),
MemberFunctions=CodeConstructor.CppExpressionsAsFunctions(2),
Computations='\n'.join(Evaluations) ))
print("All done")