- #1
maves
- 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}];
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?