Three-Body Expansion
 All Files Functions Variables Macros
Functions | Variables
threebody_plane.c File Reference
#include <stdio.h>
#include <math.h>
#include <time.h>
#include "chrlocal.h"
#include "chronos_msg.h"
#include "chronos_cmn.h"
#include "chronos_def.h"
#include "chronos_ext.h"
#include "n_corpi.h"

Go to the source code of this file.

Functions

void a2L (DMHDR *srcdmh, int eccdeg, DMHDR *dstdmh[DEG_M][DEG_E])
 This routine replaces the terms $\frac{a_1^x a_2^y}{{a_2^*}^{x+y}}$ in the function srcdmh with their expansions in $L$ and stores the result in the function dstdmh[DEG_M][eccdeg].
void add_kinetic (DMHDR *kindmh[DEG_M][DEG_E], DMHDR *hamdmh[DEG_M][DEG_E])
 This routine adds the kinetic energy, $T$, storesd in the function kindmh[DEG_M][DEG_E], to the Hamiltonian, $H$, storesd in the function hamdmh[DEG_M][DEG_E].
void compute_kinetic (DMHDR *kindmh[DEG_M][DEG_E], DMHDR *wdmh[DEG_E])
 This routine computes the expansion of the kinetic energy, $T$, in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.
void compute_utilities (FILE *ifp, FILE *ofp1, FILE *ofp2)
 This routine computes the initial conditions in Poincare' canonical heliocentric variables and some useful data.
void cuteps (DMHDR *fndmh, DOUBLE eps)
 This routine set to zero all the coefficients of the function fndmh with absolute value smaller than the threshold eps.
void cutepsrel (DMHDR *fndmh, DOUBLE eps)
 This routine set to zero all the coefficients of the function fndmh with absolute value smaller than the threshold eps times the absolute value of the maximum coefficient of the function.
void e2X (DMHDR *fdmh[DEG_E], DMHDR *e2Xdmh[DEG_E][DEG_E])
 This routine replaces the terms $e^i$ in the function fdmh[DEG_E] with their expansions in $|X|$ (storesd in the function e2Xdmh[i][DEG_E]) and stores again in the function fdmh[DEG_E].
DOUBLE fhyp (double a, double b, double c, DOUBLE x, DOUBLE eps)
 This routine computes the hypergeometric function $F(a,b,c;x)$.
INT2 flagsincos (INT2 *k, int n)
 This routine checks if the array k[n] represents a sinus or a cosinus.
void init_a2L (INT2 ka[2], DMHDR *fndmh[DEG_M])
 This routine computes the expansion of $\frac{a_1^{{\tt ka[0]}} a_2^{{\tt ka[1]}}}{{a_2^*}^{{\tt ka[0]+ka[1]}}}$ in $L$.
void init_alpha (DMHDR *alpha2dmh[DEG_M], DMHDR *alphajdmh[(3 *(DEG_E-1)+MM_CUT)/2+1][DEG_M])
 This routine computes the expansion of $\alpha^2$ and $\alpha^j$ in $L$.
void init_bij (int i, int j, DMHDR *alpha2dmh[DEG_M], DMHDR *alphajdmh[DEG_M], DMHDR *lapldmh[DEG_M])
 This routine computes the expansion of the Laplace coefficient $b_{i+\frac{1}{2}}^{(j)}(\alpha)$ in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$. If $j=0$, computes $\frac{1}{2}b_{i+\frac{1}{2}}^{(0)}(\alpha)$.
void init_e2X (DMHDR *e2Xdmh[DEG_E][DEG_E])
 This routine computes the expansions of the powers of the eccentricity in $|X|$ and $(\lambda, h)$.
void init_fhyp (double a, double b, double c, DMHDR *xydmh[DEG_M], DOUBLE reps, DMHDR *fhypdmh[DEG_M])
 This routine computes the expansion of the hypergeometric function $F(a,b,c;x+y)$ in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.
void init_h0 (DMHDR *hamdmh[DEG_M+1])
 This routine computes the expansion of the unperturbed Hamiltonian, $H_0$, in $L$, $\lambda$ and $(\xi, \eta)$.
void init_h0_split (DMHDR *ham1dmh[DEG_M+1], DMHDR *ham2dmh[DEG_M+1])
 This routine computes the expansion of the unperturbed Hamiltonian, $H_0$, in $L$, $\lambda$ and $(\xi, \eta)$, splitting the parts related to the inner and outer planet.
void init_pol2poincsec (INT2 k[2], DMHDR *fndmh[DEG_M])
 This routine computes the expansion of $|X_1|^{{\tt k[0]}} |X_2|^{{\tt k[1]}}$ in $L$, $\lambda$ and $(\xi, \eta)$. The expansion is explicit only in $L$.
void init_roveraf (DMHDR *roverafdmh[DEG_E])
 This routine computes the expansion of $\frac{\|r\|}{a}$ in $e$ and $f$. The expansion is explicit only in $f$.
void init_W (DMHDR *wdmh[DEG_E])
 This routine computes the expansion of $W(f,e)$ in $e$ and $f$. The expansion is explicit only in $f$.
void inner2outer (DMHDR *f1dmh, DMHDR *f2dmh)
 This routine shifts the indices corresponding to the first planet into the ones of the second planet.
void kin_ang (DMHDR *kinangdmh[DEG_M][DEG_E], DMHDR *wdmh[DEG_E])
 This routine computes the expansion of $T_{ang}$ in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.
void kin_ecc (DMHDR *kineccdmh[DEG_E])
 This routine computes the expansion of $T_{ecc}$ in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.
void kin_fastact (DMHDR *kinazveldmh[DEG_M])
 This routine computes the expansion of $T_{L}$ in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.
void laplace (DMHDR *svilazdmh[DEG_M][DEG_E][DEG_E], DMHDR *lapldmh[(3 *(DEG_E-1)+MM_CUT)/2+1][DEG_M], DMHDR *fndmh[DEG_M][DEG_E], int i, int ecc)
 This routine computes the expansion of the term related to the index i and of degree ecc in the eccentricity of the expansion of the potential energy, in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$..
void one2two_frag (DMHDR *srcdmh, DMHDR *dstdmh)
 This routine transforms an expansion in one trigonometric variable, into an expansion in three variables, using the second one.
void prod_D2 (DMHDR *f1dmh[DEG_E], DMHDR *f2dmh[DEG_E], DMHDR *resdmh[DEG_E], DOUBLE coef, UNS2 flag_D2)
 This routine computes the product of the functions f1dmh[DEG_E] and f2dmh[DEG_E] times the coefficient coef.
void prod_D2_cut (DMHDR *fn1dmh, DMHDR *fn2dmh, DOUBLE coef, DMHDR *resdmh)
 This routine computes the truncated product of the functions fn1dmh and fn2dmh times the coefficient coef. The limit in the truncation depends on the degree in the parameter $D_2^*$.
void rovera (DMHDR *roveradmh[DEG_E], DMHDR *wdmh[DEG_E])
 This routine computes the expansion of $\frac{\|r\|}{a}$ in $e$ and $(\lambda, h)$.
void secular2poinc (DMHDR *srcdmh[DEG_M][DEG_E], int ecc, DMHDR *dstdmh[DEG_M][DEG_E])
 This routine replaces the terms $ |X_1|^a |X_2|^b \ {\sin \atop \cos}(\alpha h_1 + \beta h_2)$ in srcdmh[DEG_M][ecc] with their expansion in $(\xi, \eta)$ and stores the result in the function dstdmh[DEG_M][ecc].
void sincosfh (DMHDR *sinfhdmh[DEG_E], DMHDR *cosfhdmh[DEG_E], DMHDR *wdmh[DEG_E])
 This routine computes the expansions of $\sin(f-h)$ and $\cos(f-h)$ in $e$ and $(\lambda, h)$.
void tri_lie (DMHDR *triliedmh[DEG_E][DEG_E], DMHDR *wdmh[DEG_E])
 This routine transforms the expansion of $F(f,e)$ in $e$ and $f$, into $G(l,e)$ in $e$ and $l$, by using the generator $W(f,e)$ of the Lie serie. The expansion is explicit only in $f$ or $l$.
void trig2poinc (INT2 k5[9], INT2 kpol[2], INT2 ktrg[2], DOUBLE coef, DMHDR *fndmh)
 This routine transforms an expansion in $\sqrt{\Lambda}|X|$ and $h$, into an expansion in $(\xi, \eta)$.
void two2five (DMHDR *srcdmh, DMHDR *dstdmh)
 This routine transforms the expansion of the function srcdmh in $e$ and $(\lambda, h)$, into an expansion in $D_2^*$, $L$, $|X|$, $\lambda$, $h$, stored in dstdmh.
void two2five_a (DMHDR *srcdmh, DMHDR *dstdmh)
 This routine transforms the expansion of the function srcdmh in $e$ and $(\lambda, h)$, into an expansion in $D_2^*$, $a$, $|X|$, $\lambda$, $h$, stored in dstdmh. Moreover the routine performs the multiplication by $a$.
void xi_over_a2 (DMHDR *xiovera2dmh[DEG_M][DEG_E], DMHDR *wdmh[DEG_E])
 This routine computes the expansion of $\frac{\Xi}{{a_2^*}^2}$ in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.
int main ()
 This program computes the expansion of the Hamiltonian of the three-body problem in $L$, $\lambda$ and $(\xi, \eta)$.

Variables

DOUBLE M0
DOUBLE M1
DOUBLE M2
DOUBLE Mu1
DOUBLE Mu2
DOUBLE Beta1
DOUBLE Beta2
DOUBLE A_fix [2]
DOUBLE Alpha_fix
DOUBLE Lambda_fix [2]
DOUBLE D2_fix

Function Documentation

void a2L ( DMHDR *  srcdmh,
int  eccdeg,
DMHDR *  dstdmh[DEG_M][DEG_E] 
)

This routine replaces the terms $\frac{a_1^x a_2^y}{{a_2^*}^{x+y}}$ in the function srcdmh with their expansions in $L$ and stores the result in the function dstdmh[DEG_M][eccdeg].

Parameters:
*srcdmhsource function (input)
eccdegdegree in the eccentricity of the function srcdmh.
*dstdmhdestination function (output)

Technical details

The function srcdmh is expanded in $D_2^*$, $a$, $|X|$, $\lambda$ and $h$. It has degree eccdeg in $|X|$.

The function dstdmh[DEG_M][DEG_E] is expanded in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

Creates the function tmpvetdmh[DEG_M] for the terms $\frac{a_1^x a_2^y}{{a_2^*}^{x+y}}$, expanded in $L$.

tmpvetdmh[i] L
perm 1-2
der POL
prod POL
can UND
gen UND
idx HPOL
smin i
smax i

Computes the expansion of $\frac{a_1^x a_2^y}{{a_2^*}^{x+y}}$ by calling the routine init_a2L and stores in the function tmpvetdmh[DEG_M].

Replaces the terms $\frac{a_1^x a_2^y}{{a_2^*}^{x+y}}$ in the function srcdmh with their expansion stores in the function tmpvetdmh[DEG_M], and stores in the function dstdmh[DEG_M][eccdeg].

Definition at line 1067 of file threebody_plane.c.

{
int i, j, l;
DOUBLE coef, tmp;
UNS2 perm1[2], count1[1], type1[1];
INT2 smin1[1], smax1[1];
INT2 ka[2], ktmp[2];
INT2 k5[9];
DMHDR *idxfr, *derfr, *prodfr, *genfr, *canfr;
DMHDR *tmpvetdmh[DEG_M], *tmpdmh;
FN_HANDLERS srch, tmph, dsth;
#if DEBUG
printf("Entry a2L.\n");
#endif
for(i = 0; i < 2; ++i)
perm1[i] = i+1;
count1[0] = 2;
type1[0] = DER_POLYNOM;
derfr = fn_def_fragm(2, perm1, 1, count1, type1);
/* */
type1[0] = PROD_POLYNOM;
prodfr = fn_def_fragm(2, perm1, 1, count1, type1);
/* */
type1[0] = GEN_UNDEF;
genfr = fn_def_fragm(2, perm1, 1, count1, type1);
/* */
type1[0] = CAN_UNDEF;
canfr = fn_def_fragm(2, perm1, 1, count1, type1);
/* */
type1[0] = IDX_HOMOPOL; /* $L$ */
for(i = 0; i < DEG_M; ++i) {
smin1[0] = i; /* $L$ */
smax1[0] = i; /* $L$ */
idxfr = fn_def_idx_fragm(2, perm1, 1, count1, type1, smin1, smax1, NULL, NULL);
tmpvetdmh[i] = fn_create(2, 0, idxfr, derfr, prodfr, genfr, canfr, FN_BIN_TREE, CF_DOUBLE);
dm_release(idxfr);
}
dm_release(derfr);
dm_release(prodfr);
dm_release(genfr);
dm_release(canfr);
get_fn_handlers(srcdmh, &srch);
get_fn_handlers(dstdmh[0][0], &dsth);
get_fn_handlers(tmpvetdmh[0], &tmph);
ka[0] = -1;
ka[1] = -1;
j = srch.first_nz(&coef, k5, srcdmh);
while(j > 0) {
if(ka[0] != k5[1] || ka[1] != k5[2]) {
for(i = 0; i < DEG_M; ++i) {
tmpdmh = fn_duplicate(tmpvetdmh[i]);
fn_remove(tmpvetdmh[i]);
tmpvetdmh[i] = tmpdmh;
}
ka[0] = k5[1];
ka[1] = k5[2];
init_a2L(ka, tmpvetdmh);
}
for(i = 0; i < DEG_M; ++i) {
l = tmph.first_nz(&tmp, ktmp, tmpvetdmh[i]);
while(l > 0) {
tmp *= coef;
k5[1] = ktmp[0];
k5[2] = ktmp[1];
dsth.add_coef(&tmp, k5, dstdmh[i][eccdeg]);
l = tmph.next_nz(&tmp, ktmp, tmpvetdmh[i]);
}
}
j = srch.next_nz(&coef, k5, srcdmh);
}
for(i = 0; i < DEG_M; ++i)
fn_remove(tmpvetdmh[i]);
#if DEBUG
printf("Exit a2L.\n");
#endif
return;
}
void add_kinetic ( DMHDR *  kindmh[DEG_M][DEG_E],
DMHDR *  hamdmh[DEG_M][DEG_E] 
)

This routine adds the kinetic energy, $T$, storesd in the function kindmh[DEG_M][DEG_E], to the Hamiltonian, $H$, storesd in the function hamdmh[DEG_M][DEG_E].

Parameters:
*kindmhkinetic energy, $T$ (input)
*hamdmhHamiltonian, $H$ (output)

Technical details

The functions kindmh[DEG_M][DEG_E] and hamdmh[DEG_M][DEG_E] are expanded in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

Adds to the Hamiltonian, $H$, only the terms up to degree MM_CUT in $\lambda$ contained in the expansion of the kinetic energy, $T$.

Substitutes the fixed value $D_2^*$ in the expansion of the kinetic energy, $T$.

Definition at line 1198 of file threebody_plane.c.

{
int i, j, l;
DOUBLE coef;
INT2 k5[9], trgdeg;
FN_HANDLERS fh;
UNS2 fnstat;
extern DOUBLE D2_fix;
#if DEBUG
printf("Entry add_kinetic.\n");
#endif
get_fn_handlers(hamdmh[0][0], &fh);
for(i = 0; i < DEG_M; ++i) {
for(l = 0; l < DEG_E; ++l) {
fnstat = fn_open_write(hamdmh[i][l]);
j = fh.first_nz(&coef, k5, kindmh[i][l]);
while(j > 0) {
trgdeg = ABS(k5[5])+ABS(k5[6]);
if(trgdeg <= MM_CUT) {
coef *= pow(D2_fix, k5[0]);
k5[0] = 0;
fh.add_coef(&coef, k5, hamdmh[i][l]);
}
j = fh.next_nz(&coef, k5, kindmh[i][l]);
}
if(fnstat != 0) fn_restore_status(hamdmh[i][l], fnstat);
}
}
#if DEBUG
printf("Exit add_kinetic.\n");
#endif
return;
}
void compute_kinetic ( DMHDR *  kindmh[DEG_M][DEG_E],
DMHDR *  wdmh[DEG_E] 
)

This routine computes the expansion of the kinetic energy, $T$, in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

Parameters:
*kindmhkinetic energy, $T$ (output)
*wdmhgenerator $W(f,e)$ (input)

Technical details

\begin{align*} T=-\frac{\beta_1 n_1 a_1 \beta_2 n_2 a_2}{m_0}\frac{1}{\sqrt{1-e_1^2}\sqrt{1-e_2^2}} \Bigl(&(\sin(f_1-h_1)+e_1\sin(-h_1)) (\sin(f_2-h_2)+e_2\sin(-h_2))\\ &+(\cos(f_1-h_1)+e_1\cos(-h_1)) (\cos(f_2-h_2)+e_2\cos(-h_2))\Bigr) \end{align*}

The function kindmh[DEG_M][DEG_E] is expanded in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

The function wdmh[DEG_E] is expanded in $f$. To keep track of the expansion in $e$, we use the index in the array of the function.

Creates the function kinangdmh[DEG_M][DEG_E], for the term

\begin{align*} T_{ang}=&(\sin(f_1-h_1)+e_1\sin(-h_1)) (\sin(f_2-h_2)+e_2\sin(-h_2))\\ &+(\cos(f_1-h_1)+e_1\cos(-h_1)) (\cos(f_2-h_2)+e_2\cos(-h_2))\,, \end{align*}

expanded in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

kinangdmh[i][j] D2* L X lam h
perm 7 8-9 1-2 5-6 3-4
der POL POL POL FOU FOU
prod POL POL POL FOU FOU
can UND UND UND UND UND
gen UND UND UND UND UND
idx HPOL HPOL DAL1 FOU DAL2
smin 0 i j 0 0
smax 0/1 i j j+2 j

Computes $T_{ang}$ by calling the routine kin_ang and stores in the function kinangdmh[DEG_M][DEG_E].

Creates the function kineccdmh[DEG_E], for the term

\[ T_{ecc}=\frac{1}{\sqrt{1-e_1^2}\sqrt{1-e_2^2}}\,, \]

expanded in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

kineccdmh[i] D2* L X lam h
perm 7 8-9 1-2 5-6 3-4
der POL POL POL FOU FOU
prod POL POL POL FOU FOU
can UND UND UND UND UND
gen UND UND UND UND UND
idx HPOL HPOL DAL1 FOU DAL2
smin 0 0 i 0 0
smax 0 0 i 0 i

Computes $T_{ecc}$ by calling the routine kin_ecc and stores in the function kineccdmh[DEG_E].

Creates the function kinLdmh[DEG_M], for the term

\[ T_{L}=-\frac{\beta_1 n_1 a_1 \beta_2 n_2 a_2}{m_0}\,, \]

expanded in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

kinLdmh[i] D2* L X lam h
perm 7 8-9 1-2 5-6 3-4
der POL POL POL FOU FOU
prod POL POL POL FOU FOU
can UND UND UND UND UND
gen UND UND UND UND UND
idx HPOL HPOL DAL1 FOU DAL2
smin 0 i 0 0 0
smax 0 i 0 0 0

Computes $T_{L}$ by calling kin_fastact and stores in the function kinLdmh[DEG_M].

Computes $T=T_{ang}\cdot T_{ecc}\cdot T_{L}$ and stores in the function kindmh[DEG_M][DEG_E].

Definition at line 1265 of file threebody_plane.c.

{
int i, j, l, m;
INT2 smin5[5], smax5[5];
DMHDR *kinangdmh[DEG_M][DEG_E];
DMHDR *kineccdmh[DEG_E];
DMHDR *kinLdmh[DEG_M];
DMHDR *tmpdmh;
#if DEBUG
printf("Entry compute_kinetic.\n");
#endif
smin5[1] = 0; /* $\lambda$ */
smin5[2] = smax5[2] = (D2_FLAG < 0.e0) ? 1 : 0; /* $D_2^*$ */
for(i = 0; i < DEG_M; ++i) {
smin5[3] = i; /* $L$ */
smax5[3] = i; /* $L$ */
for(j = 0; j < DEG_E; ++j) {
smin5[0] = j; /* $|X|$, $h$ */
smax5[0] = j; /* $|X|$, $h$ */
smax5[1] = j+2; /* $\lambda$ */
kinangdmh[i][j] = fn_clone(kindmh[0][0], 4, smin5, smax5, NULL, NULL);
}
}
kin_ang(kinangdmh, wdmh);
for(i = 0; i < DEG_M; ++i)
for(j = 0; j < DEG_E; ++j)
fn_trim(kinangdmh[i][j]);
smin5[1] = 0; /* $\lambda$ */
smax5[1] = 0; /* $\lambda$ */
smin5[2] = 0; /* $D_2^*$ */
smax5[2] = 0; /* $D_2^*$ */
smin5[3] = 0; /* $L$ */
smax5[3] = 0; /* $L$ */
for(i = 0; i < DEG_E; ++i) {
smin5[0] = i; /* $|X|$, $h$ */
smax5[0] = i; /* $|X|$, $h$ */
kineccdmh[i] = fn_clone(kindmh[0][0], 4, smin5, smax5, NULL, NULL);
}
kin_ecc(kineccdmh);
for(i = 0; i < DEG_E; ++i)
fn_trim(kineccdmh[i]);
smin5[0] = 0; /* $|X|$, $h$ */
smax5[0] = 0; /* $|X|$, $h$ */
smin5[1] = 0; /* $\lambda$ */
smax5[1] = 0; /* $\lambda$ */
smin5[2] = 0; /* $D_2^*$ */
smax5[2] = 0; /* $D_2^*$ */
for(i = 0; i < DEG_M; ++i) {
smin5[3] = i; /* $L$ */
smax5[3] = i; /* $L$ */
kinLdmh[i] = fn_clone(kindmh[0][0], 4, smin5, smax5, NULL, NULL);
}
kin_fastact(kinLdmh);
for(i = 0; i < DEG_M; ++i)
fn_trim(kinLdmh[i]);
for(i = 0; i < DEG_E; ++i) {
for(j = 0; j < DEG_E-i; ++j) {
for(l = 0; l < DEG_M; ++l) {
tmpdmh = prod(kinangdmh[l][i], kineccdmh[j]);
for(m = 0; m < DEG_M-l; ++m)
add_prod_to(tmpdmh, kinLdmh[m], kindmh[l+m][i+j]);
fn_remove(tmpdmh);
}
}
}
for(i = 0; i < DEG_M; ++i)
for(l = 0; l < DEG_E; ++l)
fn_trim(kindmh[i][l]);
for(i = 0; i < DEG_M; ++i)
for(l = 0; l < DEG_E; ++l)
fn_remove(kinangdmh[i][l]);
for(i = 0; i < DEG_E; ++i)
fn_remove(kineccdmh[i]);
for(i = 0; i < DEG_M; ++i)
fn_remove(kinLdmh[i]);
#if DEBUG
printf("Exit compute_kinetic.\n");
#endif
return;
}
void compute_utilities ( FILE *  ifp,
FILE *  ofp1,
FILE *  ofp2 
)

