<p>I would like to calculate Mercury's orbit in Mathematica. My initial conditions, in AU and AU/year, are:</p>
<pre><code>x = (-3.610096904077567*10^-1); y = (-2.285553971112067*10^-1); z = (-8.467301297840119*10^-2);

vx = (9.95457152429289*10^-3*365.25); vy = (-1.918192442772782*10^-2*365.25); vz =(-1.127871068908956*10^-2*365.25);
</code></pre>
<p>I express M,G,c in Msun,AU,years:</p>
<pre><code>M := 1;
G := 4 \Pi^2;
c := 63197.7909261;
</code></pre>
<p>If I solve Newton equations, everything works and I obtain the correct orbit. But If I try with relativistic eqns:</p>
<pre><code>sol = NDSolve[{(
c^2 r[\[Tau]] (-c^2 + G M r[\[Tau]]) Derivative[1][r][\[Tau]]^2 + 
 G M (c^2 - 2 G M r[\[Tau]])^2 Derivative[1][t][\[Tau]]^2)/(
c^2 - 2 G M r[\[Tau]]) == 
c^2 r[\[Tau]]^2 (r^\[Prime]\[Prime])[\[Tau]],  Sin[\[Theta][\[Tau]]] (2 Cos[\[Theta][\[Tau]]] Derivative[
     1][\[Theta]][\[Tau]] Derivative[1][\[Phi]][\[Tau]] + 
   Sin[\[Theta][\[Tau]]] (\[Phi]^\[Prime]\[Prime])[\[Tau]]) == 0,Cos[\[Theta][\[Tau]]] Sin[\[Theta][\[Tau]]] Derivative[
  1][\[Phi]][\[Tau]]^2 == (\[Theta]^\[Prime]\[Prime])[\[Tau]], 2 G M Derivative[1][r][\[Tau]] Derivative[1][
  t][\[Tau]] == (c^2 - 2 G M r[\[Tau]]) (
  t^\[Prime]\[Prime])[\[Tau]], r[0] == 0.4355859102907327`, r'[0] == 1.46359, \[Phi][0] == -2.5771894914147078`, \[Phi]'[0] == 
18.406053400410222`, \[Theta][0] == 
1.766430630766984`, \[Theta]'[0] == 8.975543771872271`,  t[0] == 0.00000000, t'[0] == 0.0003188110212866474`}, {r, \[Phi], \[Theta], t}, {\[Tau], 0, 1}]
</code></pre>
<p>where initial conditions are converted from cartesian to spherical, I don't obtain the right orbit:</p>
<pre><code>\[Phi]s[tt_] = \[Phi][tt] /. sol[[1]] /. tt -&gt; \[Tau]
rs[tt_] = r[tt] /. sol[[1]] /. tt -&gt; \[Tau]
\[Theta]s[tt_] = \[Theta][tt] /. sol[[1]] /. tt -&gt; \[Tau]

ParametricPlot[{(rs[tt]*Sin[\[Theta]s[tt]]*Cos[\[Phi]s[tt]]),    (rs[tt]*
Sin[\[Theta]s[tt]]*Sin[\[Phi]s[tt]])}, {\[Tau], 0, 1}]
</code></pre>
<p>I don't understand If there is something wrong with the initial conditions, because I used Newton equations with the same procedure and units and everything worked.</p>

CC BY-SA 4.0
2021-05-03T13:13:09.740Z
Hi Mark8991, the way you phrase your question, it is not really suited for Astronomy.SE, as it seems purely a matter of coding, and might thus be better for StackOverflow. If you believe your problem has to do with astronomy, you should show your _equations_, rather than your code. You should also elaborate on what you mean by "not obtaining the right orbit", I think.
pela
CC BY-SA 4.0
2021-05-03T13:16:49.047Z
There's also a Stack site dedicated to Mathematica: https://mathematica.stackexchange.com/
PM 2Ring
CC BY-SA 4.0
2021-05-03T13:24:22.727Z
Thank you, so, I checked my equations and I used the same code for other bodied and it worked. Maybe I'm writing the wrong initial conditions to pass from cartesian to spherical coordinates? With wrong orbit I mean that I obtain a sort of spiraling orbit, If I put r'=0 it closes. But the right r'[0] should be different from 0.
Mark8991
CC BY-SA 4.0
2021-05-03T13:26:18.440Z
You should also explain what your goal is here. I suppose you're wanting to see Mercury's anomalous apsidal precession of 43 arc-seconds per century. Note that the observed precession is around 574"/century, so if you want an accurate Mercury orbit, you need to account for things like the oblateness of the Sun, and perturbations from planets (esp, Venus, Earth & Jupiter).
PM 2Ring
CC BY-SA 4.0
2021-05-03T13:30:29.870Z
For now I just put relativistic equations coming from Schwarzschild Lagrangian
Mark8991
CC BY-SA 4.0
2021-05-03T17:11:00.603Z
What do you get if you put the speed of light much higher in your relativistic equation (like a factor 100 or more)?
Thomas
CC BY-SA 4.0
2021-05-03T20:26:27.737Z
I see changes if I put r'[0]=0, an orbit which seems the right one
Mark8991
<orbit>
Wrong planet orbit
Mark8991
# Wrong planet orbit

(Asked by: username_0 on 2021-05-03)

I would like to calculate Mercury's orbit in Mathematica. My initial conditions, in AU and AU/year, are:
x = (-3.610096904077567*10^-1); y = (-2.285553971112067*10^-1); z = (-8.467301297840119*10^-2);

vx = (9.95457152429289*10^-3*365.25); vy = (-1.918192442772782*10^-2*365.25); vz =(-1.127871068908956*10^-2*365.25);

I express M,G,c in Msun,AU,years:
M := 1;
G := 4 \Pi^2;
c := 63197.7909261;

If I solve Newton equations, everything works and I obtain the correct orbit. But If I try with relativistic eqns:
sol = NDSolve[{(
c^2 r[\[Tau]] (-c^2 + G M r[\[Tau]]) Derivative[1][r][\[Tau]]^2 + 
 G M (c^2 - 2 G M r[\[Tau]])^2 Derivative[1][t][\[Tau]]^2)/(
c^2 - 2 G M r[\[Tau]]) == 
c^2 r[\[Tau]]^2 (r^\[Prime]\[Prime])[\[Tau]],  Sin[\[Theta][\[Tau]]] (2 Cos[\[Theta][\[Tau]]] Derivative[
     1][\[Theta]][\[Tau]] Derivative[1][\[Phi]][\[Tau]] + 
   Sin[\[Theta][\[Tau]]] (\[Phi]^\[Prime]\[Prime])[\[Tau]]) == 0,Cos[\[Theta][\[Tau]]] Sin[\[Theta][\[Tau]]] Derivative[
  1][\[Phi]][\[Tau]]^2 == (\[Theta]^\[Prime]\[Prime])[\[Tau]], 2 G M Derivative[1][r][\[Tau]] Derivative[1][
  t][\[Tau]] == (c^2 - 2 G M r[\[Tau]]) (
  t^\[Prime]\[Prime])[\[Tau]], r[0] == 0.4355859102907327`, r'[0] == 1.46359, \[Phi][0] == -2.5771894914147078`, \[Phi]'[0] == 
18.406053400410222`, \[Theta][0] == 
1.766430630766984`, \[Theta]'[0] == 8.975543771872271`,  t[0] == 0.00000000, t'[0] == 0.0003188110212866474`}, {r, \[Phi], \[Theta], t}, {\[Tau], 0, 1}]

where initial conditions are converted from cartesian to spherical, I don't obtain the right orbit:
\[Phi]s[tt_] = \[Phi][tt] /. sol[[1]] /. tt -> \[Tau]
rs[tt_] = r[tt] /. sol[[1]] /. tt -> \[Tau]
\[Theta]s[tt_] = \[Theta][tt] /. sol[[1]] /. tt -> \[Tau]

ParametricPlot[{(rs[tt]*Sin[\[Theta]s[tt]]*Cos[\[Phi]s[tt]]),    (rs[tt]*
Sin[\[Theta]s[tt]]*Sin[\[Phi]s[tt]])}, {\[Tau], 0, 1}]

I don't understand If there is something wrong with the initial conditions, because I used Newton equations with the same procedure and units and everything worked.

**Tags:** orbit
> **Comments:**
> **username_1**: Hi Mark8991, the way you phrase your question, it is not really suited for Astronomy.SE, as it seems purely a matter of coding, and might thus be better for StackOverflow. If you believe your problem has to do with astronomy, you should show your _equations_, rather than your code. You should also elaborate on what you mean by "not obtaining the right orbit", I think.  
> **username_2**: There's also a Stack site dedicated to Mathematica: https://mathematica.stackexchange.com/  
> **username_0**: Thank you, so, I checked my equations and I used the same code for other bodied and it worked. Maybe I'm writing the wrong initial conditions to pass from cartesian to spherical coordinates? With wrong orbit I mean that I obtain a sort of spiraling orbit, If I put r'=0 it closes. But the right r'[0] should be different from 0.  
> **username_2**: You should also explain what your goal is here. I suppose you're wanting to see Mercury's anomalous apsidal precession of 43 arc-seconds per century. Note that the observed precession is around 574"/century, so if you want an accurate Mercury orbit, you need to account for things like the oblateness of the Sun, and perturbations from planets (esp, Venus, Earth & Jupiter).  
> **username_0**: For now I just put relativistic equations coming from Schwarzschild Lagrangian  
> **username_3**: What do you get if you put the speed of light much higher in your relativistic equation (like a factor 100 or more)?  
> **username_0**: I see changes if I put r'[0]=0, an orbit which seems the right one  

---

## 0 answers

