# Orbits in gravitational field of a toroidal planet

## 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: $\rho=const.$
2) Potential on a unit of mass (as the mass of the circulating object is irrelevant), $\vec{r}$ marks an arbitrary point in space and $\vec{R}$ a point on torus):
$$U = -G\rho \int_{Torus}\frac{1}{\left| {\vec{r}-\vec{R}} \right|}\,dV,$$

3) Acceleration is gradient of potential:
$$\ddot{x} =-\frac{\partial U}{\partial x} \\ \ddot{y} =-\frac{\partial U}{\partial y} \\ \ddot{z} =-\frac{\partial U}{\partial z}$$

4) Total energy:
$$E=T+V=\frac{1}{2}(\dot{x}^2+\dot{y}^2+\dot{z}^2)+U\,(x,y,z)$$

## The Attempt at a Solution

My coordinates are: arbitrary point on torus: $\vec R=(R\cos \phi, R\sin \phi, Z)$, arb. point in space: $\vec r=(x,y,z)$.
In my case, the expression for potential is
$$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} },$$
where $b$ is the inner radius and $b+c$ the outer radius of the torus.

And I will also need the derivatives of the potential, which are:
$$\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}}$$

So first of all, I draw the contour plot of potential in planes $xy$ and $xz$. 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 $xy$ and $xz$, 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}];
And then I try to find the solution of these DE by using NDSolve:
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]]
And draw a plot of what I get:
Code:
ParametricPlot[Evaluate[{x[t], y[t]} /. s], {t, 0, 10},
PlotRange -> All, AxesOrigin -> {0, 0},
AxesLabel -> {"x[t]", "y[t]"}]
I set my inner&outer radius for instance as 6 and 7, that's why I started with initial x (x[0]) close to 7 and then I started to look for a y'[0], which would give me a nice orbit. But it seems like something here is wrong, as I always get some kind of a ellipse (in best case), which even goes through the torus, which is, of course, physically impossible. Only when I make a very large y'[0], i get something, that looks like a circular orbit, but it seems like a "forced" solution, also when I calculate its energy, it is extremely large (over 3000 etc.), but I believe I should get sth small and negative (bounded states?). I also tried with a torus with a simpler symmetry (cross-section is a square), the results were similar. Where could the problem be?

Related Introductory Physics Homework Help News on Phys.org
I am not sure I understand the geometry of the problem, the bit about the cross-section. I do not understand your torroidal coordinates, either.

I am not sure I understand the geometry of the problem, the bit about the cross-section. I do not understand your torroidal coordinates, either.
I see. I don't think that the geometry here is so important, as I get similar results for square cross-section. I will append a few pictures for better understanding of the geometry.

The first picture shows you, how the coordinate system is set. The point inside a torus has coordinates (X,Y,Z) and I change them in into cylindrical, as described above, as it is easier to integrate over torus then. The point outside the torus, where the imaginary object is, has coordinates (x,y,z) and I do not put them into cylindrical coord., as it is easier to calculate the acceleration in the cartesian coord.sys.

The second picture shows the xz plane - the mentioned cross-section of the torus.

The third picture shows, how the torus looks.

#### Attachments

• 11.2 KB Views: 275
• 4.5 KB Views: 219
• 10.6 KB Views: 261
Last edited:
You should use cylindrical coordinates for the outer space, too. Why? Because the potential energy is axisymmetrical here, so it won't depend on $\phi$. This means you can represent the vector to the orbiting body this way: $\mathbf{r} = r\mathbf{e}_r + z\mathbf{e}_z$, where $\mathbf{e}_r$ and $\mathbf{e}_z$ are unit vectors along the r and z coordinates. Then if you do the math, you will get $$\ddot{r} - r \dot{\phi}^2 = - \frac {\partial U} {\partial r} \\ 2\dot{r}\dot {\phi} + r \ddot{\phi} = 0 \\ \ddot{z} = - \frac {\partial U} {\partial z}$$The first two equations are very similar to the equations of "normal" two-body problem, except that the potential energy term is different. The second equation is in fact the area law, $r^2\dot{\phi} = const$. The third equation is the harmonic oscillator. You also have conservation of energy for the entire system.

What I think you should do is the find the minimal points and curves of the potential, linearize the potential in the vicinity of the those locations, and solve the equations of motion with the linearized potentials.