This routine computes the initial conditions in Poincare' canonical heliocentric variables and some useful data.

Parameters:
*ifpASCII file with the initial orbital parameters (input)
*ofp1ASCII file for the useful data (output)
*ofp2ASCII file for the Poincare' initial conditions (output)

Technical details

Reads from the file *ifp the initial orbital parameters: $a$, $e$, $l$ and $\omega$.

Computes the initial conditions in Poincare' canonical heliocentric variables, $(\Lambda, \lambda, \xi, \eta)$, the fixed fast actions, $\Lambda^*$, and the translated fast actions, $L=\Lambda-\Lambda^*$.

Stores some useful data in the ASCII file *ofp1.

Stores the initial conditions in Poincare' canonical heliocentric variables, followed by the fixed fast actions, in the ASCII file *ofp2.

Definition at line 1461 of file threebody_plane.c.

{
int i;
DOUBLE coef;
DOUBLE a[2], l[2], e[2], om[2];
DOUBLE poinc_var[8];
char line[26];
extern DOUBLE M0, M1, M2, A_fix[2];
extern DOUBLE Mu1, Mu2, Beta1, Beta2;
extern DOUBLE Alpha_fix, Lambda_fix[2];
#if DEBUG
printf("Entry compute_utilities.\n");
#endif
for(i = 0; fgets(line, 26, ifp) != NULL; ++i) {
sscanf(line, "%lf", &coef);
if(i < 2) a[i%2] = coef; /* $a$ */
else if(i < 4) l[i%2] = coef; /* $e$ */
else if(i < 6) e[i%2] = coef; /* $e$ */
else om[i%2] = coef; /* $\lambda$ */
}
Lambda_fix[0] = Beta1*sqrt(Mu1*A_fix[0]);
Lambda_fix[1] = Beta2*sqrt(Mu2*A_fix[1]);
poinc_var[0] = Beta1*sqrt(Mu1*a[0]);
poinc_var[1] = Beta2*sqrt(Mu2*a[1]);
for(i = 0; i < 2; ++i) {
poinc_var[2+i] = fmod(l[i]+om[i], num.dupig);
if(poinc_var[2+i] < 0.e0) poinc_var[2+i] += num.dupig;
poinc_var[4+i] = sqrt(2.e0*poinc_var[i]*(1.e0-sqrt(1.e0-e[i]*e[i])));
poinc_var[6+i] = -poinc_var[4+i]*sin(om[i]);
poinc_var[4+i] *= cos(om[i]);
poinc_var[i] -= Lambda_fix[i];
}
fprintf(ofp1, "Mass of the star:\n");
fprintf(ofp1, "m_0 = %24.16le\n\n", M0);
fprintf(ofp1, "Mass of the planets:\n");
fprintf(ofp1, "m_1 = %24.16le\n", M1);
fprintf(ofp1, "m_2 = %24.16le\n\n", M2);
fprintf(ofp1, "Mu and Beta:\n");
fprintf(ofp1, "mu_1 = %24.16le\n", Mu1);
fprintf(ofp1, "mu_2 = %24.16le\n", Mu2);
fprintf(ofp1, "beta_1 = %24.16le\n", Beta1);
fprintf(ofp1, "beta_2 = %24.16le\n\n", Beta2);
fprintf(ofp1, "Fixed semi-major axes:\n");
fprintf(ofp1, "a_1^* = %24.16le\n", A_fix[0]);
fprintf(ofp1, "a_2^* = %24.16le\n\n", A_fix[1]);
fprintf(ofp1, "Ratio of the fixed semi-major axes:\n");
fprintf(ofp1, "alpha^* = %24.16le\n\n", Alpha_fix);
fprintf(ofp1, "Fixed fast actions:\n");
fprintf(ofp1, "Lambda_1^* =%24.16le\n", Lambda_fix[0]);
fprintf(ofp1, "Lambda_2^* =%24.16le\n\n", Lambda_fix[1]);
fprintf(ofp1, "Initial Poincare' canonical variables:\n");
fprintf(ofp1, "L_1 = %24.16le\n", poinc_var[0]);
fprintf(ofp1, "L_2 = %24.16le\n", poinc_var[1]);
fprintf(ofp1, "lambda_1 = %24.16le\n", poinc_var[2]);
fprintf(ofp1, "lambda_2 = %24.16le\n", poinc_var[3]);
fprintf(ofp1, "xi_1 = %24.16le\n", poinc_var[4]);
fprintf(ofp1, "xi_2 = %24.16le\n", poinc_var[5]);
fprintf(ofp1, "eta_1 = %24.16le\n", poinc_var[6]);
fprintf(ofp1, "eta_2 = %24.16le\n", poinc_var[7]);
for(i = 0; i < 8; ++i)
fprintf(ofp2, "%24.16le\n", poinc_var[i]);
fprintf(ofp2, "%24.16le\n", Lambda_fix[0]);
fprintf(ofp2, "%24.16le\n", Lambda_fix[1]);
#if DEBUG
printf("Exit compute_utilities.\n");
#endif
return;
}
void cuteps ( DMHDR *  fndmh,
DOUBLE  eps 
)

This routine set to zero all the coefficients of the function fndmh with absolute value smaller than the threshold eps.

Parameters:
*fndmhfunction to be cut (input/output)
epsthreshold

Technical details

The function fndmh is expanded in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

Definition at line 1561 of file threebody_plane.c.

{
int j;
DOUBLE coef, zero = 0.0e0;
INT2 k5[9];
FN_HANDLERS fh;
#if DEBUG
printf("Entry cuteps.\n");
#endif
get_fn_handlers(fndmh, &fh);
j = fh.first_nz(&coef, k5, fndmh);
while(j > 0) {
if((ABS(coef)) < eps) fh.put_coef(&zero, k5, fndmh);
j = fh.next_nz(&coef, k5, fndmh);
}
#if DEBUG
printf("Exit cuteps.\n");
#endif
return;
}
void cutepsrel ( DMHDR *  fndmh,
DOUBLE  eps 
)

This routine set to zero all the coefficients of the function fndmh with absolute value smaller than the threshold eps times the absolute value of the maximum coefficient of the function.

Parameters:
*fndmhfunction to be cut (input/output)
epsrelative threshold

Technical details

The function fndmh is expanded in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

The cut is done by calling the routine cuteps.

Definition at line 1602 of file threebody_plane.c.

{
int j;
DOUBLE coef, abscoef, maxcoef = 0.0e0;
INT2 k5[9];
FN_HANDLERS fh;
#if DEBUG
printf("Entry cutepsrel.\n");
#endif
get_fn_handlers(fndmh, &fh);
j = fh.first_nz(&coef, k5, fndmh);
while(j > 0) {
abscoef = ABS(coef);
if(abscoef > maxcoef) maxcoef = abscoef;
j = fh.next_nz(&coef, k5, fndmh);
}
eps = ABS(eps*maxcoef);
cuteps(fndmh, eps);
#if DEBUG
printf("Exit cutepsrel.\n");
#endif
return;
}
void e2X ( DMHDR *  fdmh[DEG_E],
DMHDR *  e2Xdmh[DEG_E][DEG_E] 
)

This routine replaces the terms $e^i$ in the function fdmh[DEG_E] with their expansions in $|X|$ (storesd in the function e2Xdmh[i][DEG_E]) and stores again in the function fdmh[DEG_E].

Parameters:
*fdmhfunction to be transformed (input/output)
*e2Xdmhexpansion of the powers of $e$ (input)

Technical details

The function fdmh[DEG_E] is originally expanded in $e$ and $(\lambda, h)$. At the end is expanded in $|X|$ and $(\lambda, h)$.

The function e2Xdmh[DEG_E][DEG_E] is expanded in $|X|$.

Definition at line 1649 of file threebody_plane.c.

{
int i, j, l;
DOUBLE coef, coeftmp;
INT2 k2[3], k2tmp[3], kecc;
DMHDR *tmpdmh[DEG_E];
FN_HANDLERS fh, ecch;
#if DEBUG
printf("Entry e2X.\n");
#endif
get_fn_handlers(fdmh[0], &fh);
get_fn_handlers(e2Xdmh[0][0], &ecch);
for(i = 0; i < DEG_E; ++i) {
tmpdmh[i] = fdmh[i];
fdmh[i] = fn_duplicate(tmpdmh[i]);
}
for(i = 0; i < DEG_E; ++i) {
j = fh.first_nz(&coef, k2, tmpdmh[i]);
while(j > 0) {
kecc = k2[0];
for(k2[0] = 0; k2[0] < DEG_E; ++k2[0]) {
l = ecch.first_nz(&coeftmp, k2tmp, e2Xdmh[kecc][k2[0]]);
while(l > 0) {
coeftmp *= coef;
fh.add_coef(&coeftmp, k2, fdmh[k2[0]]);
l = ecch.next_nz(&coeftmp, k2tmp, e2Xdmh[kecc][k2[0]]);
}
}
j = fh.next_nz(&coef, k2, tmpdmh[i]);
}
fn_remove(tmpdmh[i]);
}
#if DEBUG
printf("Exit e2X.\n");
#endif
return;
}
DOUBLE fhyp ( double  a,
double  b,
double  c,
DOUBLE  x,
DOUBLE  eps 
)

This routine computes the hypergeometric function $F(a,b,c;x)$.

Parameters:
aparameter
bparameter
cparameter
xpoint of evaluation
epsthreshold
Return values:
sum$F(a,b,c;x)$

Technical details

\[ F(a,b,c;x) = \sum_{i\geq0} \frac{(a)_i(b)_i}{(c)_i} \frac{x^i}{(1)_i}\,. \]

Definition at line 1715 of file threebody_plane.c.

{
int i;
DOUBLE sum = 1.0e0, tmp = 1.0e0;
#if DEBUG
printf("Entry fhyp.\n");
#endif
for(i = 0; i < 100 && ABS(tmp) > eps; ++i) {
tmp *= (a+i)*(b+i)/((c+i)*(1.e0+i))*x;
sum += tmp;
}
if(i == 100) {
printf("!!! ERROR: the expansion is not convergent!\n");
printf("!!! Last term, F(%le %le %le %16.8le %16.8le): %16.8le\n", a, b, c, x, eps, tmp);
}
#if DEBUG
printf("Exit fhyp.\n");
#endif
return(sum);
}
INT2 flagsincos ( INT2 *  k,
int  n 
)

This routine checks if the array k[n] represents a sinus or a cosinus.

Attention:
No check of any type on the array is made!
Parameters:
*ktrigonometric array
nlength of the array
Return values:
1cosinus
01
-1sinus

Definition at line 1758 of file threebody_plane.c.

{
int i;
INT2 flag = 0;
#if DEBUG
printf("Entry flagsincos.\n");
#endif
for(i = 0; i < n && flag == 0; ++i) {
if(k[i] > 0) flag = 1;
else if(k[i] < 0) flag = -1;
}
#if DEBUG
printf("Exit flagsincos.\n");
#endif
return(flag);
}
void init_a2L ( INT2  ka[2],
DMHDR *  fndmh[DEG_M] 
)

This routine computes the expansion of $\frac{a_1^{{\tt ka[0]}} a_2^{{\tt ka[1]}}}{{a_2^*}^{{\tt ka[0]+ka[1]}}}$ in $L$.

Parameters:
kaarray of exponents
*fndmhexpansion of $\frac{a_1^{{\tt ka[0]}} a_2^{{\tt ka[1]}}}{{a_2^*}^{{\tt ka[0]+ka[1]}}}$ (output)

Technical details

\begin{align*} \frac{a_1^{{\tt ka[0]}}}{{a_2^*}^{{\tt ka[0]}}} \frac{a_2^{{\tt ka[1]}}}{{a_2^*}^{{\tt ka[1]}}}=& \left( \frac{a_1^*}{a_2^*} \right)^{{\tt ka[0]}} \left( 1+\frac{L_1}{\Lambda_1^*} \right)^{2 {\tt ka[0]}} \left( 1+\frac{L_2}{\Lambda_2^*} \right)^{2 {\tt ka[1]}}\\ &=\sum_{i=0}^{{\tt DEG\_M}}{\tt fndmh[i]}\,. \end{align*}

The function fndmh[DEG_M] is expanded in $L$.

Definition at line 1790 of file threebody_plane.c.

{
int i, j, jmax1, jmax2;
DOUBLE coef, tmp1[DEG_M], tmp2[DEG_M];
INT2 k1[2];
FN_HANDLERS fh;
extern DOUBLE Alpha_fix, Lambda_fix[2];
#if DEBUG
printf("Entry init_a2L.\n");
#endif
get_fn_handlers(fndmh[0], &fh);
tmp1[0] = pow(Alpha_fix, ka[0]);
jmax1 = MIN(DEG_M, 2*ka[0]+1);
for(i = 1; i < jmax1; ++i)
tmp1[i] = tmp1[i-1]*(2.e0*ka[0]-i+1.e0)/(i*Lambda_fix[0]);
tmp2[0] = 1.e0;
jmax2 = MIN(DEG_M, 2*ka[1]+1);
for(i = 1; i < jmax2; ++i)
tmp2[i] = tmp2[i-1]*(2.e0*ka[1]-i+1.e0)/(i*Lambda_fix[1]);
for(i = 0; i < jmax1; ++i) {
for(j = 0; j < jmax2 && j < DEG_M-i; ++j) {
coef = tmp1[i]*tmp2[j];
k1[0] = i;
k1[1] = j;
fh.add_coef(&coef, k1, fndmh[i+j]);
}
}
#if DEBUG
printf("Exit init_a2L.\n");
#endif
return;
}
void init_alpha ( DMHDR *  alpha2dmh[DEG_M],
DMHDR *  alphajdmh[(3 *(DEG_E-1)+MM_CUT)/2+1][DEG_M] 
)

This routine computes the expansion of $\alpha^2$ and $\alpha^j$ in $L$.

Parameters:
*alpha2dmhexpansion of $\alpha^2$ (output)
*alphajdmhexpansion of $\alpha^j$ for $j=0$ to $\frac{3({\tt DEG\_E}-1)+{\tt MM\_CUT}}{2}$ (output)

Technical details:

The functions alpha2dmh[DEG_M] and alphajdmh[(3*(DEG_E-1)+MM_CUT)/2+1][DEG_M] are expanded in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

Computes $\alpha = \frac{a_1^*}{a_2^*} \frac{a_1}{a_1^*} \frac{a_2^*}{a_2}$ and stores in the function alphadmh[DEG_M].

Computes $\alpha^2$ and stores in the function alpha2dmh[DEG_M].

Computes $\alpha^j$ and stores in the function alphajdmh[j][DEG_M], for $j=0$ to $\frac{3({\tt DEG\_E}-1)+{\tt MM\_CUT}}{2}$.

Definition at line 1854 of file threebody_plane.c.

{
int i, j, l;
DOUBLE Lambda1[DEG_M], Lambda2[DEG_M], coef;
INT2 k5[9];
DMHDR *alphadmh[DEG_M];
FN_HANDLERS fh;
extern DOUBLE Alpha_fix, Lambda_fix[2];
#if DEBUG
printf("Entry init_alpha.\n");
#endif
for(i = 0; i < DEG_M; ++i)
alphadmh[i] = fn_duplicate(alpha2dmh[i]);
get_fn_handlers(alpha2dmh[0], &fh);
Lambda1[0] = Alpha_fix;
if(DEG_M > 1) Lambda1[1] = 2.e0*Alpha_fix/Lambda_fix[0];
if(DEG_M > 2) Lambda1[2] = Alpha_fix/(Lambda_fix[0]*Lambda_fix[0]);
Lambda2[0] = coef = 1.0e0;
for(i = 1; i < DEG_M; ++i) {
coef *= -Lambda_fix[1];
Lambda2[i] = (i+1.e0)/coef;
}
for(i = 0; i < 9; ++i)
k5[i] = 0;
for(i = 0; i < MIN(DEG_M, 3); ++i) {
k5[1] = i;
for(j = 0; j < DEG_M-i; ++j) {
k5[2] = j;
coef = Lambda1[i]*Lambda2[j];
fh.add_coef(&coef, k5, alphadmh[i+j]);
}
}
for(i = 0; i < DEG_M; ++i) {
for(j = 0; j < DEG_M-i; ++j) {
add_prod_to(alphadmh[i], alphadmh[j], alpha2dmh[i+j]);
}
}
for(i = 0; i < 9; ++i)
k5[i] = 0;
coef = 1.e0;
fh.put_coef(&coef, k5, alphajdmh[0][0]);
for(j = 1; j <= (3*(DEG_E-1)+MM_CUT)/2; ++j) {
for(i = 0; i < DEG_M; ++i) {
for(l = 0; l < DEG_M-i; ++l) {
add_prod_to(alphajdmh[j-1][i], alphadmh[l], alphajdmh[j][i+l]);
}
}
}
for(i = 0; i < DEG_M; ++i)
fn_remove(alphadmh[i]);
#if DEBUG
printf("Exit init_alpha.\n");
#endif
return;
}
void init_bij ( int  i,
int  j,
DMHDR *  alpha2dmh[DEG_M],
DMHDR *  alphajdmh[DEG_M],
DMHDR *  lapldmh[DEG_M] 
)

This routine computes the expansion of the Laplace coefficient $b_{i+\frac{1}{2}}^{(j)}(\alpha)$ in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$. If $j=0$, computes $\frac{1}{2}b_{i+\frac{1}{2}}^{(0)}(\alpha)$.

Parameters:
isubscript of the Laplace coefficient
jsuperscript of the Laplace coefficient
*alpha2dmhexpansion of $\alpha^2$ (input)
*alphajdmhexpansion of $\alpha^j$ (input)
*lapldmhexpansion of $b_{i+\frac{1}{2}}^{(j)}(\alpha)$ (output)

Technical details

\begin{align*} \frac{1}{2} b_{i+\frac{1}{2}}^{(j)}(\alpha) &= \frac{(i+\frac{1}{2})_{j}}{(1)_{j}} \alpha^{j} F(i+\frac{1}{2},i+\frac{1}{2}+j,j+1;\alpha^2)\\ &= \sum_{k=0}^{{\tt DEG\_M}-1} {\tt lapldmh[k]}\,. \end{align*}

The functions alpha2dmh[DEG_M], alphajdmh[(3*(DEG_E-1)+MM_CUT)/2+1][DEG_M] and lapldmh[DEG_M] are expanded in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

Computes the expansion of the hypergeometric function $F(i+\frac{1}{2}, i+\frac{1}{2}+j, j+1; \alpha^{2})$ by calling the routine init_fhyp and stores in the function fhypdmh[DEG_M].

Computes the expansion of $\frac{(i+\frac{1}{2})_{j}}{(1)_{j}} \alpha^{j} F(i+\frac{1}{2}, i+\frac{1}{2}+j, j+1; \alpha^{2})$ and stores in the function lapldmh[DEG_M].

Definition at line 1955 of file threebody_plane.c.

{
int l, m;
DOUBLE coef;
DMHDR *fhypdmh[DEG_M];
FN_HANDLERS fh;
#if DEBUG
printf("Entry init_bij.\n");
#endif
get_fn_handlers(lapldmh[0], &fh);
for(l = 0; l < DEG_M; ++l)
fhypdmh[l] = fn_duplicate(lapldmh[l]);
init_fhyp(i+0.5e0, i+j+0.5e0, j+1.0e0, alpha2dmh, 1.e-17, fhypdmh);
for(l = 0; l < DEG_M; ++l) {
for(m = 0; m < DEG_M-l; ++m) {
add_prod_to(alphajdmh[l], fhypdmh[m], lapldmh[l+m]);
}
}
for(l = 0; l < DEG_M; ++l)
fn_remove(fhypdmh[l]);
coef = 1.0e0;
for(l = 0; l < j; ++l)
coef *= (i+0.5e0+l)/(1.e0+l);
if(j != 0) coef *= 2.e0;
for(l = 0; l < DEG_M; ++l)
scalar_multiply(lapldmh[l], &coef);
#if DEBUG
printf("Exit init_bij.\n");
#endif
return;
}
void init_e2X ( DMHDR *  e2Xdmh[DEG_E][DEG_E])

