Available with Zymplectic v.0.1.0

The Two-body problem is about predicting the motion of two interacting objects. Here we will consider two heavenly bodies referred to as objects with masses **m**_{0} and **m**_{1} gravitationally bound by the potential V. The energy of the system is denoted by the constant Hamiltonian H equal to the sum of the kinetic energy T and potential energy V:
$$H(\mathbf{q},\mathbf{p}) = T(\mathbf{p}) + V(\mathbf{q})$$
$$T(\mathbf{p}) = \frac{\mathbf{p}_{object 0}^2}{2m_0} + \frac{\mathbf{p}_{object 1}^2}{2m_1}$$
$$V(\mathbf{q}) = -G\frac{m_0m_1}{|\mathbf{r}|}$$
where **q** and **p** are vectors of canonical coordinates for position and momenta respectively, G is the gravitational constant and **r** = **q**_{2} - **q**_{1}. Note that this Hamiltonian system is separable because T is only a function of **p** while V is only a function of **q**.
Four implementations of the two-dimensionel two-body problem are described in the following section.
#### Two-body problem in Cartesian coordinates (system 1)

First obtain:
$$\frac{\partial V}{\partial \mathbf{q}} = G\frac{m_0m_1\mathbf{r}}{|\mathbf{r}|^3}$$
$$\frac{\partial T}{\partial \mathbf{p}} = \frac{\mathbf{p}_{object0}}{m_0} + \frac{\mathbf{p}_{object1}}{m_1}$$
Given a set of appropriate constants, we consider the practical implementation of this system in the code below.
The kinetic energy of the system can be split into two parts. The contribution of the first object 0 is introduced by the connector "Z.T(T,dT,0,0)" and the second object 1 by "Z.T(T,dT,1,1)", where T is the energy function and dT is the gradient function. The calculation of the potential energy requires object 0 and object 1 and therefore requires the connector Z.V(V,dV,0,1); The function T (ignore the argument) is a function returning the kinetic energy of either object whereas dT calculates the gradient with components d1 and d2. Note that obj1 refers to the first connected object in Z.T(T,dT,obj1,obj2). The V and dV functions utilize the coordinates of both objects with q and d referring to obj1 and Q and D referring to obj2.

```
//Copy-paste into Zymplectic or load as file to run the simulation
const double G = 1.0; //gravitational constant
const double m[] = {1.0,5.0}; //mass of object 0 and 1
#include "math.h"
double T(Arg) {
return 0.5*(p1*p1 + p2*p2)/m[obj1];
}
void dT(Arg) {
d1 = p1/m[obj1];
d2 = p2/m[obj1];
}
double V(Arg) {
double dx = q1 - Q1;
double dy = q2 - Q2;
return -G*m[obj1]*m[obj2]/sqrt(dx*dx + dy*dy);
}
void dV(Arg) {
double dx = q1 - Q1;
double dy = q2 - Q2;
double d = sqrt(dx*dx + dy*dy);
d = G*m[obj1]*m[obj2]/(d*d*d);
d1 = dx*d; d2 = dy*d; //object 0
D1 =-dx*d; D2 =-dy*d; //object 1
}
double _q_init[] = { //initial conditions
5.0, 0.0, //object 0
0.0, 0.0}; //object 1
double _p_init[] = {
0.0,-1.0,
0.0, 1.0};
void main() {
Z.V(V,dV,0,1);
Z.T(T,dT,0,0);
Z.T(T,dT,1,1);
}
```

```
//Copy-paste into Zymplectic or load as file to run the simulation
const double G = 1.0; //gravitational constant
const double m[] = {1.0,5.0}; //mass of object 0 and 1
const double mu = m[0]*m[1]/(m[0] + m[1]); //reduced mass
const double M = m[0] + m[1];
#include "math.h"
double T(Arg) {
return 0.5*(p1*p1 + p2*p2)/mu;
}
void dT(Arg) {
d1 = p1/mu;
d2 = p2/mu;
}
double V(Arg) {
return -G*mu*M/sqrt(q1*q1 + q2*q2);
}
void dV(Arg) {
double d = sqrt(q1*q1 + q2*q2);
d = G*mu*M/(d*d*d);
d1 = q1*d;
d2 = q2*d;
}
double _q_init[] = {5.0, 0.0}; //initial conditions
double _p_init[] = {0.0,-1.0};
void main() {
Z.V(V,dV,0);
Z.T(T,dT,0);
}
```

```
//Copy-paste into Zymplectic or load as file to run the simulation
const double G = 1.0; //gravitational constant
const double m0 = 1.0; //mass of object 1
const double m1 = 5.0; //mass of object 2
const double mu = m0*m1/(m0 + m1); //reduced mass
double H(Arg) {
return 0.5*(p1*p1 + (p2*p2/(q1*q1)))/mu - G*mu*(m0 + m1)/q1; //for positive q1
}
void dH(Arg) {
dq1 = (G*q1*mu*mu*(m0 + m1) - p2*p2)/(q1*q1*q1*mu);
dq2 = 0;
dp1 = p1/mu;
dp2 = p2/(q1*q1*mu);
}
double _q_init[] = {2.0,0.0}; //radius,angle
double _p_init[] = {0.0,2.0};
double GLX[2],GLY[2];
#include "math.h"
void DrawFunc(){
GLX[0] = 0;
GLY[0] = 0;
GLX[1] = obj_q[0][0]*cos(obj_q[0][1]);
GLY[1] = obj_q[0][0]*sin(obj_q[0][1]);
draw.points(GLX,GLY,2);
}
void main() {
Z.OnDraw = DrawFunc; //replaces coordinate plots of polar coordinates with customized Cartesian coordinates
Z.xmin = -2; //set visuel range
Z.xmax = 4;
Z.ymin = -3;
Z.ymax = 3;
Z.H(H,dH,0);
}
```

```
//Copy-paste into Zymplectic or load as file to run the simulation
#include "math.h"
const double G = 1.0;
const double m0 = 1.0;
const double m1 = 1.0;
const double mu = m0*m1/(m0 + m1); //reduced mass
double energy_H(Arg) {
return -G*mu*(m0 + m1)/sqrt(q1*q1 + q2*q2) + 0.5*(p1*p1 + p2*p2)/mu;
}
void dH(Arg) {
double d = sqrt(q1*q1 + q2*q2);
double f = q1*q1 + q2*q2; //f is the time scale factor
d = G*mu*(m0 + m1)/(d*d*d);
dq1 = q1*d*f; //time factor is 1 at apoapsis (based on _q_init)
dq2 = q2*d*f; //time factor becomes small at close approaches
dp1 = p1/mu*f;
dp2 = p2/mu*f;
}
double _q_init[] = {1.0,0.0};
double _p_init[] = {0.0,0.2};
void main() { //this function is launched when executing Zymplectic.exe
Z.H(energy_H,dH,0,0);
Z.delay = 2.0;
Z.track = 1;
}
```