^Thank you for your answer, I will try doing it that way, and ask again, if another problem occurs.
And by the way, a little problem with the sign, that bothers me, and I don't understand, where I have missed it:

If the potential is $$U\,(r,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{r^2+R^2-2Rr\cos{\phi} +(z-Z)^2 } }$$

and the acceleration is $-\nabla U$, so wouldn't for instance the radial component be
$$\ddot{r}-r\dot{\varphi}^2 = + G\rho \int_{0}^{2\pi}\int_{b}^{b+c}\int_{R-b-c}^{-R+b+c} \frac{R(r-R\cos{\phi}) \,d\phi \,dR\,dZ}{\sqrt[3]{r^2+R^2-2Rr\cos{\phi} +(z-Z)^2 }},$$
as the "-" before gradient "destroys" the one in potential? But the force must be attractive, which means negative (as regards r) What have I overlooked?

Last edited:
I think you forgot the -1/2 factor when differentiating the integral.

I think you forgot the -1/2 factor when differentiating the integral.
Aaaargh... of course! :) Thanks.

Hello again,
maybe another stupid question, but still - what exactly should I consider as orbits? For instance, I found some nice rosettes in xy plane (see in attachment), but when I calculate their energy (E=T+U), they don't seem to be constant in every point.

Furthermore, I am positive there should be some orbits, that only circulate around "one side of the planet", when looking at the xz plane, but no matter how I try to find them, the trajectory always tends to also go to the other side and forms all kinds of ellipses, which is in my opinion more energy consuming.

And, something that also bothers me, is the definition of energy - to be negative when we have bound states ... all my calculated energies are positive ... :S

#### Attachments

• 19.5 KB Views: 294
I am unsure how you could get non-constant energy trajectories. The force is potential, so energy must be conserved.

As for the negative energy. Consider the orbit where z = 0 and R = const. Is that a possible orbit? It seems it should be, given the symmetries in the problem.

Yes, I have found an orbit like that, but the professor also advised me to look for the orbits, that go around the torus like a ring (xz plane) - to compare their energies. But my trajectories always go around the whole torus.

Edit: I attached a file of an example, that seems pretty illogical to me. I deliberately started on the inner side of the torus, to "force" the trajectory to go over the triangle, but it went the other way ... or is there something wrong with my understanding?

#### Attachments

• 7 KB Views: 230
Last edited:
Have you tried to plot the potential? What does it look like?

Potential in xy plane are concentric circles, potential in xz plane looks like this ... (the second picture is zoomed and shows, where the torus is located)

#### Attachments

• 16.6 KB Views: 225
• 22.5 KB Views: 224
In the xz plane, due to the symmetry, you could treat the problem as 2D (forces in the y direction should cancel out). Then what you would have is a special kind of the three-body problem. It goes by many names, one of them is the problem of two fixed centers. This problem is known to admit closed orbits only when initial conditions are specially selected, otherwise the trajectory fills some portion of the plane densely. It is also known that motion about just one center is possible (which should be intuitively clear, as the Earth - Sun system can be regarded as an example of two fixed centers).

Thank you very much for your help, I will consider this idea.
And about the energies ... not that they are not constant, they actually fluctuate around one value, but the amplitude is mostly quite big ... could this be only due to numerical errors? :S

I also figured out, that if I enlarge the inner radius of the torus, the potential reaches much larger (negative) values, which is kinda logical, so i thought maybe then i would get negative full energies (which was a childish attempt to avoid unwanted results), but with these dimensions, the evaluation takes muuuch more time and i think it's practically impossible to find sth practical in a decent period of time :S Ohh, numerical methods ...

#### Attachments

• 8.1 KB Views: 269
Oscillating total energy is definitely not right. The force is potential, so total energy must constant. There is probably an error in the computation of the derivatives, or in the equations of motion.

Have you considered studying some approximation to the problem?

I tried to approximate the integrands with a series and integrate them after that ... and it turned out, that Mathematica 8 integrates in a different way than the earlier versions of the program (I tested this on an expression I evaluated some time ago when I had M.6 or 7 and it put out a completely different expression, including some wild complex numbers ... )
:S :S :S

You should try the multipole approximation. This is a standard method in astrodynamics.

What you could also do, is approximate your torus with an infinitely thin ring. In this case, I think, the potential could be integrated analytically.

Regardless, you should figure out what error results in oscillating energy and fix that.