This routine computes the expansions of the powers of the eccentricity in $|X|$ and $(\lambda, h)$.

Parameters:
*e2Xdmhexpansion of the powers of the eccentricity (output)

Technical details

\begin{align*} e &= |X| \left( 1-\frac{|X|^2}{4} \right)^{\frac{1}{2}}\\ &= \sum_{i\geq0} \left( \prod_{j=0}^{i-1}\left(\frac{1}{2}-j\right) \right) \frac{1}{i!}\left(\frac{-1}{4}\right)^i |X|^{2i+1}\\ &= \sum_{i=0}^{{\tt DEG\_E}-1} {\tt e2Xdmh[1][i]}\,. \end{align*}

The function e2Xdmh[DEG_E][DEG_E] is expanded in $|X|$ and $(\lambda, h)$.

Remarks:
The coefficents of the expansion are computed recursively using the factor $(i-2)/(4(i+1))$.

Definition at line 2029 of file threebody_plane.c.

{
int i, j, k;
DOUBLE coef;
INT2 k2[3];
FN_HANDLERS fh;
#if DEBUG
printf("Entry init_e2X.\n");
#endif
get_fn_handlers(e2Xdmh[0][0], &fh);
for(i = 0; i < 3; ++i)
k2[i] = 0;
coef = 1.e0;
fh.put_coef(&coef, k2, e2Xdmh[0][0]);
for(i = 1; i < DEG_E; i += 2) {
k2[0] = i;
fh.put_coef(&coef, k2, e2Xdmh[1][i]);
coef *= (i-2.e0)/(4.e0*(i+1.e0));
}
for(k = 2; k < DEG_E; ++k) {
for(i = 0; i < DEG_E; ++i) {
for(j = 0; j < DEG_E-i; ++j) {
add_prod_to(e2Xdmh[k-1][i], e2Xdmh[1][j], e2Xdmh[k][i+j]);
}
}
}
#if DEBUG
printf("Exit init_e2X.\n");
#endif
return;
}
void init_fhyp ( double  a,
double  b,
double  c,
DMHDR *  xydmh[DEG_M],
DOUBLE  reps,
DMHDR *  fhypdmh[DEG_M] 
)

This routine computes the expansion of the hypergeometric function $F(a,b,c;x+y)$ in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

Parameters:
aparameter
bparameter
cparameter
xydmhexpansion of $x+y$
repsthreshold
fhypdmhexpansion of $F(a,b,c;x+y)$ (output)

Technical details

\begin{align*} F(a,b,c;x+y) =& \sum_{i\geq0} \frac{(a)_i(b)_i}{(c)_i} \frac{y^i}{(1)_i} F(a+i,b+i,c+i;x)\\ =& \sum_{i=0}^{{\tt DEG_M}-1} {\tt fhypdmh[i]}\,,\\ \end{align*}

The function fhypdmh[DEG_M] is expanded in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

Splits the constant term, $x$, and the polynomial part, $y$, of xydmh[DEG_M] and stores in x and ydmh[DEG_M] respectively.

Computes the expansion of $F(a,b,c;x+y)$ by calling the function fhyp for the computation of $F(a+i,b+i,c+i;x)$ and stores in the function fhypdmh[DEG_M].

Definition at line 2100 of file threebody_plane.c.

{
int i, j, l;
DOUBLE x, tmp, fact;
INT2 k5[9];
DMHDR *ydmh[DEG_M], *ymdmh[DEG_M], *tmpdmh[DEG_M];
FN_HANDLERS fh;
#if DEBUG
printf("Entry init_fhyp.\n");
#endif
get_fn_handlers(xydmh[0], &fh);
for(i = 0; i < 9; ++i)
k5[i] = 0;
fh.get_coef(&x, k5, xydmh[0]);
ydmh[0] = fn_duplicate(xydmh[0]);
ymdmh[0] = fn_duplicate(xydmh[0]);
for(i = 1; i < DEG_M; ++i) {
ydmh[i] = duplicate(xydmh[i]);
ymdmh[i] = fn_duplicate(xydmh[i]);
}
fact = 1.e0;
fh.put_coef(&fact, k5, ymdmh[0]);
for(i = 0; i < DEG_M; ++i) {
tmp = fact*fhyp(a+i, b+i, c+i, x, reps);
for(j = 0; j < DEG_M; ++j)
add_scalarmultiply_to(ymdmh[j], &tmp, fhypdmh[j]);
if(i < DEG_M-1) {
for(j = 0; j < DEG_M; ++j)
tmpdmh[j] = fn_duplicate(ymdmh[j]);
for(j = 0; j < DEG_M; ++j) {
for(l = 0; l < DEG_M-j; ++l) {
add_prod_to(ymdmh[j], ydmh[l], tmpdmh[j+l]);
}
}
for(j = 0; j < DEG_M; ++j) {
fn_remove(ymdmh[j]);
ymdmh[j] = duplicate(tmpdmh[j]);
fn_remove(tmpdmh[j]);
}
fact *= (a+i)*(b+i)/((c+i)*(1.0e0+i));
}
}
for(i = 0; i < DEG_M; ++i) {
fn_remove(ydmh[i]);
fn_remove(ymdmh[i]);
}
#if DEBUG
printf("Exit init_fhyp.\n");
#endif
return;
}
void init_h0 ( DMHDR *  hamdmh[DEG_M+1])

This routine computes the expansion of the unperturbed Hamiltonian, $H_0$, in $L$, $\lambda$ and $(\xi, \eta)$.

Parameters:
*hamdmhunperturbed Hamiltonian, $H_0$ (output)

Technical details

\[ H_0 = -\frac{\mu_1^2 \beta_1^3}{2 \Lambda_1^2} -\frac{\mu_2^2 \beta_2^3}{2 \Lambda_2^2}\,. \]

The expansion of these terms is easy to compute,

\begin{align*} -\frac{\mu^2 \beta^3}{2 \Lambda^2} &= -\frac{\mu^2 \beta^3}{2 {\Lambda^*}^2}\left(1+\frac{L}{\Lambda^*}\right)^{-2}\\ &= - \frac{m_0 m}{2 a^*} \sum_{i\geq0} (i+1)\left(-\frac{1}{\Lambda^*}\right)^i L^i\\ &= \sum_{i=0}^{{\tt DEG\_M}-1}{\tt hamdmh[i]}\,. \end{align*}

The function ham0dmh[DEG_M+1] is expanded in $L$, $\lambda$ and $(\xi, \eta)$.

Definition at line 2194 of file threebody_plane.c.

{
int i;
DOUBLE coef, tmp, fact;
INT2 k3[8];
FN_HANDLERS fh;
extern DOUBLE A_fix[2], Lambda_fix[2];
#if DEBUG
printf("Entry init_h0.\n");
#endif
get_fn_handlers(hamdmh[0], &fh);
for(i = 0; i < 8; ++i)
k3[i] = 0;
fact = -M0*M1/(2.e0*A_fix[0]);
tmp = 1.0e0;
for(i = 0; i <= DEG_M; ++i) {
coef = fact*(i+1.0e0)*tmp;
tmp /= -Lambda_fix[0];
k3[0] = i;
fh.add_coef(&coef, k3, hamdmh[i]);
}
for(i = 0; i < 8; ++i)
k3[i] = 0;
fact = -M0*M2/(2.e0*A_fix[1]);
tmp = 1.0e0;
for(i = 0; i <= DEG_M; ++i) {
coef = fact*(i+1.0e0)*tmp;
tmp /= -Lambda_fix[1];
k3[1] = i;
fh.add_coef(&coef, k3, hamdmh[i]);
}
#if DEBUG
printf("Exit init_h0.\n");
#endif
return;
}
void init_h0_split ( DMHDR *  ham1dmh[DEG_M+1],
DMHDR *  ham2dmh[DEG_M+1] 
)

This routine computes the expansion of the unperturbed Hamiltonian, $H_0$, in $L$, $\lambda$ and $(\xi, \eta)$, splitting the parts related to the inner and outer planet.

Parameters:
ham1dmhunperturbed Hamiltonian, $H_0$ (inner planet)
ham2dmhunperturbed Hamiltonian, $H_0$ (outer planet)
See also:
init_h0

Definition at line 2266 of file threebody_plane.c.

{
int i;
DOUBLE coef, tmp, fact;
INT2 k3[8];
FN_HANDLERS fh;
extern DOUBLE A_fix[2], Lambda_fix[2];
#if DEBUG
printf("Entry init_h0_split.\n");
#endif
get_fn_handlers(ham1dmh[0], &fh);
for(i = 0; i < 8; ++i)
k3[i] = 0;
fact = -M0*M1/(2.e0*A_fix[0]);
tmp = 1.0e0;
for(i = 0; i <= DEG_M; ++i) {
coef = fact*(i+1.0e0)*tmp;
tmp /= -Lambda_fix[0];
k3[0] = i;
fh.add_coef(&coef, k3, ham1dmh[i]);
}
get_fn_handlers(ham2dmh[0], &fh);
for(i = 0; i < 8; ++i)
k3[i] = 0;
fact = -M0*M2/(2.e0*A_fix[1]);
tmp = 1.0e0;
for(i = 0; i <= DEG_M; ++i) {
coef = fact*(i+1.0e0)*tmp;
tmp /= -Lambda_fix[1];
k3[1] = i;
fh.add_coef(&coef, k3, ham2dmh[i]);
}
#if DEBUG
printf("Exit init_h0_split.\n");
#endif
return;
}
void init_pol2poincsec ( INT2  k[2],
DMHDR *  fndmh[DEG_M] 
)

This routine computes the expansion of $|X_1|^{{\tt k[0]}} |X_2|^{{\tt k[1]}}$ in $L$, $\lambda$ and $(\xi, \eta)$. The expansion is explicit only in $L$.

Parameters:
karray of exponents
*fndmhexpansion of $|X_1|^a |X_2|^b$

Technical details

\begin{align*} |X_1|^{{\tt k[0]}} |X_2|^{{\tt k[1]}} &= \frac{|X_1|^{{\tt k[0]}} \Lambda_1^{{{\tt k[0]}}/2}}{\Lambda_1^{{{\tt k[0]}}/2}} \frac{|X_2|^{{\tt k[1]}} \Lambda_2^{{{\tt k[1]}}/2}}{\Lambda_2^{{{\tt k[1]}}/2}}\\ &= \Bigl( |X_1|^{{\tt k[0]}} \Lambda_1^{{{\tt k[0]}}/2} \Bigr) {\Lambda_1^*}^{-{{\tt k[0]}}/2} \left( 1+\frac{L_1}{\Lambda_1^*} \right)^{-{{\tt k[0]}}/2} \Bigl( |X_2|^{{\tt k[1]}} \Lambda_2^{{\tt k[1]}{/2}} \Bigr) {\Lambda_2^*}^{-{{\tt k[1]}}/2} \left( 1+\frac{L_2}{\Lambda_2^*} \right)^{-{{\tt k[1]}}/2}\\ &=\sum_{i=0}^{{\tt DEG\_M}-1} \sqrt\frac{\xi_1^2+\eta_1^2}{2}^{{\tt k[0]}} \sqrt\frac{\xi_2^2+\eta_2^2}{2}^{{\tt k[1]}} {\tt fndmh[i]} \end{align*}

The function fndmh[DEG_M] is expanded in $L$.

Definition at line 2323 of file threebody_plane.c.

{
int i, j;
DOUBLE coef, tmp1coef[DEG_M], tmp2coef[DEG_M];
INT2 k1[2];
FN_HANDLERS fh;
extern DOUBLE Lambda_fix[2];
#if DEBUG
printf("Entry init_pol2poincsec.\n");
#endif
tmp1coef[0] = pow(Lambda_fix[0], -k[0]/2.0e0);
for(i = 1; i < DEG_M; ++i)
tmp1coef[i] = tmp1coef[i-1]*(-k[0]/2.0e0-i+1.0e0)/(i*Lambda_fix[0]);
tmp2coef[0] = pow(Lambda_fix[1], -k[1]/2.0e0);
for(i = 1; i < DEG_M; ++i)
tmp2coef[i] = tmp2coef[i-1]*(-k[1]/2.0e0-i+1.0e0)/(i*Lambda_fix[1]);
get_fn_handlers(fndmh[0], &fh);
for(i = 0; i < DEG_M; ++i) {
for(j = 0; j < DEG_M-i; ++j) {
coef = tmp1coef[i]*tmp2coef[j];
k1[0] = i;
k1[1] = j;
fh.add_coef(&coef, k1, fndmh[i+j]);
}
}
#if DEBUG
printf("Exit init_pol2poincsec.\n");
#endif
return;
}
void init_roveraf ( DMHDR *  roverafdmh[DEG_E])

This routine computes the expansion of $\frac{\|r\|}{a}$ in $e$ and $f$. The expansion is explicit only in $f$.

Parameters:
*roverafdmhexpansion of $\frac{\|r\|}{a}$ (output)

Technical details

\begin{align*} \frac{\|r\|}{a} &= \frac{1-e^2}{1+e\cos(f)}\\ &= 1+ (-\cos(f))e + (-\sin^2(f))e^2 \sum_{i\geq0}(-\cos(f))^i e^i\,.\\ &= \sum_{i=0}^{{\tt DEG\_E}-1} \frac{{\tt roverafdmh[i]}}{i!}e^{i}\,. \end{align*}

The function roverafdmh[DEG_E] is expanded in $f$. To keep track of the expansion in $e$ we use the index in the array of the function.

Remarks:
The routine tri_lie requires an expansion containing the factorials in the denominators.

Definition at line 2388 of file threebody_plane.c.

{
int i;
DOUBLE coef;
INT2 k[1];
FN_HANDLERS fh;
#if DEBUG
printf("Entry init_roveraf.\n");
#endif
get_fn_handlers(roverafdmh[0], &fh);
for(i = 0; i < DEG_E; ++i) {
if(i == 0) { /* $1$ */
k[0] = 0;
coef = 1.e0;
fh.put_coef(&coef, k, roverafdmh[0]);
}
else if(i == 1) { /* $-\cos(f)$ */
k[0] = 1;
coef = -1.e0;
fh.put_coef(&coef, k, roverafdmh[1]);
}
else if(i == 2) { /* $-\sin^{2}(f) = -0.5+0.5\cos(2f)$ */
k[0] = 0;
coef = -1.e0;
fh.put_coef(&coef, k, roverafdmh[2]);
k[0] = 2;
coef = 1.e0;
fh.put_coef(&coef, k, roverafdmh[2]);
}
else { /* $(-\sin^2(f))(-\cos(f))^i$ */
add_prod_to(roverafdmh[i-1], roverafdmh[1], roverafdmh[i]);
coef = (DOUBLE) i;
scalar_multiply(roverafdmh[i], &coef);
}
}
#if DEBUG
printf("Exit init_roveraf.\n");
#endif
return;
}
void init_W ( DMHDR *  wdmh[DEG_E])

This routine computes the expansion of $W(f,e)$ in $e$ and $f$. The expansion is explicit only in $f$.

Parameters:
*wdmhexpansion of $W(f,e)$ (output)

Technical details

\begin{align*} W(f,e) &= \frac{\partial f}{\partial e} = \frac{4\sin(f) + \sin(2f)e}{2(1-e^2)}\\ &= 2 \sin(f) \sum_{i\geq0} e^{2i} +\frac{1}{2} \sin(2f) \sum_{i\geq0} e^{2i+1}\\ &= \sum_{i=0}^{{\tt DEG\_E}-1} \frac{{\tt wdmh[i]}}{i!}e^{i}\,. \end{align*}

The function wdmh[DEG_E] is expanded in $f$. To keep track of the expansion in $e$ we use the index in the array of the function.

Remarks:
The routine tri_lie requires an expansion containing the factorials in the denominators.

Definition at line 2457 of file threebody_plane.c.

{
int i;
DOUBLE fact;
DOUBLE coef;
INT2 k[1];
FN_HANDLERS wh;
#if DEBUG
printf("Entry init_W.\n");
#endif
get_fn_handlers(wdmh[0], &wh);
fact = 1.0e0;
for(i = 0; i < DEG_E; ++i) {
if(i%2 == 0) { /* $-2\sin(-f)$ */
k[0] = -1;
coef = -fact*2.e0;
}
else { /* $-0.5\sin(-2f)$ */
k[0] = -2;
coef = -fact/2.e0;
}
wh.put_coef(&coef, k, wdmh[i]);
fact *= (DOUBLE) i+1;
}
#if DEBUG
printf("Exit init_W.\n");
#endif
return;
}
void inner2outer ( DMHDR *  f1dmh,
DMHDR *  f2dmh 
)

This routine shifts the indices corresponding to the first planet into the ones of the second planet.

Parameters:
*f1dmhthe source function (input)
*f2dmhthe destination function (output)

Technical details

The function f1dmh and f2dmh are expanded in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

Definition at line 2514 of file threebody_plane.c.

{
int i, j;
DOUBLE coef;
INT2 k5[9];
FN_HANDLERS fh;
#if DEBUG
printf("Entry inner2outer.\n");
#endif
get_fn_handlers(f1dmh, &fh);
j = fh.first_nz(&coef, k5, f1dmh);
while(j > 0) {
for(i = 1; i < 9; i += 2) {
k5[i+1] = k5[i];
k5[i] = 0;
}
fh.put_coef(&coef, k5, f2dmh);
j = fh.next_nz(&coef, k5, f1dmh);
}
#if DEBUG
printf("Exit inner2outer.\n");
#endif
return;
}
void kin_ang ( DMHDR *  kinangdmh[DEG_M][DEG_E],
DMHDR *  wdmh[DEG_E] 
)

This routine computes the expansion of $T_{ang}$ in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

Parameters:
*kinangdmhexpansion of $T_{ang}$ (output)
*wdmhexpansion of $W(f,e)$ (input)

Technical details

\begin{align*} T_{ang}=&\left(\sin(f_1-h_1)+e_1\sin(-h_1)\right) \left(\sin(f_2-h_2)+e_2\sin(-h_2)\right)\\ &+\left(\cos(f_1-h_1)+e_1\cos(-h_1)\right) \left(\cos(f_2-h_2)+e_2\cos(-h_2)\right)\\ =& \sum_{i=0}^{{\tt DEG\_M}-1} \sum_{j=0}^{{\tt DEG\_E}-1}{\tt kinangdmh[i][j]}\,. \end{align*}

The function kinangdmh[DEG_M][DEG_E] is expanded in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

Creates the function e2Xdmh[DEG_E][DEG_E] for the powers of the eccentricity, $e^i$, expanded in $|X|$ and $(\lambda, h)$.

e2Xdmh[i][j] X lam & h
perm 1 2-3
der POL FOU
prod POL FOU
can UND UND
gen UND UND
idx HPOL FOU
smin j 0
smax j 0

Computes the expansions of the powers of the eccentricity, $e^i$, by calling the routine init_e2X and stores in the function e2Xdmh[i][DEG_E].

Creates the functions sinfhdmh[DEG_E] and cosfhdmh[DEG_E] for the terms $\sin(f-h)$ and $\cos(f-h)$ originally expanded in $e$ and $(\lambda, h)$. Later are expanded in $|X|$ and $(\lambda, h)$.

(sin/cos)fhdmh[i] e / X lam & h
perm 1 2-3
der POL FOU
prod POL FOU
can UND UND
gen UND UND
idx HPOL FOU
smin i 0
smax i 2i+3

Computes the expansions of $\sin(f-h)$ and $\cos(f-h)$ by calling the routine sincosfh and stores in the functions sinfhdmh[DEG_E] and cosfhdmh[DEG_E].

Computes the expansion of $\sin(f-h)$ and $\cos(f-h)$ in $|X|$ and $(\lambda, h)$ by calling the routine e2X.

Creates the functions sinfh1dmh[DEG_E], cosfh1dmh[DEG_E], sinfh2dmh[DEG_E] and cosfh2dmh[DEG_E] for the terms $\sin(f_1-h_1)$, $\cos(f_1-h_1)$, $\sin(f_2-h_2)$ and $\cos(f_2-h_2)$, expanded in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

(sin/cos)fh(1/2)dmh[i] D2* L X lam h
perm 7 8-9 1-2 5-6 3-4
der POL POL POL FOU FOU
prod POL POL POL FOU FOU
can UND UND UND UND UND
gen UND UND UND UND UND
idx HPOL HPOL DAL1 FOU DAL2
smin 0 0 i 0 0
smax 0 0 i i+1 i

Computes the expansion of $\sin(f_1-h_1)$ and $\cos(f_1-h_1)$ by calling the routine two2five and stores in the functions sinfh1dmh[DEG_E] and cosfh1dmh[DEG_E].

Computes the expansions of $\sin(f_2-h_2)$ and $\cos(f_2-h_2)$ by calling the routine inner2outer and stores in the functions sinfh2dmh[DEG_E] and cosfh2dmh[DEG_E].

Creates the functions esinhdmh[DEG_E] and ecoshdmh[DEG_E] for the quantities $e\sin(-h)$ and $e\cos(-h)$, originally expanded in $e$ and $(\lambda, h)$. Later are expanded in $|X|$ and $(\lambda, h)$.

e(sin/cos)hdmh[i] e / X lam & h
perm 1 2-3
der POL FOU
prod POL FOU
can UND UND
gen UND UND
idx HPOL FOU
smin i 0
smax i 1

Computes the expansions of $e\sin(-h)$ and $e\cos(-h)$ and stores in the functions esinhdmh[DEG_E] and ecoshdmh[DEG_E].

Computes the expansion of $e\sin(-h)$ and $e\cos(-h)$ in $|X|$ by calling the routine e2X.

