- #1

- 17

- 0

## Homework Statement

We have a homogenous toroidal planet. Its cross-section is half of a square, so that the hypotenuse is parallel to the z-axis. Find the orbits of constant (especially lowest) energy around that planet.

## Homework Equations

1) Homogenous planet: [itex]\rho=const. [/itex]

2) Potential on a unit of mass (as the mass of the circulating object is irrelevant), [itex]\vec{r}[/itex] marks an arbitrary point in space and [itex]\vec{R}[/itex] a point on torus):

[tex]U = -G\rho \int_{Torus}\frac{1}{\left| {\vec{r}-\vec{R}} \right|}\,dV, [/tex]

3) Acceleration is gradient of potential:

[tex] \ddot{x} =-\frac{\partial U}{\partial x} \\

\ddot{y} =-\frac{\partial U}{\partial y} \\

\ddot{z} =-\frac{\partial U}{\partial z} [/tex]

4) Total energy:

[tex]E=T+V=\frac{1}{2}(\dot{x}^2+\dot{y}^2+\dot{z}^2)+U\,(x,y,z)[/tex]

## The Attempt at a Solution

My coordinates are: arbitrary point on torus: [itex]\vec R=(R\cos \phi, R\sin \phi, Z)[/itex], arb. point in space: [itex]\vec r=(x,y,z)[/itex].

In my case, the expression for potential is

[tex]U\,(x,y,z) = -G\rho \int_{0}^{2\pi}\int_{b}^{b+c}\int_{R-b-c}^{-R+b+c} \frac{R\,d\phi \,dR\,dZ}{\sqrt{(x-R\cos{\phi} )^2+( y-R\sin{\phi} )^2+( z-Z )^2} },[/tex]

where [itex]b[/itex] is the inner radius and [itex]b+c[/itex] the outer radius of the torus.

And I will also need the derivatives of the potential, which are:

[tex]\ddot{x}=- G\rho \int_{0}^{2\pi}\int_{b}^{b+c}\int_{R-b-c}^{-R+b+c} \frac{R(x-R\cos{\phi}) \,d\phi \,dR\,dZ}{\sqrt[3]{(x-R\cos{\phi})^2+( y-R\sin{\phi} )^2+( z-Z )^2}} \\ \\

\ddot{y}= -G\rho \int_{0}^{2\pi}\int_{b}^{b+c}\int_{R-b-c}^{-R+b+c} \frac{R(y-R\sin{\phi}) \,d\phi \,dR\,dZ}{\sqrt[3]{(x-R\cos{\phi})^2+( y-R\sin{\phi} )^2+( z-Z )^2}} \\ \\

\ddot{z} = - G\rho \int_{0}^{2\pi}\int_{b}^{b+c}\int_{R-b-c}^{-R+b+c} \frac{R(z-Z)\,d\phi \,dR\,dZ}{\sqrt[3]{(x-R\cos{\phi} )^2+( y-R\sin{\phi} )^2+( z-Z )^2}} [/tex]

So first of all, I draw the contour plot of potential in planes [itex]xy[/itex] and [itex]xz[/itex]. I discovered, that the potential has lowest values near to torus, so if I want lowest energy, it is best to have initial conditions somewhere near torus. So what I have to do next is to solve the differential equations under 3), draw a graph of my solutions to see, how they look, and see, which orbit has lowest energy (eq.4)).

Presumably, orbits of lowest energy will lie in planes (it is less "energy consuming"-? ), so all I have to do is to examine planes [itex]xy[/itex] and [itex]xz[/itex], which means, that I will always have a set of 2 DE, and the third coordinate will be zero.

OK, so I tried solving this in Mathematica.

Firstly, I defined the integrals. I had to integrate numerically.

This "?NumericQ" must be here to keep the program from evaluating the integral before you give him any values (which I will get later, from NDSolve).

For xy plane:

Code:

```
Clear[x, y, z, t];
int1[x_?NumericQ, y_?NumericQ] :=
NIntegrate[(-R (x - R Cos[\[Phi]]))/((x - R Cos[\[Phi]])^2 + (y - R Sin[\[Phi]])^2 + (0 - Z)^2)^(3/2), {\[Phi], 0, 2 Pi}, {R, b, b + c}, {Z, R - b - c, -R + b + c}];
int2[x_?NumericQ, y_?NumericQ] :=
NIntegrate[(-R (y - R Sin[\[Phi]]))/((x - R Cos[\[Phi]])^2 + (y - R Sin[\[Phi]])^2 + (0 - Z)^2)^(3/2), {\[Phi], 0, 2 Pi}, {R, b, b + c}, {Z, R - b - c, -R + b + c}];
```

Code:

```
s = NDSolve[{
x''[t] == int1[x[t], y[t]],
y''[t] == int2[x[t], y[t]],
x[0] == 7.1, x'[0] == 0, y[0] == 0, y'[0] == 1.5}, {x, x', y,
y'}, {t, 0, 10}][[1]]
```

Code:

```
ParametricPlot[Evaluate[{x[t], y[t]} /. s], {t, 0, 10},
PlotRange -> All, AxesOrigin -> {0, 0},
AxesLabel -> {"x[t]", "y[t]"}]
```