Creates the functions esinh1dmh[DEG_E], ecosh1dmh[DEG_E], esinh2dmh[DEG_E] and ecosh2dmh[DEG_E] for the terms $e_1\sin(-h_1)$, $e_1\cos(-h_1)$, $e_2\sin(-h_2)$ and $e_2\cos(-h_2)$, expanded in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

e(sin/cos)h(1/2)dmh[i] D2* L X lam h
perm 7 8-9 1-2 5-6 3-4
der POL POL POL FOU FOU
prod POL POL POL FOU FOU
can UND UND UND UND UND
gen UND UND UND UND UND
idx HPOL HPOL DAL1 FOU DAL2
smin 0 0 i 0 0
smax 0 0 i 0 i

Computes the expansion of $e_1\sin(-h_1)$ and $e_1\cos(-h_1)$ by calling the routine two2five and stores in esinh1dmh[DEG_E] and ecosh1dmh[DEG_E].

Computes the expansions of $e_2\sin(-h_2)$ and $e_2\cos(-h_2)$ by calling the routine inner2outer and stores in esinh2dmh[DEG_E] and ecosh2dmh[DEG_E].

Compute the expansion $\sin(f_1-h_1) \sin(f_2-h_2)$ and add to the function kinangdmh[0][DEG_E].

Compute the expansion $\sin(f_1-h_1) e_2\sin(-h_2)$ and add to the function kinangdmh[0][DEG_E].

Compute the expansion $e_1\sin(-h_1) \sin(f_2-h_2)$ and add to the function kinangdmh[0][DEG_E].

Compute the expansion $e_1\sin(-h_1) e_2\sin(-h_2)$ and add to the function kinangdmh[0][DEG_E].

Compute the expansion $e_1\sin(-h_1) e_2\sin(-h_2)$ and add to the function kinangdmh[0][DEG_E].

Compute the expansion $\cos(f_1-h_1) e_2\cos(-h_2)$ and add to the function kinangdmh[0][DEG_E].

Compute the expansion $e_1\cos(-h_1) \cos(f_2-h_2)$ and add to the function kinangdmh[0][DEG_E].

Compute the expansion $e_1\cos(-h_1) e_2\cos(-h_2)$ and add to the function kinangdmh[0][DEG_E].

Definition at line 2558 of file threebody_plane.c.

{
int i, j;
DOUBLE coef;
UNS2 perm2[3], count2[2], type2[2];
INT2 smin2[2], smax2[2], k2[3];
INT2 smin5[9], smax5[9];
DMHDR *idxfr, *derfr, *prodfr, *genfr, *canfr;
DMHDR *e2Xdmh[DEG_E][DEG_E];
DMHDR *sinfhdmh[DEG_E], *cosfhdmh[DEG_E];
DMHDR *sinfh1dmh[DEG_E], *cosfh1dmh[DEG_E];
DMHDR *sinfh2dmh[DEG_E], *cosfh2dmh[DEG_E];
DMHDR *ecoshdmh[DEG_E], *esinhdmh[DEG_E];
DMHDR *esinh1dmh[DEG_E], *ecosh1dmh[DEG_E];
DMHDR *esinh2dmh[DEG_E], *ecosh2dmh[DEG_E];
FN_HANDLERS fh;
extern DOUBLE Alpha_fix;
#if DEBUG
printf("Entry kin_ang.\n");
#endif
for(i = 0; i < 3; ++i)
perm2[i] = i+1;
count2[0] = 1;
count2[1] = 2;
type2[0] = DER_POLYNOM;
type2[1] = DER_FOURIER;
derfr = fn_def_fragm(3, perm2, 2, count2, type2);
/* */
type2[0] = PROD_POLYNOM;
type2[1] = PROD_FOURIER;
prodfr = fn_def_fragm(3, perm2, 2, count2, type2);
/* */
count2[0] = 3;
type2[0] = GEN_UNDEF;
genfr = fn_def_fragm(3, perm2, 1, count2, type2);
/* */
type2[0] = CAN_UNDEF;
canfr = fn_def_fragm(3, perm2, 1, count2, type2);
/* */
count2[0] = 1;
count2[1] = 2;
type2[0] = IDX_HOMOPOL; /* $e$ */
type2[1] = IDX_FOURIER; /* $\lambda$, $h$ */
smin2[1] = 0; /* $\lambda$, $h$ */
smax2[1] = 0; /* $\lambda$, $h$ */
for(i = 0; i < DEG_E; ++i) {
for(j = 0; j < DEG_E; ++j) {
smin2[0] = j; /* $e$ */
smax2[0] = j; /* $e$ */
idxfr = fn_def_idx_fragm(3, perm2, 2, count2, type2, smin2, smax2, NULL, NULL);
e2Xdmh[i][j] = fn_create(3, 0, idxfr, derfr, prodfr, genfr, canfr, FN_BIN_TREE, CF_DOUBLE);
dm_release(idxfr);
}
}
dm_release(derfr);
dm_release(prodfr);
dm_release(genfr);
dm_release(canfr);
init_e2X(e2Xdmh);
smin2[1] = 0; /* $\lambda$, $h$ */
for(i = 0; i < DEG_E; ++i) {
smin2[0] = i; /* $e$ */
smax2[0] = i; /* $e$ */
smax2[1] = 2*i+3; /* $\lambda$, $h$ */
sinfhdmh[i] = fn_clone(e2Xdmh[0][0], 2, smin2, smax2, NULL, NULL);
cosfhdmh[i] = fn_duplicate(sinfhdmh[i]);
}
sincosfh(sinfhdmh, cosfhdmh, wdmh);
e2X(sinfhdmh, e2Xdmh);
e2X(cosfhdmh, e2Xdmh);
smin5[1] = 0; /* $\lambda$ */
smin5[2] = 0; /* $D_2^*$ */
smax5[2] = 0; /* $D_2^*$ */
smin5[3] = 0; /* $L$ */
smax5[3] = 0; /* $L$ */
for(i = 0; i < DEG_E; ++i) {
smin5[0] = i; /* $|X|$, $h$ */
smax5[0] = i; /* $|X|$, $h$ */
smax5[1] = i+1; /* $\lambda$ */
sinfh1dmh[i] = fn_clone(kinangdmh[0][0], 4, smin5, smax5, NULL, NULL);
cosfh1dmh[i] = fn_duplicate(sinfh1dmh[i]);
sinfh2dmh[i] = fn_duplicate(sinfh1dmh[i]);
cosfh2dmh[i] = fn_duplicate(sinfh1dmh[i]);
}
for(i = 0; i < DEG_E; ++i) {
two2five(sinfhdmh[i], sinfh1dmh[i]);
two2five(cosfhdmh[i], cosfh1dmh[i]);
fn_remove(sinfhdmh[i]);
fn_remove(cosfhdmh[i]);
inner2outer(sinfh1dmh[i], sinfh2dmh[i]);
inner2outer(cosfh1dmh[i], cosfh2dmh[i]);
}
smin2[1] = 0; /* $\lambda$, $h$ */
for(i = 0; i < DEG_E; ++i) {
smin2[0] = i; /* $e$ */
smax2[0] = i; /* $e$ */
smax2[1] = 1; /* $\lambda$, $h$ */
esinhdmh[i] = fn_clone(cosfhdmh[0], 2, smin2, smax2, NULL, NULL);
ecoshdmh[i] = fn_duplicate(esinhdmh[i]);
}
get_fn_handlers(esinhdmh[1], &fh);
/* $e\sin(-h)$ */
k2[0] = 1;
k2[1] = 0;
k2[2] = -1;
coef = 1.e0;
fh.put_coef(&coef, k2, esinhdmh[1]);
/* $e\cos(-h)$ */
k2[0] = 1;
k2[1] = 0;
k2[2] = 1;
coef = 1.e0;
fh.put_coef(&coef, k2, ecoshdmh[1]);
e2X(esinhdmh, e2Xdmh);
e2X(ecoshdmh, e2Xdmh);
for(i = 0; i < DEG_E; ++i)
for(j = 0; j < DEG_E; ++j)
fn_remove(e2Xdmh[i][j]);
smin5[1] = 0; /* $\lambda$ */
smax5[1] = 0; /* $\lambda$ */
smin5[2] = 0; /* $D_2^*$ */
smax5[2] = 0; /* $D_2^*$ */
smin5[3] = 0; /* $L$ */
smax5[3] = 0; /* $L$ */
for(i = 0; i < DEG_E; ++i) {
smin5[0] = i; /* $|X|$, $h$ */
smax5[0] = i; /* $|X|$, $h$ */
esinh1dmh[i] = fn_clone(kinangdmh[0][0], 4, smin5, smax5, NULL, NULL);
ecosh1dmh[i] = fn_duplicate(esinh1dmh[i]);
esinh2dmh[i] = fn_duplicate(esinh1dmh[i]);
ecosh2dmh[i] = fn_duplicate(esinh1dmh[i]);
}
for(i = 0; i < DEG_E; ++i) {
two2five(esinhdmh[i], esinh1dmh[i]);
two2five(ecoshdmh[i], ecosh1dmh[i]);
fn_remove(esinhdmh[i]);
fn_remove(ecoshdmh[i]);
inner2outer(esinh1dmh[i], esinh2dmh[i]);
inner2outer(ecosh1dmh[i], ecosh2dmh[i]);
}
for(i = 0; i < DEG_E; ++i) {
for(j = 0; j < DEG_E-i; ++j) {
add_prod_to(sinfh1dmh[j], sinfh2dmh[i], kinangdmh[0][i+j]);
}
}
for(i = 0; i < DEG_E; ++i) {
for(j = 0; j < DEG_E-i; ++j) {
add_prod_to(sinfh1dmh[j], esinh2dmh[i], kinangdmh[0][i+j]);
}
}
for(i = 0; i < DEG_E; ++i) {
for(j = 0; j < DEG_E-i; ++j) {
add_prod_to(esinh1dmh[i], sinfh2dmh[j], kinangdmh[0][i+j]);
}
}
for(i = 0; i < DEG_E; ++i) {
for(j = 0; j < DEG_E-i; ++j) {
add_prod_to(esinh1dmh[i], esinh2dmh[j], kinangdmh[0][i+j]);
}
}
for(i = 0; i < DEG_E; ++i) {
for(j = 0; j < DEG_E-i; ++j) {
add_prod_to(cosfh1dmh[j], cosfh2dmh[i], kinangdmh[0][i+j]);
}
}
for(i = 0; i < DEG_E; ++i) {
for(j = 0; j < DEG_E-i; ++j) {
add_prod_to(cosfh1dmh[j], ecosh2dmh[i], kinangdmh[0][i+j]);
}
}
for(i = 0; i < DEG_E; ++i) {
for(j = 0; j < DEG_E-i; ++j) {
add_prod_to(ecosh1dmh[i], cosfh2dmh[j], kinangdmh[0][i+j]);
}
}
for(i = 0; i < DEG_E; ++i) {
for(j = 0; j < DEG_E-i; ++j) {
add_prod_to(ecosh1dmh[i], ecosh2dmh[j], kinangdmh[0][i+j]);
}
}
for(i = 0; i < DEG_E; ++i) {
fn_remove(esinh1dmh[i]);
fn_remove(ecosh1dmh[i]);
fn_remove(sinfh1dmh[i]);
fn_remove(cosfh1dmh[i]);
fn_remove(esinh2dmh[i]);
fn_remove(ecosh2dmh[i]);
fn_remove(sinfh2dmh[i]);
fn_remove(cosfh2dmh[i]);
}
#if DEBUG
printf("Exit kin_ang.\n");
#endif
return;
}
void kin_ecc ( DMHDR *  kineccdmh[DEG_E])

This routine computes the expansion of $T_{ecc}$ in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

Parameters:
*kineccdmhexpansion of $T_{ecc}$ (output)

Technical details

\begin{align*} T_{ecc}=&\sqrt{1-e_1^2}^{-1} \sqrt{1-e_2^2}^{-1}\\ =&\left(1-\frac{|X_1|^2}{2}\right)^{-1}\left(1-\frac{|X_2|^2}{2}\right)^{-1}\\ =& \sum_{i=0}^{{\tt DEG\_E}-1}{\tt kineccdmh[i]}\,. \end{align*}

The function kineccdmh[DEG_E] is expanded in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

Definition at line 2953 of file threebody_plane.c.

{
int i, j;
DOUBLE coef, coef1[DEG_E], coef2[DEG_E];
INT2 k5[9];
FN_HANDLERS fh;
#if DEBUG
printf("Entry kin_ecc.\n");
#endif
coef1[0] = coef2[0] = 1.0e0;
for(i = 2; i < DEG_E; i += 2) {
coef1[i] = 0.5*coef1[i-2];
coef2[i] = coef1[i];
}
get_fn_handlers(kineccdmh[0], &fh);
for(i = 0; i < 9; ++i)
k5[i] = 0;
for(i = 0; i < DEG_E; i += 2) {
for(j = 0; j < DEG_E-i; j += 2) {
k5[3] = i;
k5[4] = j;
coef = coef1[i]*coef2[j];
fh.put_coef(&coef, k5, kineccdmh[i+j]);
}
}
#if DEBUG
printf("Exit kin_ecc.\n");
#endif
return;
}
void kin_fastact ( DMHDR *  kinazveldmh[DEG_M])

This routine computes the expansion of $T_{L}$ in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

Parameters:
*kinazveldmhexpansion of $T_{L}$ (output)

Technical details

\begin{align*} T_{L} =& -\frac{\beta_1 \beta_2}{m_0} (n_1 a_1) (n_2 a_2)\\ =& \sum_{i=0}^{{\tt DEG\_M}-1}{\tt kinazveldmh[i]}\,. \end{align*}

The function kinazveldmh[DEG_M] is expanded in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

Using the relation $n=\sqrt{\mu/a^2}$, we can express the term $na$, as

\begin{align*} na &= \frac{1}{\sqrt{a}} \bigl(\sqrt{n^2 a^3}\bigr)\\ &= \sqrt{\frac{\mu}{a}} = \sqrt{\frac{\mu}{a^*}} \sqrt{\frac{a^*}{a}} = \sqrt{\frac{\mu}{a^*}} \left( 1+\frac{L}{\Lambda^*} \right)\\ &= \sqrt{\frac{\mu}{a^*}} \sum_{i\geq0} \left(-\frac{1}{\Lambda^*}\right)^i L^i\,. \end{align*}

Definition at line 3009 of file threebody_plane.c.

{
int i, j;
DOUBLE coef, coef1[DEG_M], coef2[DEG_M];
INT2 k5[9];
FN_HANDLERS fh;
extern DOUBLE A_fix[2], Lambda_fix[2];
#if DEBUG
printf("Entry kin_fastact.\n");
#endif
get_fn_handlers(kinazveldmh[0], &fh);
coef = Beta1*Beta2/M0*sqrt(Mu1 *Mu2/(A_fix[0]*A_fix[1]));
coef1[0] = coef;
for(i = 1; i < DEG_M; ++i)
coef1[i] = coef1[i-1]/(-Lambda_fix[0]);
coef2[0] = 1.0e0;
for(i = 1; i < DEG_M; ++i)
coef2[i] = coef2[i-1]/(-Lambda_fix[1]);
for(i = 0; i < 9; ++i)
k5[i] = 0;
for(i = 0; i < DEG_M; ++i) {
for(j = 0; j < DEG_M-i; ++j) {
coef = coef1[i]*coef2[j];
k5[1] = i;
k5[2] = j;
fh.add_coef(&coef, k5, kinazveldmh[i+j]);
}
}
#if DEBUG
printf("Exit kin_fastact.\n");
#endif
return;
}
void laplace ( DMHDR *  svilazdmh[DEG_M][DEG_E][DEG_E],
DMHDR *  lapldmh[(3 *(DEG_E-1)+MM_CUT)/2+1][DEG_M],
DMHDR *  fndmh[DEG_M][DEG_E],
int  i,
int  ecc 
)

This routine computes the expansion of the term related to the index i and of degree ecc in the eccentricity of the expansion of the potential energy, in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$..

Parameters:
*svilazdmhexpansion of $\left( \frac{\Xi}{{a_1^*}^2} \right)^{i}$
*lapldmhexpansions of $b_{i+\frac{1}{2}}^{(|n|)}(\alpha)$ for $n=1$ to $\frac{3({\tt DEG\_E}-1)+{\tt MM\_CUT}}{2}$. For $n=0$ $\frac{1}{2} b_{i+\frac{1}{2}}^{(0)}(\alpha)$ (input)
*fndmhexpansion of the term of the potential energy
iindex of the term to be computes
eccdegree in $|X|$

Technical details

\[ -\frac{m_1 m_2}{a_2^*}\frac{(-1)^i (\frac{1}{2})_i}{(1)_i} \left( \frac{\Xi}{{a_2^*}^2} \right)^{i} \left( \frac{a_2^*}{{a_2}} \right)^{2i+1} \left( \sum_{n=-\infty}^{+\infty} \frac{1}{2} b_{i+\frac{1}{2}}^{|n|}(\alpha) \cos(n\lambda_1-n\lambda_2) \right) \]

The functions svilazdmh[DEG_M][DEG_E][DEG_E], lapldmh[(3*(DEG_E-1)+MM_CUT)/2+1][DEG_M] and fndmh[DEG_M][DEG_E] are expanded in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

Multiplies by $-\frac{m_1 m_2}{a_2^*}$ and evaluates the power of $D_2^*$ in the expansion svilazdmh[DEG_M][i][ecc].

\[ \left({\frac{a_2^*}{a_2}}\right)^{2i+1}= \left(1+\frac{L_2}{\Lambda_2}\right)^{-(4i+2)}\,, \]

Remarks:
The coefficents of the expansion are computed recursively using the factor $-\frac{4i+2+m}{(m+1)L_2^*}$, $m$ being the degree in the fast action $L_2$.

Definition at line 3089 of file threebody_plane.c.

{
int j, l, m, n, r, s, sign;
DOUBLE coef, coeftmp, coeflapl;
INT2 k5[9], ktmp[9], klapl[9];
INT2 tmptrgdeg, trgdeg;
INT2 sincosold, sincosnew;
FN_HANDLERS svh, lph, fh;
UNS2 fnstat, fnstatvet[DEG_M];
extern DOUBLE A_fix[2], Lambda_fix[2], D2_fix;
#if DEBUG
printf("Entry laplace.\n");
#endif
get_fn_handlers(svilazdmh[0][i][ecc], &svh);
get_fn_handlers(lapldmh[0][0], &lph);
get_fn_handlers(fndmh[0][ecc], &fh);
for(j = 0; j < DEG_M; ++j)
fnstatvet[j] = fn_open_write(fndmh[j][ecc]);
/* loop over $\left( \frac{\Xi}{{a_1^*}^2} \right)^{i}$ */
for(j = 0; j < DEG_M; ++j) {
fnstat = fn_open_read(svilazdmh[j][i][ecc]);
l = svh.first_nz(&coef, k5, svilazdmh[j][i][ecc]);
while(l > 0) {
coef *= -1.e0*M1*M2/A_fix[1];
coef *= pow(D2_fix, k5[0]);
ktmp[0] = 0;
for(m = 1; m < 9; ++m)
ktmp[m] = k5[m];
sincosold = flagsincos(k5+5, 2);
tmptrgdeg = ABS(k5[5])+ABS(k5[6]);
/* loop over $({a_2^*}/{a_2})^{2i+1}$ */
for(m = 0; m < DEG_M-j; ++m) {
/* loop in $\cos(n\lambda_1-n\lambda_2)$ */
for(n = 0; n <= (tmptrgdeg+MM_CUT)/2; ++n) {
for(sign = -1; sign < 2; sign += 2) {
if(n == 0) {
trgdeg = tmptrgdeg;
coeftmp = coef;
ktmp[5] = k5[5];
ktmp[6] = k5[6];
if(sign > 0) coeftmp = 0.0e0;
}
else if(sincosold == 0) { /* terms independent of $\lambda$ */
trgdeg = n+n;
coeftmp = coef;
ktmp[5] = n;
ktmp[6] = -n;
if(sign > 0) coeftmp = 0.0e0;
}
else {
ktmp[5] = k5[5]+sign*n;
ktmp[6] = k5[6]-sign*n;
coeftmp = 0.5e0*coef;
trgdeg = ABS(ktmp[5])+ABS(ktmp[6]);
sincosnew = flagsincos(ktmp+5, 2);
if(sincosold > 0 && sincosnew < 0) {
ktmp[5] = -ktmp[5];
ktmp[6] = -ktmp[6];
}
else if(sincosold < 0 && sincosnew > 0) {
ktmp[5] = -ktmp[5];
ktmp[6] = -ktmp[6];
coeftmp = -coeftmp;
}
else if(sincosold < 0 && sincosnew == 0) {
coeftmp = 0.0e0;
}
}
if(trgdeg <= MM_CUT) {
/* loop over $1/2\ b_{i+1/2}^{|n|}(\alpha)$ */
for(r = 0; r < DEG_M-j-m; ++r) {
s = lph.first_nz(&coeflapl, klapl, lapldmh[n][r]);
while(s > 0) {
coeflapl *= coeftmp;
ktmp[1] = k5[1]+klapl[1];
ktmp[2] = k5[2]+klapl[2];
fh.add_coef(&coeflapl, ktmp, fndmh[r+j+m][ecc]);
s = lph.next_nz(&coeflapl, klapl, lapldmh[n][r]);
}
}
}
}
}
coef *= (-1.e0*(4*i+2+m))/((m+1)*Lambda_fix[1]);
++k5[2];
}
l = svh.next_nz(&coef, k5, svilazdmh[j][i][ecc]);
}
if(fnstat != 0) fn_restore_status(svilazdmh[j][i][ecc], fnstat);
}
for(j = 0; j < DEG_M; ++j)
if(fnstatvet[j] != 0) fn_restore_status(fndmh[j][ecc], fnstatvet[j]);
#if DEBUG
printf("Exit laplace.\n");
#endif
return;
}
int main ( )

This program computes the expansion of the Hamiltonian of the three-body problem in $L$, $\lambda$ and $(\xi, \eta)$.

Technical details

Reads from the file parameters.dat the values of the masses $(m_0, m_1, m_2)$, in Jupiter mass, and of the fixed semi-major axes $(a_1^*, a_2^*)$, in A.U..

Computes $\mu_1$, $\mu_2$, $\beta_1$, $\beta_2$ and $\alpha^*$.

Computes the initial conditions in $L$, $\lambda$, $\xi$ and $\eta$; and some useful data by calling the routine compute_utilities.

Creates the function wdmh[DEG_E] for the generator $W(f,e)$ expanded in $e$ and $f$. The expansion is explicit only in $f$, to keep track of the expansion in $e$, we use the index in the array of the function.

wdmh[i] f
perm 1
der FOU
prod FOU
can UND
gen UND
idx FOU
smin 0
smax i%2+1

Initializes the generator $W(f,e)$ by calling the routine init_W and stores in the function wdmh[DEG_E].

Creates the function xiovera2dmh[DEG_M][DEG_E] for $\frac{\Xi}{{a_1^*}^2}$ expanded in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

xiovera2dmh[i][j] D2* L X lam h
perm 7 8-9 1-2 5-6 3-4
der POL POL POL FOU FOU
prod POL POL POL FOU FOU
can UND UND UND UND UND
gen UND UND UND UND UND
idx POL HPOL DAL1 FOU DAL2
smin 0 i j 0 0
smax 0/1 i j j+2 j

Computes the expansion of $\frac{\Xi}{{a_1^*}^2}$ by calling the routine xi_over_a2 and stores in the function xiovera2dmh[DEG_M][DEG_E].

Creates the function xiovera2powdmh[DEG_M][0][DEG_E] for the term $\left(\frac{\Xi}{{a_1^*}^2}\right)^0=1$, expanded in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

xiovera2powdmh[i][0][j] D2* L X lam h
perm 7 8-9 1-2 5-6 3-4
der POL POL POL FOU FOU
prod POL POL POL FOU FOU
can UND UND UND UND UND
gen UND UND UND UND UND
idx HPOL HPOL DAL1 FOU DAL2
smin 0 i j 0 0
smax 0 i j j j

Computes the expansion of $\left( \frac{\Xi}{{a_1^*}^2} \right)^0 = 1$ and stores in the function xiovera2powdmh[0][0][0].

Create the functions lapldmh[(3*(DEG_E-1)+MM_CUT)/2][DEG_M] for the Laplace coefficients, $\frac{1}{2} b_{0+\frac{1}{2}}^{(|j|)}$, alpha2dmh[DEG_M] for $\alpha^2$ and alphajdmh[(3*(DEG_E-1)+MM_CUT+1)/2][DEG_M] for $\alpha^j$ ( $j=0$ to $\frac{3{\tt DEG\_E-1}+{\tt MM\_CUT}}{2}$), expanded in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

lapldmh[i][j]/alpha(2/j)dmh([i])[j] D2* L X lam h
perm 7 8-9 1-2 5-6 3-4
der POL POL POL FOU FOU
prod POL POL POL FOU FOU
can UND UND UND UND UND
gen UND UND UND UND UND
idx HPOL HPOL DAL1 FOU DAL2
smin 0 j 0 0 0
smax 0 j 0 0 0

Computes the expansions of the terms $\alpha^2$ and $\alpha^j$ by calling the routine init_alpha and stores in the functions alpha2dmh[DEG_M] and alphajdmh[(3*(DEG_E-1)+MM_CUT)/2][DEG_M].

Computes the Laplace coefficients $\frac{1}{2} b_{0+\frac{1}{2}}^{(|j|)}$ by calling the routine init_bij and stores in the function lapldmh[MM_CUT/2][DEG_M].

Creates the function hamtmpdmh[DEG_M][DEG_E] for the expansion of the Hamiltonian in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

hamtmpdmh[i][j] D2* L X lam h
perm 7 8-9 1-2 5-6 3-4
der POL POL POL FOU FOU
prod POL POL POL FOU FOU
can UND UND UND UND UND
gen UND UND UND UND UND
idx HPOL HPOL DAL1 FOU DAL2
smin 0 i j 0 0
smax 0 i j MM_CUT j

Computes the expansion of

\[ -\frac{m_1 m_2}{a_2^*} \frac{(-1)^0 (\frac{1}{2})_0}{(1)_0} \left( \frac{\Xi}{{a_2^*}^2} \right)^{0} \left( \frac{a_2^*}{{a_2}} \right)^{2\cdot0+1} \left( \sum_{n=-\infty}^{+\infty} \frac{1}{2} b_{0+\frac{1}{2}}^{|n|}(\alpha) \cos(n\lambda_1^*-n\lambda_2^*) \right)\,, \]

by calling the routine laplace and stores in the function hamtmpdmh[DEG_M][DEG_E].

Remarks:
The computation of the potential energy assumes that the parameter DEG_E is bigger than DEG_M.

The term $\left(\frac{\Xi}{{a_1^*}^2}\right)^l$ has degree ${\tt DEG\_E}-1+2l$ in the fast angle, $\lambda$, so we compute the Laplace coefficients $b_{l+\frac{1}{2}}^{|i|} (\alpha)$ with i up to $\frac{{\tt DEG\_E}-1+2l+{\tt MM\_CUT}}{2}$. The limits are due to the fact that the Laplace coefficients must be multiplyed by $\cos(i\lambda_1-i\lambda_2)$, that has trignometric degree 2i, so the terms of higher degree will exceed the limit.

Computes the Laplace coefficients $b_{l+\frac{1}{2}}^{(|i|)}$ by calling the routine init_bij and stores in the function lapldmh[(DEG_E-1+2*l+MM_CUT)/2][DEG_M]

Creates the function xiovera2powdmh[DEG_M][l][DEG_E] for the term $\left(\frac{\Xi}{{a_1^*}^2}\right)^l$, expanded in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

xiovera2powdmh[i][l][j] D2* L X lam h
perm 7 8-9 1-2 5-6 3-4
der POL POL POL FOU FOU
prod POL POL POL FOU FOU
can UND UND UND UND UND
gen UND UND UND UND UND
idx HPOL HPOL DAL1 FOU DAL2
smin 0 i j 0 0
smax 0/MIN(l, (DEG_E-1-j)/2) i j j+2l j

Computes the term $(-1)^l \frac{(\frac{1}{2})_l}{(1)_l} \left(\frac{\Xi}{{a_1^*}^2} \right)^l$ by calling the routine prod_D2 (with the flag set to TRUE)and stores in xiovera2powdmh[DEG_M][l][DEG_E].

Computes the expansion of

\[ -\frac{m_1 m_2}{a_2^*} \frac{(-1)^l (\frac{1}{2})_l}{(1)_l} \left( \frac{\Xi}{{a_2^*}^2} \right)^{l} \left( \frac{a_2^*}{{a_2}} \right)^{2l+1} \left( \sum_{n=-\infty}^{+\infty} \frac{1}{2} b_{l+\frac{1}{2}}^{|n|}(\alpha) \cos(n\lambda_1^*-n\lambda_2^*) \right) \]

by calling the routine laplace and stores in the function hamtmpdmh[DEG_M][DEG_E].

Creates the function kindmh[DEG_M][DEG_E] for the kinetic energy, $T$, expanded in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

kindmh[i][j] D2* L X lam h
perm 7 8-9 1-2 5-6 3-4
der POL POL POL FOU FOU
prod POL POL POL FOU FOU
can UND UND UND UND UND
gen UND UND UND UND UND
idx HPOL HPOL DAL1 FOU DAL2
smin 0 i j 0 0
smax 0/1 i j j+2 j

Computes the kinetic energy, $T$, by calling the routine compute_kinetic and stores in kindmh[DEG_M][DEG_E].

Adds the kinetic energy, $T$, to the Hamiltonian, $H$, by calling the routine add_kinetic and stores in the function hamtmpdmh[DEG_M][DEG_E].

Creates the function ham0dmh[DEG_M+1] for the unperturbed Hamiltonian, $H_0$, expanded in $L$, $\lambda$ and $(\xi, \eta)$.

ham0dmh[i] L lam xi/eta
perm 7-8 5-6 1,3,2,4
der POL FOU POL
prod POL FOU POL
can ACT1 ANG1 GEN
gen UND UND UND
idx HPOL FOU HPOL
smin i 0 0
smax i 0 0

Computes the expansion of $H_0$ by calling the routine init_h0 and stores in the function ham0dmh[DEG_M+1].

Creates the function ham0dmh_1 and ham0dmh_2 for the "single planet" unperturbed Hamiltonian, $H_0$, expanded in $L$, $\lambda$ and $(\xi, \eta)$.

ham0dmh_(1/2)[i] L lam xi/eta
perm 7-8 5-6 1,3,2,4
der POL FOU POL
prod POL FOU POL
can ACT1 ANG1 GEN
gen UND UND UND
idx HPOL FOU HPOL
smin i 0 0
smax i 0 0

Computes the expansions of the "single planet" $H_0$, splitting the inner and outer planet, by calling the routine init_h0_split and stores in the functions ham0dmh_1 and ham0dmh_2.

Create the function ham1dmh[DEG_M][DEG_E] for the perturbation, $H_1$, expanded in $L$, $\lambda$ and $(\xi, \eta)$.

ham1dmh[i][j] L lam xi/eta
perm 7-8 5-6 1,3,2,4
der POL FOU POL
prod POL FOU POL
can ACT1 ANG1 GEN
gen UND UND UND
idx HPOL FOU HPOL
smin i 0 j
smax i MM_CUT j

Transforms the expansion of $H_1$ contained in hamtmp[DEG_M][DEG_E], by calling the routine secular2poinc and stores in the function ham1dmh[DEG_M][DEG_E].

Definition at line 285 of file threebody_plane.c.

{
int i, j, l, m, n;
DOUBLE coef;
UNS2 perm1[1], count1[1], type1[1];
INT2 smin1[1], smax1[1];
UNS2 perm3[8], count3[3], type3[3];
INT2 smin3[3], smax3[3], k3[8];
UNS2 perm5[9], count5[5], type5[5];
INT2 smin5[5], smax5[5], k5[9];
DMHDR *idxfr, *derfr, *prodfr, *genfr, *canfr;
DMHDR *tmpdmh;
DMHDR *wdmh[DEG_E];
DMHDR *alpha2dmh[DEG_M];
DMHDR *alphajdmh[(3*(DEG_E-1)+MM_CUT)/2+1][DEG_M];
DMHDR *xiovera2dmh[DEG_M][DEG_E];
DMHDR *xiovera2powdmh[DEG_M][DEG_E][DEG_E];
DMHDR *lapldmh[(3*(DEG_E-1)+MM_CUT)/2+1][DEG_M];
DMHDR *kindmh[DEG_M][DEG_E];
DMHDR *hamtmpdmh[DEG_M][DEG_E];
DMHDR *ham0dmh[DEG_M+1];
DMHDR *ham0dmh_1[DEG_M+1];
DMHDR *ham0dmh_2[DEG_M+1];
DMHDR *ham1dmh[DEG_M][DEG_E];
FN_HANDLERS fh;
UNS2 fnstat;
char line[26], filename[256];
FILE *ifp;
FILE *ofp1, *ofp2;
extern DOUBLE M0, M1, M2, A_fix[2];
extern DOUBLE Alpha_fix, Mu1, Mu2, Beta1, Beta2;
chr_init();
fn_setflag_sleep_default(0);
fn_set_current_type(FN_BIN_TREE);
ifp = fopen("parameters.dat", "r");
for(i = 0; fgets(line, 26, ifp) != NULL; ++i) {
sscanf(line, "%lf", &coef);
if(i == 0) M0 = MAS_SUN*coef;
else if(i == 1) M1 = MAS_JUP*coef;
else if(i == 2) M2 = MAS_JUP*coef;
else if(i == 3) A_fix[0] = coef;
else if(i == 4) A_fix[1] = coef;
}
fclose(ifp);
Mu1 = M0+M1;
Mu2 = M0+M2;
Beta1 = M0/(M0+M1)*M1;
Beta2 = M0/(M0+M2)*M2;
ifp = fopen("orbitals.dat", "r");
ofp1 = fopen("utilities.dat", "w");
ofp2 = fopen("poinc_ini.dat", "w");
compute_utilities(ifp, ofp1, ofp2);
fclose(ifp);
fclose(ofp1);
fclose(ofp2);
perm1[0] = 1;
count1[0] = 1;
type1[0] = DER_FOURIER;
derfr = fn_def_fragm(1, perm1, 1, count1, type1);
/* */
type1[0] = PROD_FOURIER;
prodfr = fn_def_fragm(1, perm1, 1, count1, type1);
/* */
type1[0] = GEN_UNDEF;
genfr = fn_def_fragm(1, perm1, 1, count1, type1);
/* */
type1[0] = CAN_UNDEF;
canfr = fn_def_fragm(1, perm1, 1, count1, type1);
/* */
type1[0] = IDX_FOURIER; /* $f$ */
smin1[0] = 0; /* $f$ */
for(i = 0; i < DEG_E; ++i) {
smax1[0] = i%2+1;/* $f$ */
idxfr = fn_def_idx_fragm(1, perm1, 1, count1, type1, smin1, smax1, NULL, NULL);
wdmh[i] = fn_create(1, 0, idxfr, derfr, prodfr, genfr, canfr, FN_BIN_TREE, CF_DOUBLE);
dm_release(idxfr);
}
dm_release(prodfr);
dm_release(derfr);
dm_release(genfr);
dm_release(canfr);
init_W(wdmh);
for(i = 0; i < 9; ++i)
perm5[i] = i+1;
count5[0] = 1;
for(i = 1; i < 5; ++i)
count5[i] = 2;
type5[0] = DER_POLYNOM;
type5[1] = DER_POLYNOM;
type5[2] = DER_POLYNOM;
type5[3] = DER_FOURIER;
type5[4] = DER_FOURIER;
derfr = fn_def_fragm(9, perm5, 5, count5, type5);
/* */
type5[0] = PROD_POLYNOM;
type5[1] = PROD_POLYNOM;
type5[2] = PROD_POLYNOM;
type5[3] = PROD_FOURIER;
type5[4] = PROD_FOURIER;
prodfr = fn_def_fragm(9, perm5, 5, count5, type5);
/* */
count5[0] = 9;
type5[0] = GEN_UNDEF;
genfr = fn_def_fragm(9, perm5, 1, count5, type5);
/* */
type5[0] = CAN_UNDEF;
canfr = fn_def_fragm(9, perm5, 1, count5, type5);
/* */
perm5[0] = 7;
perm5[1] = 8;
perm5[2] = 9;
perm5[3] = 1;
perm5[4] = 2;
perm5[5] = 5;
perm5[6] = 6;
perm5[7] = 3;
perm5[8] = 4;
count5[0] = 4;
count5[1] = 2;
count5[2] = 1;
count5[3] = 2;
type5[0] = IDX_DALEMBERT; /* $|X|$, $h$ */
type5[1] = IDX_FOURIER; /* $\lambda$ */
type5[2] = IDX_POLYNOM; /* $D_2^*$ */
type5[3] = IDX_HOMOPOL; /* $L$ */
smin5[1] = 0; /* $\lambda$ */
smin5[2] = 0; /* $D_2^*$ */
smax5[2] = (D2_FLAG < 0.e0) ? 1 : 0; /* $D_2^*$ */
for(i = 0; i < DEG_M; ++i) {
smin5[3] = i; /* $L$ */
smax5[3] = i; /* $L$ */
for(j = 0; j < DEG_E; ++j) {
smin5[0] = j; /* $|X|$, $h$ */
smax5[0] = j; /* $|X|$, $h$ */
smax5[1] = j+2; /* $\lambda$ */
idxfr = fn_def_idx_fragm(9, perm5, 4, count5, type5, smin5, smax5, NULL, NULL);
xiovera2dmh[i][j] = fn_create(9, 0, idxfr, derfr, prodfr, genfr, canfr, FN_BIN_TREE, CF_DOUBLE);
dm_release(idxfr);
}
}
dm_release(derfr);
dm_release(prodfr);
dm_release(genfr);
dm_release(canfr);
xi_over_a2(xiovera2dmh, wdmh);
for(i = 0; i < DEG_M; ++i) {
for(j = 0; j < DEG_E; ++j) {
fn_trim(xiovera2dmh[i][j]);
fn_set_sleep(xiovera2dmh[i][j]);
}
}
smin5[1] = 0; /* $\lambda$ */
smin5[2] = 0; /* $D_2^*$ */
smax5[2] = 0; /* $D_2^*$ */
for(i = 0; i < DEG_M; ++i) {
smin5[3] = i; /* $L$ */
smax5[3] = i; /* $L$ */
for(j = 0; j < DEG_E; ++j) {
smin5[0] = j; /* $|X|$, $h$ */
smax5[0] = j; /* $|X|$, $h$ */
smax5[1] = j; /* $\lambda$ */
xiovera2powdmh[i][0][j] = fn_clone(xiovera2dmh[0][0], 4, smin5, smax5, NULL, NULL);
fn_set_sleep(xiovera2powdmh[i][0][j]);
}
}
fnstat = fn_open_write(xiovera2powdmh[0][0][0]);
get_fn_handlers(xiovera2powdmh[0][0][0], &fh);
for(i = 0; i < 9; ++i)
k5[i] = 0;
coef = 1.e0;
fh.put_coef(&coef, k5, xiovera2powdmh[0][0][0]);
if(fnstat != 0) fn_restore_status(xiovera2powdmh[0][0][0], fnstat);
for(i = 0; i < DEG_M; ++i) {
for(j = 0; j < DEG_E; ++j) {
fn_trim(xiovera2powdmh[i][0][j]);
}
}
smin5[0] = 0; /* $|X|$, $h$ */
smax5[0] = 0; /* $|X|$, $h$ */
smin5[1] = 0; /* $\lambda$ */
smax5[1] = 0; /* $\lambda$ */
smin5[2] = 0; /* $D_2^*$ */
smax5[2] = 0; /* $D_2^*$ */
for(i = 0; i <= (3*(DEG_E-1)+MM_CUT)/2; ++i) {
for(j = 0; j < DEG_M; ++j) {
smin5[3] = j; /* $L$ */
smax5[3] = j; /* $L$ */
lapldmh[i][j] = fn_clone(xiovera2dmh[0][0], 4, smin5, smax5, NULL, NULL);
}
}
for(i = 0; i < DEG_M; ++i) {
alpha2dmh[i] = fn_duplicate(lapldmh[0][i]);
}
for(j = 0; j <= (3*(DEG_E-1)+MM_CUT)/2; ++j) {
for(i = 0; i < DEG_M; ++i) {
alphajdmh[j][i] = fn_duplicate(lapldmh[0][i]);
}
}
init_alpha(alpha2dmh, alphajdmh);
for(j = 0; j <= MM_CUT/2; ++j)
init_bij(0, j, alpha2dmh, alphajdmh[j], lapldmh[j]);
smin5[1] = 0; /* $\lambda$ */
smax5[1] = MM_CUT; /* $\lambda$ */
smin5[2] = 0; /* $D_2^*$ */
smax5[2] = 0; /* $D_2^*$ */
for(i = 0; i < DEG_M; ++i) {
smin5[3] = i; /* $L$ */
smax5[3] = i; /* $L$ */
for(j = 0; j < DEG_E; ++j) {
smin5[0] = j; /* $|X|$, $h$ */
smax5[0] = j; /* $|X|$, $h$ */
hamtmpdmh[i][j] = fn_clone(xiovera2dmh[0][0], 4, smin5, smax5, NULL, NULL);
fn_set_sleep(hamtmpdmh[i][j]);
}
}
for(i = 0; i < DEG_E; ++i)
laplace(xiovera2powdmh, lapldmh, hamtmpdmh, 0, i);
for(l = 1; l < DEG_E; ++l) {
for(i = 0; i <= (DEG_E-1+2*l+MM_CUT)/2; ++i) {
for(j = 0; j < DEG_M; ++j) {
tmpdmh = fn_duplicate(lapldmh[i][j]);
fn_remove(lapldmh[i][j]);
lapldmh[i][j] = fn_duplicate(tmpdmh);
fn_remove(tmpdmh);
}
}
for(i = 0; i <= (DEG_E-1+2*l+MM_CUT)/2; ++i)
init_bij(l, i, alpha2dmh, alphajdmh[i], lapldmh[i]);
smin5[1] = 0; /* $\lambda$ */
for(i = 0; i < DEG_M; ++i) {
smin5[3] = i; /* $L$ */
smax5[3] = i; /* $L$ */
for(j = 0; j < DEG_E; ++j) {
smin5[0] = j; /* $|X|$, $h$ */
smax5[0] = j; /* $|X|$, $h$ */
smax5[1] = j+2*l; /* $\lambda$ */
smin5[2] = smax5[2] = (D2_FLAG < 0.e0) ? MIN(l, (DEG_E-1-j)/2) : 0; /* $D_2^*$ */
xiovera2powdmh[i][l][j] = fn_clone(xiovera2dmh[0][0], 4, smin5, smax5, NULL, NULL);
fn_set_sleep(xiovera2powdmh[i][l][j]);
}
}
for(j = 0; j < DEG_M; ++j) {
for(m = 0; m < DEG_M-j; ++m) {
prod_D2(xiovera2powdmh[j][l-1], xiovera2dmh[m], xiovera2powdmh[j+m][l],
(DOUBLE) (0.5e0-l)/(DOUBLE) l, TRUE);
}
}
for(j = 0; j < DEG_M; ++j) {
for(i = 0; i < DEG_E; ++i) {
fn_remove(xiovera2powdmh[j][l-1][i]);
fn_trim(xiovera2powdmh[j][l][i]);
}
}
dm_shuffle(0);
for(i = 0; i < DEG_E; ++i)
laplace(xiovera2powdmh, lapldmh, hamtmpdmh, l, i);
}
for(j = 0; j < DEG_M; ++j) {
for(i = 0; i < DEG_E; ++i) {
fn_remove(xiovera2dmh[j][i]);
fn_remove(xiovera2powdmh[j][DEG_E-1][i]);
}
}
for(i = 0; i < DEG_M; ++i) {
fn_remove(alpha2dmh[i]);
}
for(j = 0; j <= (3*(DEG_E-1)+MM_CUT)/2; ++j) {
for(i = 0; i < DEG_M; ++i) {
fn_remove(alphajdmh[j][i]);
fn_remove(lapldmh[j][i]);
}
}
dm_shuffle(0);
smin5[1] = 0; /* $\lambda$ */
smin5[2] = smax5[2] = (D2_FLAG < 0.e0) ? 1 : 0; /* $D_2^*$ */
for(i = 0; i < DEG_M; ++i) {
smin5[3] = i; /* $L$ */
smax5[3] = i; /* $L$ */
for(j = 0; j < DEG_E; ++j) {
smin5[0] = j; /* $|X|$, $h$ */
smax5[0] = j; /* $|X|$, $h$ */
smax5[1] = j+2; /* $\lambda$ */
kindmh[i][j] = fn_clone(hamtmpdmh[0][0], 4, smin5, smax5, NULL, NULL);
}
}
compute_kinetic(kindmh, wdmh);
for(i = 0; i < DEG_E; ++i) {
fn_remove(wdmh[i]);
}
add_kinetic(kindmh, hamtmpdmh);
for(i = 0; i < DEG_M; ++i) {
for(j = 0; j < DEG_E; ++j) {
fn_remove(kindmh[i][j]);
}
}
for(i = 0; i < DEG_M; ++i) {
for(j = 0; j < DEG_E; ++j) {
fn_trim(hamtmpdmh[i][j]);
}
}
for(i = 0; i < 8; ++i)
perm3[i] = i+1;
count3[0] = 2;
count3[1] = 2;
count3[2] = 4;
type3[0] = DER_POLYNOM;
type3[1] = DER_FOURIER;
type3[2] = DER_POLYNOM;
derfr = fn_def_fragm(8, perm3, 3, count3, type3);
/* */
type3[0] = PROD_POLYNOM;
type3[1] = PROD_FOURIER;
type3[2] = PROD_POLYNOM;
prodfr = fn_def_fragm(8, perm3, 3, count3, type3);
/* */
count3[0] = 8;
type3[0] = GEN_UNDEF;
genfr = fn_def_fragm(8, perm3, 1, count3, type3);
/* */
perm3[0] = 3;
perm3[1] = 4;
perm3[2] = 1;
perm3[3] = 2;
perm3[4] = 7;
perm3[5] = 8;
perm3[6] = 5;
perm3[7] = 6;
count3[0] = 4;
count3[1] = 4;
type3[0] = CAN_R_ACTION_ANGLE;
type3[1] = CAN_R_GENERIC;
canfr = fn_def_fragm(8, perm3, 2, count3, type3);
/* */
perm3[0] = 7;
perm3[1] = 8;
perm3[2] = 5;
perm3[3] = 6;
perm3[4] = 1;
perm3[5] = 3;
perm3[6] = 2;
perm3[7] = 4;
count3[0] = 4;
count3[1] = 2;
count3[2] = 2;
type3[0] = IDX_HOMOPOL; /* $xi$, $eta$ */
type3[1] = IDX_FOURIER; /* $\lambda$ */
type3[2] = IDX_HOMOPOL; /* $L$ */
smin3[0] = 0; /* $xi$, $eta$ */
smax3[0] = 0; /* $xi$, $eta$ */
smin3[1] = 0; /* $\lambda$ */
smax3[1] = 0; /* $\lambda$ */
for(i = 0; i <= DEG_M; ++i) {
smin3[2] = i; /* $L$ */
smax3[2] = i; /* $L$ */
idxfr = fn_def_idx_fragm(8, perm3, 3, count3, type3, smin3, smax3, NULL, NULL);
ham0dmh[i] = fn_create(8, 0, idxfr, derfr, prodfr, genfr, canfr, FN_BIN_TREE, CF_DOUBLE);
dm_release(idxfr);
}
init_h0(ham0dmh);
for(i = 0; i <= DEG_M; ++i) {
ham0dmh_1[i] = fn_duplicate(ham0dmh[i]);
ham0dmh_2[i] = fn_duplicate(ham0dmh[i]);
}
init_h0_split(ham0dmh_1, ham0dmh_2);
smin3[1] = 0; /* $\lambda$ */
smax3[1] = MM_CUT; /* $\lambda$ */
for(i = 0; i < DEG_M; ++i) {
smin3[2] = i; /* $L$ */
smax3[2] = i; /* $L$ */
for(j = 0; j < DEG_E; ++j) {
smin3[0] = j; /* $xi$, $eta$ */
smax3[0] = j; /* $xi$, $eta$ */
idxfr = fn_def_idx_fragm(8, perm3, 3, count3, type3, smin3, smax3, NULL, NULL);
ham1dmh[i][j] = fn_create(8, 0, idxfr, derfr, prodfr, genfr, canfr, FN_BIN_TREE, CF_DOUBLE);
dm_release(idxfr);
}
}
for(j = 0; j < DEG_E; ++j) {
secular2poinc(hamtmpdmh, j, ham1dmh);
for(i = 0; i < DEG_M; ++i)
fn_remove(hamtmpdmh[i][j]);
for(i = 0; i < DEG_M; ++i)
fn_trim(ham1dmh[i][j]);
}
dm_release(derfr);
dm_release(prodfr);
dm_release(genfr);
dm_release(canfr);
for(i = 0; i <= DEG_M; ++i) {
sprintf(filename, "ham0dmh_1_%d.bin", i);
fn_export_bin(ham0dmh_1[i], filename);
fn_remove(ham0dmh_1[i]);
sprintf(filename, "ham0dmh_2_%d.bin", i);
fn_export_bin(ham0dmh_2[i], filename);
fn_remove(ham0dmh_2[i]);
}
for(i = i; i < DEG_M; ++i) {
for(j = 0; j < DEG_E; ++j) {
sprintf(filename, "ham1dmh_%d_%d.bin", i, j);
fn_export_bin(ham1dmh[i][j], filename);
}
}
ofp1 = fopen("expansion.out", "w");
get_fn_handlers(ham0dmh[0], &fh);
for(i = 0; i <= DEG_M; ++i) {
j = fh.first_nz(&coef, k3, ham0dmh[i]);
while(j > 0) {
for(n = 0; n < 8; n++)
fprintf(ofp1, "%3d", k3[n]);
fprintf(ofp1, "%24.16le\n", coef);
j = fh.next_nz(&coef, k3, ham0dmh[i]);
}
}
get_fn_handlers(ham1dmh[0][0], &fh);
for(i = 0; i < DEG_M; ++i) {
for(j = 0; j < DEG_E; ++j) {
l = fh.first_nz(&coef, k3, ham1dmh[i][j]);
while(l > 0) {
for(n = 0; n < 8; n++)
fprintf(ofp1, "%3d", k3[n]);
fprintf(ofp1, "%24.16le\n", coef);
l = fh.next_nz(&coef, k3, ham1dmh[i][j]);
}
}
}
fclose(ofp1);
for(i = 0; i <= DEG_M; ++i)
fn_remove(ham0dmh[i]);
for(i = 0; i < DEG_M; ++i)
for(j = 0; j < DEG_E; ++j)
fn_remove(ham1dmh[i][j]);
dm_shuffle(0);
db_list_summary();
chr_finish();
printf("Test completed.\n");
return(0);
}
void one2two_frag ( DMHDR *  srcdmh,
DMHDR *  dstdmh 
)

This routine transforms an expansion in one trigonometric variable, into an expansion in three variables, using the second one.

Parameters:
*srcdmhsource function (input)
*dstdmhdestination function (output)

Technical details

The function srcdmh is expanded in $l$.

The function dstdmh is expanded in $e$ and $(l, h)$.

Definition at line 3242 of file threebody_plane.c.

{
int j;
DOUBLE coef;
INT2 k1[1];
INT2 k2[3];
FN_HANDLERS srch, dsth;
#if DEBUG
printf("Entry one2two_frag.\n");
#endif
get_fn_handlers(srcdmh, &srch);
get_fn_handlers(dstdmh, &dsth);
k2[0] = 0;
k2[2] = 0;
j = srch.first_nz(&coef, k1, srcdmh);
while(j > 0) {
k2[1] = k1[0];
dsth.put_coef(&coef, k2, dstdmh);
j = srch.next_nz(&coef, k1, srcdmh);
}
#if DEBUG
printf("Exit one2two_frag.\n");
#endif
return;
}
void prod_D2 ( DMHDR *  f1dmh[DEG_E],
DMHDR *  f2dmh[DEG_E],
DMHDR *  resdmh[DEG_E],
DOUBLE  coef,
UNS2  flag_D2 
)

This routine computes the product of the functions f1dmh[DEG_E] and f2dmh[DEG_E] times the coefficient coef.

Parameters:
*f1dmhsource function (input)
*f2dmhsource function (input)
*resdmhdestination function (output)
coefcoefficient
flag_D2truncation flag

Technical details

The functions fn1dmh[DEG_E], fn2dmh[DEG_E] and resdmh[DEG_E] are expanded in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

If the flag flag_D2 is set to TRUE the product is computed by calling the routine prod_D2_cut.

Cleans roundoff errors and removes neglegible terms in resdmh[DEG_E] by calling the routine cutepsrel.

Definition at line 3292 of file threebody_plane.c.

{
int l, j;
UNS2 f1stat, f2stat, rstat;
DMHDR *tmpdmh;
#if DEBUG
printf("Entry prod_D2.\n");
#endif
for(l = 0; l < DEG_E; ++l) {
rstat = fn_open_write(resdmh[l]);
for(j = 0; j <= l; ++j) {
f1stat = fn_open_read(f1dmh[l-j]);
f2stat = fn_open_read(f2dmh[j]);
if(flag_D2) {
prod_D2_cut(f1dmh[l-j], f2dmh[j], coef, resdmh[l]);
}
else {
tmpdmh = prod(f1dmh[l-j], f2dmh[j]);
add_scalarmultiply_to(tmpdmh, &coef, resdmh[l]);
fn_remove(tmpdmh);
}
if(f1stat != 0) fn_restore_status(f1dmh[l-j], f1stat);
if(f2stat != 0) fn_restore_status(f2dmh[j], f2stat);
}
cutepsrel(resdmh[l], 1.e-15);
if(rstat != 0) fn_restore_status(resdmh[l], rstat);
}
#if DEBUG
printf("Exit prod_D2.\n");
#endif
return;
}
void prod_D2_cut ( DMHDR *  fn1dmh,
DMHDR *  fn2dmh,
DOUBLE  coef,
DMHDR *  resdmh 
)

This routine computes the truncated product of the functions fn1dmh and fn2dmh times the coefficient coef. The limit in the truncation depends on the degree in the parameter $D_2^*$.

Parameters:
*fn1dmhsource function (input)
*fn2dmhsource function (input)
coefcoefficient
*resdmhdestination function (output)

Technical details

The functions fn1dmh, fn2dmh and resdmh are expanded in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

The products

\[ \left({\sin \atop \cos}(\alpha_1\lambda_1+\alpha_2\lambda_2) {\sin \atop \cos}(\beta_1 h_1+\beta_2 h_2)\right) \left({\sin \atop \cos}(\alpha'_1\lambda_1+\alpha'_2\lambda_2) {\sin \atop \cos}(\beta'_1 h_1+\beta'_2 h_2)\right) \]

are computed using the elementay trigonometric relations.

Computes the products only if

\[ 2 {\tt k({D_2^*})} + {\tt k({|X_1|})} + {\tt k({|X_2|})} < {\tt DEG\_E}\,. \]

Indeed, the parameter $D_2^*$ is of order 2 in the eccentricities.

Definition at line 3358 of file threebody_plane.c.

{
int i1, i2, j, l, m;
DOUBLE coef1, coef2, tmp1, tmp2, coefres;
INT2 k1_5[9], k2_5[9], ktmp[9];
INT2 kfoures[4][4], sign12, sign22;
INT2 flag11, flag12, flag1;
INT2 flag21, flag22, flag2;
INT2 flagtmp1, flagtmp2, flagtmp;
FN_HANDLERS fh1, fh2, resh;
#if DEBUG
printf("Entry prod_D2_cut.\n");
#endif
get_fn_handlers(fn1dmh, &fh1);
get_fn_handlers(fn2dmh, &fh2);
get_fn_handlers(resdmh, &resh);
coef /= 4.e0;
i1 = fh1.first_nz(&coef1, k1_5, fn1dmh);
while(i1 > 0) {
tmp1 = coef*coef1;
flag11 = flagsincos(k1_5+5, 2);
flag12 = flagsincos(k1_5+7, 2);
if(flag11 == 0) flag11 = 1;
if(flag12 == 0) flag12 = 1;
i2 = fh2.first_nz(&coef2, k2_5, fn2dmh);
while(i2 > 0) {
for(j = 0; j < 5; ++j)
ktmp[j] = k1_5[j]+k2_5[j];
if(2*ktmp[0]+ktmp[3]+ktmp[4] < DEG_E) {
flag21 = flagsincos(k2_5+5, 2);
flag22 = flagsincos(k2_5+7, 2);
if(flag21 == 0) flag21 = 1;
if(flag22 == 0) flag22 = 1;
flag1 = flag11*flag21;
flag2 = flag12*flag22;
tmp2 = tmp1*coef2;
for(j = 0; j < 2; ++j) {
sign12 = 1-2*j;
for(l = 0; l < 2; ++l) {
sign22 = 1-2*l;
for(m = 0; m < 2; ++m)
kfoures[j*2+l][m] = k1_5[5+m]+sign12*k2_5[5+m];
for(m = 2; m < 4; ++m)
kfoures[j*2+l][m] = k1_5[5+m]+sign22*k2_5[5+m];
}
}
for(j = 0; j < 2; ++j) {
for(l = 0; l < 2; ++l) {
flagtmp1 = flagsincos(kfoures[0]+(j*2+l)*4, 2);
if(flagtmp1 != 0 || flag1 == 1) {
flagtmp2 = flagsincos(kfoures[0]+(j*2+l)*4+2, 2);
if(flagtmp2 != 0 || flag2 == 1) {
coefres = tmp2;
if(flag21 == -1) coefres *= flag11*(1-2*j);
if(flagtmp1 == 0) flagtmp1 = 1;
flagtmp = flagtmp1*flag1;
if(flag1 == -1 && flagtmp != 1) coefres *= flagtmp;
for(m = 0; m < 2; ++m)
ktmp[5+m] = flagtmp*kfoures[(j*2+l)][m];
if(flag22 == -1) coefres *= flag12*(1-2*l);
if(flagtmp2 == 0) flagtmp2 = 1;
flagtmp = flagtmp2*flag2;
if(flag2 == -1 && flagtmp != 1) coefres *= flagtmp;
for(m = 2; m < 4; ++m)
ktmp[5+m] = flagtmp*kfoures[(j*2+l)][m];
resh.add_coef(&coefres, ktmp, resdmh);
}
}
}
}
}
i2 = fh2.next_nz(&coef2, k2_5, fn2dmh);
}
i1 = fh1.next_nz(&coef1, k1_5, fn1dmh);
}
#if DEBUG
printf("Exit prod_D2_cut.\n");
#endif
return;
}
void rovera ( DMHDR *  roveradmh[DEG_E],
DMHDR *  wdmh[DEG_E] 
)

This routine computes the expansion of $\frac{\|r\|}{a}$ in $e$ and $(\lambda, h)$.

Parameters:
*roveradmhexpansion of $\frac{\|r\|}{a}$ (output)
*wdmhexpansion of $W(f,e)$ (input)

Technical details

The function roveradmh[DEG_E] is expanded in $e$ and $(\lambda, h)$.

The function wdmh[DEG_E] is expanded in $f$. To keep track of the expansion in $e$ we use the index in the array of the function.

\begin{align*} \frac{\|r\|}{a} &= \sum_{i=0}^{{\tt DEG\_E}-1} {\tt triliedmh[0][i]}\,e^i\\ &= \sum_{i=0}^{{\tt DEG\_E}-1} {\tt roveradmh[i]}\,. \end{align*}

Creates the function triliedmh[DEG_E][DEG_E] for $\frac{\|r\|}{a}$, originally expanded in $f$. Later is expanded $l$. To keep track of the expansion in $e$ we use the index in the array of the function.

trilie[i][j] f / l
perm 1
der FOU
prod FOU
can UND
gen UND
idx FOU
smin 0
smax i+j

Computes the expansion of $\frac{\|r\|}{a}$ in $f$, by calling the routine init_roveraf and stores in the function triliedmh[0][DEG_E]

Computes the expansion of $\frac{\|r\|}{a}$ in $l$, by calling the routine tri_lie and stores in the function triliedmh[DEG_E][0].

Computes the expansion of $\frac{\|r\|}{a}$ in $e$ and $(\lambda, h)$, by using the relation $l=\lambda+h$ and stores in the function roveradmh[DEG_E].

Remarks:
The routine tri_lie requires an expansion containing the factorials in the denominators.

Definition at line 3485 of file threebody_plane.c.

{
int i, j;
DOUBLE coef, fact = 1.0e0;
INT2 smin1[1], smax1[1], k1[1];
INT2 k2[3];
DMHDR *triliedmh[DEG_E][DEG_E];
FN_HANDLERS rah, trih;
#if DEBUG
printf("Entry rovera.\n");
#endif
smin1[0] = 0; /* $f$ or $l$ */
for(i = 0; i < DEG_E; ++i) {
for(j = 0; j < DEG_E-i; ++j) {
smax1[0] = i+j; /* $f$ or $l$ */
triliedmh[i][j] = fn_clone(wdmh[0], 1, smin1, smax1, NULL, NULL);
}
}
init_roveraf(triliedmh[0]);
tri_lie(triliedmh, wdmh);
for(i = 0; i < DEG_E; ++i)
for(j = 1; j < DEG_E-i; ++j)
fn_remove(triliedmh[i][j]);
get_fn_handlers(triliedmh[0][0], &trih);
get_fn_handlers(roveradmh[0], &rah);
for(i = 0; i < DEG_E; ++i) {
k2[0] = i; /* $e$ */
j = trih.first_nz(&coef, k1, triliedmh[i][0]);
while(j > 0) {
k2[1] = k1[0]; /* $\lambda$ */
k2[2] = k1[0]; /* $h$ */
coef /= fact;
rah.put_coef(&coef, k2, roveradmh[i]);
j = trih.next_nz(&coef, k1, triliedmh[i][0]);
}
fn_remove(triliedmh[i][0]);
fact *= (DOUBLE) i+1.0e0;
}
#if DEBUG
printf("Exit rovera.\n");
#endif
return;
}
void secular2poinc ( DMHDR *  srcdmh[DEG_M][DEG_E],
int  ecc,
DMHDR *  dstdmh[DEG_M][DEG_E] 
)

This routine replaces the terms $ |X_1|^a |X_2|^b \ {\sin \atop \cos}(\alpha h_1 + \beta h_2)$ in srcdmh[DEG_M][ecc] with their expansion in $(\xi, \eta)$ and stores the result in the function dstdmh[DEG_M][ecc].

Parameters:
*srcdmhsource function (input)
eccdegree in the eccentricities
*dstdmhthe destination function (output)
Returns:
void

Technical details

The function srcdmh[DEG_M][DEG_E] is expanded in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

The function dstdmh[DEG_M][DEG_E] is expanded in $L$, $\lambda$ and $(\xi, \eta)$.

Creates the temporary function tmpvetdmh[DEG_M] for $ |X_1|^a |X_2|^b$, expanded in $L$.

tmpvetdmh[i] L
perm 1-2
der POL
prod POL
can UND
gen UND
idx HPOL
smin i
smax i

Computes the expansion of $|X_1|^{\tt kpol[0]} |X_2|^{\tt kpol[1]}$ in $L$ and $\sqrt{\Lambda}|X|$ by calling init_pol2poincsec and stores in the function tmpvetdmh[DEG_M]. The expansion is explicit only in $L$.

Replaces the terms $ (\sqrt{\Lambda_1}|X_1|)^a (\sqrt{\Lambda_2}|X_2|)^b \ {\sin \atop \cos}(\alpha h_1 + \beta h_2)$ in srcdmh[DEG_M][ecc] with their expansions in $\xi$ and $\eta$ by calling the routine trig2poinc and stores in the function dstdmh[DEG_M][DEG_E].

Definition at line 3599 of file threebody_plane.c.

{
int i, j, l, m;
DOUBLE coef, tmp;
UNS2 perm1[2], count1[1], type1[1];
INT2 smin1[1], smax1[1], k1[2];
INT2 k5[9], k5new[9], kpol[2], ktrg[2];
DMHDR *idxfr, *derfr, *prodfr, *genfr, *canfr;
DMHDR *tmpvetdmh[DEG_M], *tmpdmh;
UNS2 fnstat, fnstatvet[DEG_M];
FN_HANDLERS srch, fh;
extern DOUBLE D2_fix;
#if DEBUG
printf("Entry secular2poinc.\n");
#endif
for(i = 0; i < 2; ++i)
perm1[i] = i+1;
count1[0] = 2;
type1[0] = FN_DER_POLYNOM;
derfr = fn_def_fragm(2, perm1, 1, count1, type1);
/* */
type1[0] = FN_PROD_POLYNOM;
prodfr = fn_def_fragm(2, perm1, 1, count1, type1);
/* */
type1[0] = FN_GEN_UNDEF;
genfr = fn_def_fragm(2, perm1, 1, count1, type1);
/* */
type1[0] = FN_CAN_UNDEF;
canfr = fn_def_fragm(2, perm1, 1, count1, type1);
type1[0] = IDX_HOMOPOL; /* $L$ */
for(i = 0; i < DEG_M; ++i) {
smin1[0] = i; /* $L$ */
smax1[0] = i; /* $L$ */
idxfr = fn_def_idx_fragm(2, perm1, 1, count1, type1, smin1, smax1, NULL, NULL);
tmpvetdmh[i] = fn_create(2, 0, idxfr, derfr, prodfr, genfr, canfr, FN_BIN_TREE, CF_DOUBLE);
dm_release(idxfr);
}
dm_release(derfr);
dm_release(prodfr);
dm_release(genfr);
dm_release(canfr);
get_fn_handlers(srcdmh[0][0], &srch);
get_fn_handlers(tmpvetdmh[0], &fh);
for(i = 0; i < DEG_M; ++i)
fnstatvet[i] = fn_open_write(dstdmh[i][ecc]);
for(i = 0; i < DEG_M; ++i) {
fnstat = fn_open_read(srcdmh[i][ecc]);
kpol[0] = -1;
kpol[1] = -1;
j = srch.first_nz(&coef, k5, srcdmh[i][ecc]);
while(j > 0) {
if(k5[3] != kpol[0] || k5[4] != kpol[1]) {
kpol[0] = k5[3];
kpol[1] = k5[4];
for(l = 0; l < DEG_M; ++l) {
tmpdmh = fn_duplicate(tmpvetdmh[l]);
fn_remove(tmpvetdmh[l]);
tmpvetdmh[l] = fn_duplicate(tmpdmh);
fn_remove(tmpdmh);
}
init_pol2poincsec(kpol, tmpvetdmh);
}
for(l = 0; l < 9; ++l)
k5new[l] = k5[l];
for(l = 0; l < DEG_M-i; ++l) {
m = fh.first_nz(&tmp, k1, tmpvetdmh[l]);
while(m > 0) {
tmp *= coef*pow(D2_fix, k5[0]);
k5new[1] = k5[1]+k1[0];
k5new[2] = k5[2]+k1[1];
if(k5new[7] == 0 || k5new[8] == 0) {
trig2poinc(k5new, k5new+3, k5new+7, tmp, dstdmh[i+l][ecc]);
}
else {
ktrg[0] = k5new[7];
ktrg[1] = ABS(k5new[8]);
trig2poinc(k5new, k5new+3, ktrg, tmp, dstdmh[i+l][ecc]);
ktrg[0] = -k5new[7];
ktrg[1] = -ktrg[1];
if(k5new[8] > 0) tmp = -tmp;
trig2poinc(k5new, k5new+3, ktrg, tmp, dstdmh[i+l][ecc]);
}
m = fh.next_nz(&tmp, k1, tmpvetdmh[l]);
}
}
j = srch.next_nz(&coef, k5, srcdmh[i][ecc]);
}
if(fnstat != 0) fn_restore_status(srcdmh[i][ecc], fnstat);
}
for(i = 0; i < DEG_M; ++i)
fn_remove(tmpvetdmh[i]);
for(i = 0; i < DEG_M; ++i)
if(fnstatvet[i] != 0) fn_restore_status(dstdmh[i][ecc], fnstatvet[i]);
#if DEBUG
printf("Exit secular2poinc.\n");
#endif
return;
}
void sincosfh ( DMHDR *  sinfhdmh[DEG_E],
DMHDR *  cosfhdmh[DEG_E],
DMHDR *  wdmh[DEG_E] 
)

This routine computes the expansions of $\sin(f-h)$ and $\cos(f-h)$ in $e$ and $(\lambda, h)$.

Parameters:
*sinfhdmhexpansion of $\sin(f-h)$ (output)
*cosfhdmhexpansion of $\cos(f-h)$ (output)
*wdmhexpansion of $W(f,e)$ (input)

Technical details

\begin{align*} \sin(f-h) &= \sin(f)\cos(h) - \cos(f)\sin(h)\\ &= \sum_{i=0}^{{\tt DEG\_E}-1} {\tt sinfhdmh[i]}\,,\\ \cos(f-h) &= \cos(f)\cos(h) + \sin(f)\sin(h)\\ &= \sum_{i=0}^{{\tt DEG\_E}-1} {\tt cosfhdmh[i]}\,. \end{align*}

The functions sinfhdmh[DEG_E] and cosfhdmh[DEG_E] are expanded in $e$ and $(\lambda, h)$.

Creates the functions sinfdmh[DEG_E][DEG_E] and cosfdmh[DEG_E][DEG_E] for $\sin(f)$ and $\cos(f)$, originally expanded in $f$. Later are expanded $l$. To keep track of the expansion in $e$ we use the index in the array of the function.

(sin/cos)fdmh[i][j] f / l
perm 1
der FOU
prod FOU
can UND
gen UND
idx FOU
smin 0
smax i+j+1
Remarks:
The terms $\sin(f)$ and $\cos(f)$ have degree 0 in $e$ and degree 1 in $f$. For this reason the terms of degree i in $e$ have degree i+1 in $f$.

Creates the functions sinhdmh and coshdmh for the terms $\sin(h)$ and $\cos(h)$, expanded in $e$ and $(\lambda, h)$.

(sin/cos)hdmh e lam & h
perm 1 2-3
der POL FOU
prod POL FOU
can UND UND
gen UND UND
idx HPOL FOU
smin 0 0
smax 0 1

Computes the expansions of $\sin(f)$ and $\cos(f)$ in $f$, and stores in sinfdmh[0][0] and cosfdmh[0][0].

Computes the expansions of $\sin(f)$ and $\cos(f)$ in $l$ by calling the routine tri_lie and stores in sinfdmh[DEG_E][0] and cosfdmh[DEG_E][0].

Computes the expansions of $\sin(h)$ and $\cos(h)$ in $e$ and $(l, h)$, and stores in sinhdmh and coshdmh.

Computes the expansions of $\sin(f)$ and $\cos(f)$ in $e$ and $(l, h)$, by calling the routine one2two_frag.

Computes the expansions of $\sin(f-h)$ and $\cos(f-h)$ in $e$ and $(\lambda, h)$, by using the relation $l=\lambda+h$, and stores in sinfhdmh[DEG_E] and cosfhdmh[DEG_E].

Remarks:
The routine tri_lie computes an expansion containing the factorials in the denominators.

Definition at line 3757 of file threebody_plane.c.

{
int i, j;
DOUBLE coef, fact = 1.e0;
INT2 smin1[1], smax1[1];
INT2 k1[1];
INT2 smin2[3], smax2[3];
INT2 k2[3];
DMHDR *sinfdmh[DEG_E][DEG_E];
DMHDR *cosfdmh[DEG_E][DEG_E];
DMHDR *sinhdmh, *coshdmh;
DMHDR *tmpsindmh, *tmpcosdmh;
DMHDR *tmpdmh;
FN_HANDLERS fh;
#if DEBUG
printf("Entry sincosfh.\n");
#endif
smin1[0] = 0; /* $f$ or $l$ */
for(i = 0; i < DEG_E; ++i) {
for(j = 0; j < DEG_E; ++j) {
smax1[0] = i+j+1; /* $f$ or $l$ */
sinfdmh[i][j] = fn_clone(wdmh[0], 1, smin1, smax1, NULL, NULL);
cosfdmh[i][j] = fn_duplicate(sinfdmh[i][j]);
}
}
smin2[0] = 0; /* $e$ */
smax2[0] = 0; /* $e$ */
smin2[1] = 0; /* $\lambda$, $h$ */
smax2[1] = 1; /* $\lambda$, $h$ */
sinhdmh = fn_clone(sinfhdmh[0], 2, smin2, smax2, NULL, NULL);
coshdmh = fn_duplicate(sinhdmh);
get_fn_handlers(sinfdmh[0][0], &fh);
/* $-sin(-f)$ */
k1[0] = -1;
coef = -1.e0;
fh.put_coef(&coef, k1, sinfdmh[0][0]);
/* $cos(f)$ */
k1[0] = 1;
coef = 1.e0;
fh.put_coef(&coef, k1, cosfdmh[0][0]);
tri_lie(sinfdmh, wdmh);
tri_lie(cosfdmh, wdmh);
for(i = 0; i < DEG_E; ++i) {
for(j = 1; j < DEG_E; ++j) {
fn_remove(sinfdmh[i][j]);
fn_remove(cosfdmh[i][j]);
}
}
get_fn_handlers(sinhdmh, &fh);
/* $-sin(-h)$ */
k2[0] = 0;
k2[1] = 0;
k2[2] = -1;
coef = -1.e0;
fh.put_coef(&coef, k2, sinhdmh);
/* $cos(h)$ */
k2[0] = 0;
k2[1] = 0;
k2[2] = 1;
coef = 1.e0;
fh.put_coef(&coef, k2, coshdmh);
get_fn_handlers(sinfhdmh[0], &fh);
smin2[0] = 0;
smax2[0] = 0;
smin2[1] = 0;
for(i = 0; i < DEG_E; ++i) {
smax2[1] = 2*i+3; /* $\lambda$, $h$ */
tmpsindmh = fn_clone(sinfhdmh[0], 2, smin2, smax2, NULL, NULL);
tmpcosdmh = fn_duplicate(tmpsindmh);
smax2[1] = i+1; /* $l$, $h$ */
tmpdmh = fn_clone(sinfhdmh[0], 2, smin2, smax2, NULL, NULL);
one2two_frag(sinfdmh[i][0], tmpdmh);
fn_remove(sinfdmh[i][0]);
add_prod_to(tmpdmh, coshdmh, tmpsindmh);
add_prod_to(tmpdmh, sinhdmh, tmpcosdmh);
fn_remove(tmpdmh);
smax2[1] = i+1; /* $l$, $h$ */
tmpdmh = fn_clone(cosfhdmh[0], 2, smin2, smax2, NULL, NULL);
one2two_frag(cosfdmh[i][0], tmpdmh);
fn_remove(cosfdmh[i][0]);
add_prod_to(tmpdmh, coshdmh, tmpcosdmh);
coef = -1.e0;
scalar_multiply(tmpdmh, &coef);
add_prod_to(tmpdmh, sinhdmh, tmpsindmh);
fn_remove(tmpdmh);
j = fh.first_nz(&coef, k2, tmpsindmh);
while(j > 0) {
k2[0] = i; /* $e$ */
k2[2] += k2[1]; /* $h$ */
coef /= fact;
fh.put_coef(&coef, k2, sinfhdmh[i]);
j = fh.next_nz(&coef, k2, tmpsindmh);
}
fn_remove(tmpsindmh);
j = fh.first_nz(&coef, k2, tmpcosdmh);
while(j > 0) {
k2[0] = i; /* $e$ */
k2[2] += k2[1]; /* $h$ */
coef /= fact;
fh.put_coef(&coef, k2, cosfhdmh[i]);
j = fh.next_nz(&coef, k2, tmpcosdmh);
}
fn_remove(tmpcosdmh);
fact *= (DOUBLE) i+1.0e0;
}
fn_remove(sinhdmh);
fn_remove(coshdmh);
#if DEBUG
printf("Exit sincosfh.\n");
#endif
return;
}
void tri_lie ( DMHDR *  triliedmh[DEG_E][DEG_E],
DMHDR *  wdmh[DEG_E] 
)

This routine transforms the expansion of $F(f,e)$ in $e$ and $f$, into $G(l,e)$ in $e$ and $l$, by using the generator $W(f,e)$ of the Lie serie. The expansion is explicit only in $f$ or $l$.

Parameters:
*triliedmhexpansion of $F(f,e)$ to be transformed in $G(l,e)$ (input/output)
*wdmhthe generator $W(f,e)$ (input)
Returns:
void

Technical details

The function triliedmh[DEG_E][DEG_E] is originally expanded in $e$ and $f$. Later is expanded in $e$ and $l$. To keep track of the expansion in $e$ we use the index in the array of the function.

The function wdmh[DEG_E] is expanded in $f$. To keep track of the expansion in $e$ we use the index in the array of the function.

\begin{align*} F(f,e) &= \sum_{i\geq0} \frac{F_i^0}{i!}e^i = \sum_{i\geq0} \frac{{\tt triliedmh[0][i]}}{i!}e^i\,.\\ W(f,e) &= \sum_{i\geq0} \frac{W_i}{i!}e^i = \sum_{i\geq0} \frac{{\tt wdmh[i]}}{i!}e^i\,.\\ G(l,e) &= F(f(l,e),e) = \sum_{i\geq0} \frac{F_0^i}{i!}e^i = \sum_{i\geq0} \frac{{\tt triliedmh[i][0]}}{i!}e^i\,. \end{align*}

The transformed function is computed by the recursive formula

\[ F_{i-j}^{j} = F_{i-j+1}^{j-1} + \sum_{k=0}^{i-j} {i-j \choose k} \langle\nabla F_{i-j-k}^{j-1},W_{k}\rangle\,. \]

Proceeding line-by-line in the Lie triangle, we have to compute only one new derivative at each step, $\nabla F_{i-j}^{j-1}$.

Definition at line 3973 of file threebody_plane.c.

{
int i, j, k;
DOUBLE comb[DEG_E][DEG_E];
DMHDR *derdmh[DEG_E][DEG_E];
DMHDR *tmpdmh;
#if DEBUG
printf("Entry tri_lie.\n");
#endif
for(i = 1; i < DEG_E; ++i) {
comb[i-1][0] = 1.e0;
for(k = 1; k <= i-2; ++k) {
comb[i-1][k] = comb[i-2][k]+comb[i-2][k-1];
}
comb[i-1][i-1] = 1.e0;
for(j = 1; j <= i; ++j) {
derdmh[j-1][i-j] = deriv(triliedmh[j-1][i-j], 0);
add_to(triliedmh[j-1][i-j+1], triliedmh[j][i-j]);
for(k = 0; k <= i-j; ++k) {
tmpdmh = prod(wdmh[k], derdmh[j-1][i-j-k]);
scalar_multiply(tmpdmh, &comb[i-j][k]);
add_to(tmpdmh, triliedmh[j][i-j]);
fn_remove(tmpdmh);
}
}
}
for(i = 1; i < DEG_E; ++i)
for(j = 1; j <= i; ++j)
fn_remove(derdmh[j-1][i-j]);
#if DEBUG
printf("Exit tri_lie.\n");
#endif
return;
}
void trig2poinc ( INT2  k5[9],
INT2  kpol[2],
INT2  ktrg[2],
DOUBLE  coef,
DMHDR *  fndmh 
)

This routine transforms an expansion in $\sqrt{\Lambda}|X|$ and $h$, into an expansion in $(\xi, \eta)$.

Parameters:
k5full array of exponents
kpolarray of polynomial exponents
ktrgarray of trigonometric exponents
coefcoefficient
*fndmhexpansion of the function to be transformed (input/output)

Technical details

The function fndmh is expanded in $L$, $\lambda$ and $(\xi, \eta)$.

\begin{align*} (\sqrt{\Lambda_1}|X_1|)^{{\tt kpol[0]}}\ &{\sin \atop \cos}\Bigl({\tt ktrg[0]} h_1\Bigr)\ (\sqrt{\Lambda_2}|X_2|)^{{\tt kpol[1]}}\ {\sin \atop \cos}\Bigl({\tt ktrg[1]} h_2\Bigr)=\\ &(\sqrt{\Lambda_1}|X_1|)^{{\tt kpol[0]}-|{\tt ktrg[0]}|}\ (\sqrt{\Lambda_2}|X_2|)^{{\tt kpol[1]}-|{\tt ktrg[1]}|} (\sqrt{\Lambda_1}|X_1|)^{|{\tt ktrg[0]}|}\ {\sin \atop \cos}\Bigl({\tt ktrg[0]} h_1\Bigr) (\sqrt{\Lambda_2}|X_2|)^{|{\tt ktrg[1]}|}\ {\sin \atop \cos}\Bigl({\tt ktrg[1]} h_2\Bigr) \end{align*}

Remark that

\begin{align*} (\sqrt{\Lambda}|X|)^{k} &= (\sqrt{\Lambda}|X|)^{k} \left( \sin^{2}(h)+\cos^{2}(h)\right)^{k/2}\\ &= (\sqrt{\Lambda}|X|)^{k} \sum_{i=0}^{k/2} {k/2 \choose i} \cos^{k-2i}(h) \sin^{2i}(h)\\ &= \sum_{i=0}^{k/2} {k/2 \choose i} \xi^{k-2i} \eta^{2i} \end{align*}

Definition at line 4056 of file threebody_plane.c.

{
int i, j, i1, i2, i3, i4;
DOUBLE tmp, sign, newcoef;
DOUBLE coef1[2][DEG_E], coef2[2][DEG_E];
INT2 k3[8], length1[2], length2[2];
FN_HANDLERS fh;
#if DEBUG
printf("Entry trig2poinc.\n");
#endif
get_fn_handlers(fndmh, &fh);
for(i = 0; i < 2; ++i) {
for(j = 0; j < DEG_E; ++j) {
coef1[i][j] = coef2[i][j] = 0.0e0;
}
}
for(i = 0; i < 2; ++i) {
length1[i] = (kpol[i]-ABS(ktrg[i]))/2;
coef1[i][0] = 1.0e0;
for(j = 1; j <= length1[i]; ++j) {
coef1[i][j] = coef1[i][j-1]*(length1[i]-j+1.0e0)/(DOUBLE) j;
}
}
for(i = 0; i < 2; ++i) {
length2[i] = ABS(ktrg[i]);
tmp = 1.0e0;
sign = 1.0e0;
if(ktrg[i] < 0) {
for(j = 1; j <= length2[i]; ++j) {
tmp *= (length2[i]-j+1.0e0)/(DOUBLE) j;
if(j%2 == 1) {
sign = -sign;
coef2[i][j] = tmp*sign;
}
}
}
else {
coef2[i][0] = 1.0e0;
for(j = 1; j <= length2[i]; ++j) {
tmp *= (length2[i]-j+1.0e0)/(DOUBLE) j;
if(j%2 == 0) {
sign = -sign;
coef2[i][j] = tmp*sign;
}
}
}
}
k3[0] = k5[1];
k3[1] = k5[2];
k3[2] = k5[5];
k3[3] = k5[6];
for(i1 = 0; i1 <= length1[0]; ++i1) {
for(i2 = 0; i2 <= length1[1]; ++i2) {
for(i3 = 0; i3 <= length2[0]; ++i3) {
k3[4] = 2*(length1[0]-i1)+length2[0]-i3;
k3[6] = 2*i1+i3;
for(i4 = 0; i4 <= length2[1]; ++i4) {
k3[5] = 2*(length1[1]-i2)+length2[1]-i4;
k3[7] = 2*i2+i4;
newcoef = coef*coef1[0][i1]*coef1[1][i2]*coef2[0][i3]*coef2[1][i4];
fh.add_coef(&newcoef, k3, fndmh);
}
}
}
}
#if DEBUG
printf("Exit trig2poinc.\n");
#endif
return;
}
void two2five ( DMHDR *  srcdmh,
DMHDR *  dstdmh 
)

This routine transforms the expansion of the function srcdmh in $e$ and $(\lambda, h)$, into an expansion in $D_2^*$, $L$, $|X|$, $\lambda$, $h$, stored in dstdmh.

Parameters:
*srcdmhsource function (input)
*dstdmhdestination function (output)

Technical details

The function srcdmh is expanded in $e$ and $(\lambda, h)$.

The function dstdmh is expanded in $D_2^*$, $L$, $|X|$, $\lambda$, $h$.

Definition at line 4168 of file threebody_plane.c.

{
int i, j;
DOUBLE coef;
INT2 k2[3], k5[9];
FN_HANDLERS sh, dh;
#if DEBUG
printf("Entry two2five.\n");
#endif
get_fn_handlers(srcdmh, &sh);
get_fn_handlers(dstdmh, &dh);
j = sh.first_nz(&coef, k2, srcdmh);
while(j > 0) {
for(i = 0; i < 9; ++i)
k5[i] = 0;
k5[3] = k2[0];
if(k2[1] == 0 || k2[2] == 0) {
k5[5] = k2[1];
k5[7] = k2[2];
dh.put_coef(&coef, k5, dstdmh);
}
else {
k5[5] = k2[1];
k5[7] = ABS(k2[2]);
dh.put_coef(&coef, k5, dstdmh);
k5[5] = -k5[5];
k5[7] = -k5[7];
if(k2[2] > 0) coef = -coef;
dh.put_coef(&coef, k5, dstdmh);
}
j = sh.next_nz(&coef, k2, srcdmh);
}
#if DEBUG
printf("Exit two2five.\n");
#endif
return;
}
void two2five_a ( DMHDR *  srcdmh,
DMHDR *  dstdmh 
)

This routine transforms the expansion of the function srcdmh in $e$ and $(\lambda, h)$, into an expansion in $D_2^*$, $a$, $|X|$, $\lambda$, $h$, stored in dstdmh. Moreover the routine performs the multiplication by $a$.

Parameters:
*srcdmhsource function (input)
*dstdmhdestination function (output)

Technical details

The function srcdmh is expanded in $e$ and $(\lambda, h)$.

The function dstdmh is expanded in $D_2^*$, $a$, $|X|$, $\lambda$, $h$.

Definition at line 4231 of file threebody_plane.c.

{
int j;
DOUBLE coef;
INT2 k2[3];
INT2 k5[9];
FN_HANDLERS sh, dh;
#if DEBUG
printf("Entry two2five_a.\n");
#endif
get_fn_handlers(srcdmh, &sh);
get_fn_handlers(dstdmh, &dh);
for(j = 0; j < 9; ++j)
k5[j] = 0;
k5[1] = 1; /* $a$ */
j = sh.first_nz(&coef, k2, srcdmh);
while(j > 0) {
k5[3] = k2[0];
if(k2[1] == 0 || k2[2] == 0) {
k5[5] = k2[1];
k5[7] = k2[2];
dh.put_coef(&coef, k5, dstdmh);
}
else {
k5[5] = k2[1];
k5[7] = ABS(k2[2]);
dh.put_coef(&coef, k5, dstdmh);
k5[5] = -k5[5];
k5[7] = -k5[7];
if(k2[2] > 0) coef = -coef;
dh.put_coef(&coef, k5, dstdmh);
}
j = sh.next_nz(&coef, k2, srcdmh);
}
#if DEBUG
printf("Exit two2five_a.\n");
#endif
return;
}
void xi_over_a2 ( DMHDR *  xiovera2dmh[DEG_M][DEG_E],
DMHDR *  wdmh[DEG_E] 
)

This routine computes the expansion of $\frac{\Xi}{{a_2^*}^2}$ in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

Parameters:
*xiovera2dmhfunction $\frac{\Xi}{{a_2^*}^2}$ (output)
*wdmhgenerator $W(f,e)$ (input)

Technical details

The function xiovera2dmh[DEG_M][DEG_E] is expanded in $D_2^*$, $L$, $|X|$, $\lambda$ and $h$.

The function wdmh[DEG_E] is expanded in $f$. To keep track of the expansion in $e$ we use the index in the array of the function.

\[ \frac{\Xi}{{a_2^*}^2} = \sum_{i=0}^{{\tt DEG\_M}-1}\sum_{j=0}^{{\tt DEG\_E}-1} {\tt xiovera2dmh[i][j]}\,. \]

The term $\frac{\|r_1-r_2\|^2}{{a_2^*}^2}$ is computed as

\begin{align*} \frac{\|r_1-r_2\|^2}{{a_2^*}^2} = &\frac{\|r_1\|^2}{{a_2^*}^2}+\frac{\|r_2\|^2}{{a_2^*}^2}\\ &-2 \frac{\|r_1\|}{{a_2^*}} \cos(f_1-h_1)\ \frac{\|r_2\|}{{a_2^*}} \cos(f_2-h_2)\\ &-2 \frac{\|r_1\|}{{a_2^*}} \sin(f_1-h_1)\ \frac{\|r_2\|}{{a_2^*}} \sin(f_2-h_2)\,. \end{align*}

The term $\frac{\|r\|}{a_2^*}$ is computed as

\[ \frac{\|r\|}{a_2^*} = \frac{\|r\|}{a}\frac{a}{a_2^*}\,. \]

Creates the function tmproveradmh[DEG_E] for the term $\frac{\|r\|}{a}$, expanded in $e$ and $(\lambda, h)$.

tmproveradm[i] e lam & h
perm 1 2-3
der POL FOU
prod POL FOU
can UND UND
gen UND UND
idx HPOL FOU
smin i 0
smax i 2i

Computes the expansion of $\frac{\|r\|}{a}$ by calling the routine rovera and stores in the function tmproveradmh[DEG_E].

Creates the functions sinfhdmh[DEG_E] and cosfhdmh[DEG_E] for the terms $\sin(f-h)$ and $\cos(f-h)$, expanded in $e$ and $(\lambda, h)$.

(sin/cos)fhdmh[i] e lam & h
perm 1 2-3
der POL FOU
prod POL FOU
can UND UND
gen UND UND
idx HPOL FOU
smin i 0
smax i 2i+3

Computes the expansions of $\sin(f-h)$ and $\cos(f-h)$ by calling the routine sincosfh and stores in the functions sinfhdmh[DEG_E] and cosfhdmh[DEG_E].

Creates the function e2Xdmh[DEG_E][DEG_E] for the powers of the eccentricity, $e^i$, expanded in $|X|$ and $(\lambda, h)$.

e2Xdmh[i][j] X lam & h
perm 1 2-3
der POL FOU
prod POL FOU
can UND UND
gen UND UND
idx HPOL FOU
smin j 0
smax j 0

Computes the expansions of the powers of the eccentricity, $e^i$, by calling the routine init_e2X and stores in the function e2Xdmh[i][DEG_E].

Creates the function roveradmh[DEG_E], for the term $\frac{\|r\|}{a}$, expanded in $D_2^*$, $a$, $|X|$, $\lambda$ and $h$.

roveradmh[i] D2* a X lam h
perm 7 8-9 1-2 5-6 3-4
der POL POL POL FOU FOU
prod POL POL POL FOU FOU
can UND UND UND UND UND
gen UND UND UND UND UND
idx HPOL HPOL DAL1 FOU DAL2
smin 0 1 i 0 0
smax 0 1 i i i

Creates the functions rasin1dmh[DEG_E], racos1dmh[DEG_E], rasin2dmh[DEG_E] and racos2dmh[DEG_E] for the terms $\frac{\|r_1\|}{a_1} \sin(f_1-h_1)$, $\frac{\|r_1\|}{a_1} \cos(f_1-h_1)$, $\frac{\|r_2\|}{a_2} \sin(f_2-h_2)$ and $\frac{\|r_2\|}{a_2} \cos(f_2-h_2)$, expranded in $D_2^*$, $a$, $|X|$, $\lambda$ and $h$.

ra(sin/cos)(1/2)dmh[i] D2* a X lam h
perm 7 8-9 1-2 5-6 3-4
der POL POL POL FOU FOU
prod POL POL POL FOU FOU
can UND UND UND UND UND
gen UND UND UND UND UND
idx HPOL HPOL DAL 1 FOU DAL2
smin 0 1 i 0 0
smax 0 1 i i+1 i

Computes the expansion of $\frac{\|r\|}{a} \sin(f-h)$ and stores in tmpvetdmh[DEG_E].

Computes the expansion of $\frac{\|r\|}{a} \sin(f-h)$ as a function of $|X|$ by calling the routine e2X.

Computes the expansion of $a_1 \frac{\|r_1\|}{a_1} \sin(f_1-h_1)$ by calling the routine two2five_a and stores in rasin1dmh[DEG_E].

Computes the expansion of $\frac{\|r\|}{a} \cos(f-h)$ and stores in tmpvetdmh[DEG_E].

Computes the expansion of $\frac{\|r\|}{a} \cos(f-h)$ as a function of $|X|$ by calling the routine e2X.

Computes the expansion of $a_1\frac{\|r_1\|}{a_1} \cos(f_1-h_1)$ by calling the routine two2five_a and stores in racos1dmh[DEG_E].

Computes the expansions of $a_2 \frac{\|r_2\|}{a_2} \sin(f_2-h_2)$ and $a_2 \frac{\|r_2\|}{a_2} \cos(f_2-h_2)$ by calling the routine inner2outer and stores in rasin2dmh[DEG_E] and racos2dmh[DEG_E].

Computes the expansion of $\frac{\|r_1\|}{a_1}$ as a function of $|X|$ by calling the routine e2X.

Computes the expansion of $a_1 \frac{\|r_1\|}{a_1}$ by calling the routine two2five_a and stores in roveradmh[DEG_E].

Creates the function xidmh[DEG_E] for the term $\Xi$, expanded in $D_2^*$, $a$, $|X|$, $\lambda$ and $h$.

xidmh[i] D2* a X lam h
perm 7 8-9 1-2 5-6 3-4
der POL POL POL FOU FOU
prod POL POL POL FOU FOU
can UND UND UND UND UND
gen UND UND UND UND UND
idx HPOL HPOL DAL1 FOU DAL2
smin 0 2 i 0 0
smax 0 2 i i+2 i

Computes the expansions of $\left(a_1 \frac{\|r_1\|}{a_1} \right)^2$ and $\left(a_2 \frac{\|r_2\|}{a_2}\right)^2$ and stores in xidmh[DEG_E].

Computes the expansions of $-2 a_1 \frac{\|r_1\|}{a_1} \sin(f_1-h_1)\ a_2 \frac{\|r_2\|}{a_2} \sin(f_2-h_2)$ and $-2 a_1 \frac{\|r_1\|}{a_1} \cos(f_1-h_1)\ a_2 \frac{\|r_2\|}{a_2} \cos(f_2-h_2)$ and stores in xidmh[DEG_E].

Replaces the terms $a_1^x a_2^y$ by their expansions in $L$, and perform the division by ${a_2^*}^2$, by calling the routine a2L and stores in xiovera2dmh[DEG_E].

Subtracts the circular approximation, $\frac{a_1^2 + a_2^2 - 2 a_1 a_2 \cos(\lambda_1-\lambda_2)}{{a_2^*}^2}$ and stores in xiovera2dmh[DEG_E].

Cleans the roundoff errors and the neglegible terms in xiovera2dmh[DEG_E] by calling the routine cutepsrel.

Definition at line 4294 of file threebody_plane.c.

{
int i, j;
DOUBLE coef;
UNS2 perm2[3], count2[2], type2[2];
INT2 smin2[2], smax2[2];
INT2 smin5[5], smax5[5];
INT2 k5[9];
DMHDR *idxfr, *derfr, *prodfr, *genfr, *canfr;
DMHDR *xidmh[DEG_E];
DMHDR *roveradmh[DEG_E];
DMHDR *sinfhdmh[DEG_E], *cosfhdmh[DEG_E];
DMHDR *e2Xdmh[DEG_E][DEG_E];
DMHDR *rasin1dmh[DEG_E], *racos1dmh[DEG_E];
DMHDR *rasin2dmh[DEG_E], *racos2dmh[DEG_E];
DMHDR *tmp0dmh, *tmp1dmh;
DMHDR *tmproveradmh[DEG_E];
DMHDR *tmpvetdmh[DEG_E];
FN_HANDLERS fh;
extern DOUBLE Alpha_fix, D2_fix;
#if DEBUG
printf("Entry xi_over_a2.\n");
#endif
for(i = 0; i < 3; ++i)
perm2[i] = i+1;
count2[0] = 1;
count2[1] = 2;
type2[0] = DER_POLYNOM;
type2[1] = DER_FOURIER;
derfr = fn_def_fragm(3, perm2, 2, count2, type2);
/* */
type2[0] = PROD_POLYNOM;
type2[1] = PROD_FOURIER;
prodfr = fn_def_fragm(3, perm2, 2, count2, type2);
/* */
count2[0] = 3;
type2[0] = GEN_UNDEF;
genfr = fn_def_fragm(3, perm2, 1, count2, type2);
/* */
type2[0] = CAN_UNDEF;
canfr = fn_def_fragm(3, perm2, 1, count2, type2);
/* */
count2[0] = 1;
count2[1] = 2;
type2[0] = IDX_HOMOPOL; /* $e$ */
type2[1] = IDX_FOURIER; /* $\lambda$, $h$ */
smin2[1] = 0; /* $\lambda$, $h$ */
for(i = 0; i < DEG_E; ++i) {
smin2[0] = i; /* $e$ */
smax2[0] = i; /* $e$ */
smax2[1] = 2*i; /* $\lambda$, $h$ */
idxfr = fn_def_idx_fragm(3, perm2, 2, count2, type2, smin2, smax2, NULL, NULL);
tmproveradmh[i] = fn_create(3, 0, idxfr, derfr, prodfr, genfr, canfr, FN_BIN_TREE, CF_DOUBLE);
dm_release(idxfr);
}
dm_release(derfr);
dm_release(prodfr);
dm_release(genfr);
dm_release(canfr);
rovera(tmproveradmh, wdmh);
smin2[1] = 0; /* $\lambda$, $h$ */
for(i = 0; i < DEG_E; ++i) {
smin2[0] = i; /* $e$ */
smax2[0] = i; /* $e$ */
smax2[1] = 2*i+3; /* $\lambda$, $h$ */
sinfhdmh[i] = fn_clone(tmproveradmh[0], 2, smin2, smax2, NULL, NULL);
cosfhdmh[i] = fn_duplicate(sinfhdmh[i]);
}
sincosfh(sinfhdmh, cosfhdmh, wdmh);
smin2[1] = 0; /* $\lambda$, $h$ */
smax2[1] = 0; /* $\lambda$, $h$ */
for(i = 0; i < DEG_E; ++i) {
for(j = 0; j < DEG_E; ++j) {
smin2[0] = j; /* $|X|$ */
smax2[0] = j; /* $|X|$ */
e2Xdmh[i][j] = fn_clone(tmproveradmh[i], 2, smin2, smax2, NULL, NULL);
}
}
init_e2X(e2Xdmh);
smin5[1] = 0; /* $\lambda */
smin5[2] = 0; /* $D_2^*$ */
smax5[2] = 0; /* $D_2^*$ */
smin5[3] = 1; /* $a$ */
smax5[3] = 1; /* $a$ */
for(i = 0; i < DEG_E; ++i) {
smin5[0] = i; /* $|X|$, $h$ */
smax5[0] = i; /* $|X|$, $h$ */
smax5[1] = i; /* $\lambda */
roveradmh[i] = fn_clone(xiovera2dmh[0][0], 4, smin5, smax5, NULL, NULL);
}
smin5[1] = 0; /* $\lambda$ */
smin5[2] = 0; /* $D_2^*$ */
smax5[2] = 0; /* $D_2^*$ */
smin5[3] = 1; /* $a$ */
smax5[3] = 1; /* $a$ */
for(i = 0; i < DEG_E; ++i) {
smin5[0] = i; /* $|X|$, $h$ */
smax5[0] = i; /* $|X|$, $h$ */
smax5[1] = i+1; /* $\lambda$ */
rasin1dmh[i] = fn_clone(xiovera2dmh[0][0], 4, smin5, smax5, NULL, NULL);
racos1dmh[i] = fn_clone(xiovera2dmh[0][0], 4, smin5, smax5, NULL, NULL);
rasin2dmh[i] = fn_clone(xiovera2dmh[0][0], 4, smin5, smax5, NULL, NULL);
racos2dmh[i] = fn_clone(xiovera2dmh[0][0], 4, smin5, smax5, NULL, NULL);
}
for(i = 0; i < DEG_E; ++i) {
for(j = 0; j < DEG_E-i; ++j) {
if(i == 0) tmpvetdmh[i+j] = prod(tmproveradmh[i], sinfhdmh[j]);
else add_prod_to(tmproveradmh[i], sinfhdmh[j], tmpvetdmh[i+j]);
}
}
for(i = 0; i < DEG_E; ++i)
fn_remove(sinfhdmh[i]);
e2X(tmpvetdmh, e2Xdmh);
for(i = 0; i < DEG_E; ++i) {
two2five_a(tmpvetdmh[i], rasin1dmh[i]);
fn_remove(tmpvetdmh[i]);
}
for(i = 0; i < DEG_E; ++i) {
for(j = 0; j < DEG_E-i; ++j) {
if(i == 0) tmpvetdmh[i+j] = prod(tmproveradmh[i], cosfhdmh[j]);
else add_prod_to(tmproveradmh[i], cosfhdmh[j], tmpvetdmh[i+j]);
}
}
for(i = 0; i < DEG_E; ++i)
fn_remove(cosfhdmh[i]);
e2X(tmpvetdmh, e2Xdmh);
for(i = 0; i < DEG_E; ++i) {
two2five_a(tmpvetdmh[i], racos1dmh[i]);
fn_remove(tmpvetdmh[i]);
}
for(i = 0; i < DEG_E; ++i) {
inner2outer(rasin1dmh[i], rasin2dmh[i]);
inner2outer(racos1dmh[i], racos2dmh[i]);
}
e2X(tmproveradmh, e2Xdmh);
for(i = 0; i < DEG_E; ++i)
for(j = 0; j < DEG_E; ++j)
fn_remove(e2Xdmh[i][j]);
for(i = 0; i < DEG_E; ++i) {
two2five_a(tmproveradmh[i], roveradmh[i]);
fn_remove(tmproveradmh[i]);
}
smin5[1] = 0; /* $\lambda$ */
smin5[2] = 0; /* $D_2^*$ */
smax5[2] = 0; /* $D_2^*$ */
smin5[3] = 2; /* $L$ */
smax5[3] = 2; /* $L$ */
for(i = 0; i < DEG_E; ++i) {
smin5[0] = i; /* $|X|$, $h$ */
smax5[0] = i; /* $|X|$, $h$ */
smax5[1] = i+2; /* $\lambda$ */
xidmh[i] = fn_clone(xiovera2dmh[0][0], 4, smin5, smax5, NULL, NULL);
}
for(i = 0; i < DEG_E; ++i) {
for(j = 0; j < DEG_E-i; ++j) {
tmp0dmh = prod(roveradmh[i], roveradmh[j]);
tmp1dmh = fn_duplicate(tmp0dmh);
inner2outer(tmp0dmh, tmp1dmh);
add_to(tmp0dmh, xidmh[i+j]);
fn_remove(tmp0dmh);
add_to(tmp1dmh, xidmh[i+j]);
fn_remove(tmp1dmh);
}
}
for(i = 0; i < DEG_E; ++i)
fn_remove(roveradmh[i]);
for(i = 0; i < DEG_E; ++i) {
for(j = 0; j < DEG_E-i; ++j) {
tmp0dmh = prod(rasin1dmh[i], rasin2dmh[j]);
coef = -2.e0;
add_scalarmultiply_to(tmp0dmh, &coef, xidmh[i+j]);
fn_remove(tmp0dmh);
}
}
for(i = 0; i < DEG_E; ++i) {
fn_remove(rasin1dmh[i]);
fn_remove(rasin2dmh[i]);
}
for(i = 0; i < DEG_E; ++i) {
for(j = 0; j < DEG_E-i; ++j) {
tmp0dmh = prod(racos1dmh[i], racos2dmh[j]);
coef = -2.e0;
add_scalarmultiply_to(tmp0dmh, &coef, xidmh[i+j]);
fn_remove(tmp0dmh);
}
}
for(i = 0; i < DEG_E; ++i) {
fn_remove(racos1dmh[i]);
fn_remove(racos2dmh[i]);
}
for(i = 0; i < DEG_E; ++i) {
a2L(xidmh[i], i, xiovera2dmh);
}
tmp0dmh = fn_duplicate(xidmh[0]);
for(i = 0; i < DEG_E; ++i)
fn_remove(xidmh[i]);
get_fn_handlers(tmp0dmh, &fh);
for(i = 0; i < 9; ++i)
k5[i] = 0;
k5[1] = 2;
coef = -1.e0;
fh.add_coef(&coef, k5, tmp0dmh);
/* */
k5[1] = 0;
k5[2] = 2;
fh.add_coef(&coef, k5, tmp0dmh);
/* */
k5[1] = 1;
k5[2] = 1;
k5[5] = 1;
k5[6] = -1;
coef = 2.e0;
fh.add_coef(&coef, k5, tmp0dmh);
/* */
a2L(tmp0dmh, 0, xiovera2dmh);
fn_remove(tmp0dmh);
for(i = 0; i < DEG_M; ++i) {
for(j = 0; j < DEG_E; ++j) {
cutepsrel(xiovera2dmh[i][j], 1.e-10);
}
}
dm_shuffle(0);
#if DEBUG
printf("Exit xi_over_a2.\n");
#endif
return;
}

Variable Documentation

DOUBLE A_fix[2]

Fixed semi-major axes: $a_1^{*}$, $a_2^{*}$

Definition at line 274 of file threebody_plane.c.

DOUBLE Alpha_fix

Fixed semi-major axes ratio: $\alpha^*={a_1^*}/{a_2^*}$

Definition at line 275 of file threebody_plane.c.

DOUBLE Beta1

$\beta_1=m_0m_1/(m_0+m_1)$

Definition at line 272 of file threebody_plane.c.

DOUBLE Beta2

$\beta_2=m_0m_2/(m_0+m_2)$

Definition at line 273 of file threebody_plane.c.

DOUBLE D2_fix

The small parameter introduced by Laskar & Robutel: $D_2^{*}$

Definition at line 277 of file threebody_plane.c.

DOUBLE Lambda_fix[2]

Fixed fast actions: $\Lambda_1^{*}$, $\Lambda_2^{*}$

Definition at line 276 of file threebody_plane.c.

DOUBLE M0

Mass of the star: $m_0$

Definition at line 267 of file threebody_plane.c.

DOUBLE M1

Mass of the inner planet: $m_1$

Definition at line 268 of file threebody_plane.c.

DOUBLE M2

Mass of the outer planet: $m_2$

Definition at line 269 of file threebody_plane.c.

DOUBLE Mu1

$\mu_1=G(m_0+m_1)$

Definition at line 270 of file threebody_plane.c.

DOUBLE Mu2

$\mu_2=G(m_0+m_2)$

Definition at line 271 of file threebody_plane.c.