<p>If I follow the code, you are using Kepler's formulae to calculate the distance and speed of the planet, and then increasing the planet's position using speed × time.  That is essentially a form of Euler's method for solving numerical differential equations, and as such it may be sensitive to initial conditions and rounding.</p>
<p>You've chosen a small step size of one second.  Small steps can reduce one source of error (the essential error in Euler methods of approximation of a continuous process by a discrete one) but can introduce another error: If you you add a very small amount to a large amount, you can lose a lot of significant figure accuracy  (eg to 10 sig fig,  12345.00000 + 0.00001234567890 = 12345.00001, you would have lost all but the first significant figure. Now consider what would happen with <code>n=12345; for i in range 10000 { n+= 1.23456789e-5}</code> compared with <code>n+= 1.23456789e-5 * 10000</code>)</p>
<p>I suspect that the issue here is a matter of programming, rather than with the formulae.  Standard good practice is to do all solar system programming  not in SI units, but in AU-SolarMass-Years units.  These units can help avoid some of the issues with very small values, They make many of the constants, such as <span class="math-container">$G$</span> or <span class="math-container">$\mu$</span> more manageable (indeed in these units <span class="math-container">$G=\mu=4\pi^2$</span>)</p>
<p>Moreover I wonder about the general method.  You are using Kepler's laws, so why not just use Kepler's law for the orbital period?  Or, if you want to investigate orbital dynamics, why not use Newton's law of Gravity?</p>

CC BY-SA 4.0
2022-01-03T06:59:31.420Z
From what I've read, solar system programming used to use AU as the unit of distance, solar mass as the unit of mass, and a day (not a year) as the unit of time. This practice dates back to Gauss. The current practice uses the AU as a convenience distance whose value is defined as exactly 149597870700 meters rather than a calculated value based on the Gaussian gravitational constant.
David Hammen
CC BY-SA 4.0
2022-01-02T17:18:36.360Z
2022-01-02T21:53:19.440Z
2022-01-02T21:53:19.440Z
James K
<p>Too long for a comment, I refer you back to <a href="https://astronomy.stackexchange.com/a/48037/7982">@JamesK's answer</a> for the answer.</p>
<p>In Space Exploration SE See:</p>
<ul>
<li><a href="https://space.stackexchange.com/q/22915/12102">Calculating the planets and moons based on Newtons's gravitational force</a></li>
<li><a href="https://space.stackexchange.com/q/23408/12102">How to calculate the planets and moons beyond Newtons's gravitational force?</a>
<ul>
<li>Use a reasonable ODE solver; <em>at least</em> <a href="https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods" rel="nofollow noreferrer">RK4</a> (the classic Runge-Kutta method) or better. Do not use just the <a href="https://en.wikipedia.org/wiki/Euler_method" rel="nofollow noreferrer">Euler Method</a></li>
</ul>
</li>
<li><a href="https://space.stackexchange.com/q/43726/12102">Two Body Kepler Propogation: Errors from using float / single precision</a></li>
<li><a href="https://space.stackexchange.com/a/55906/12102">this answer to <em>Differences between numerical propagators</em></a> which links to interesting and helpful resources!</li>
</ul>
<p>Also see:</p>
<ul>
<li><a href="https://en.wikipedia.org/wiki/Semi-implicit_Euler_method" rel="nofollow noreferrer">Semi-implicit Euler method</a> as mentioned in <a href="https://space.stackexchange.com/questions/23408/how-to-calculate-the-planets-and-moons-beyond-newtonss-gravitational-force?noredirect=1&amp;lq=1#comment103549_23409">@DavidHammen's comment</a></li>
</ul>

CC BY-SA 4.0
2022-01-03T00:04:10.183Z
2022-01-03T00:04:10.183Z
uhoh
<p>Question:</p>
<blockquote>
<p>I am trying to make a simple java program to do planetary orbital calculation like how many days it takes to complete the orbital along with each Earth day data but I end up getting a 4 days error. Can anybody explain why?</p>
</blockquote>
<p>After fixing some problems after @PM2Ring's suggestion on other site.</p>
<p>I have this code</p>
<pre><code>        double semiMajorAxis = 1.5237 * 149598000;
        double eccentricty = 0.0934;
        double current_time = 0;
        double last_speed;
        double angular_distance = 0; // unit : rad
        boolean stop = false;

        while (!stop) {
            double radius = get_radius(angular_distance, eccentricty, semiMajorAxis);
            last_speed = (get_speed(radius * 1000, semiMajorAxis * 1000) / radius);
            angular_distance += last_speed;
            distance_in_km = angular_distance * radius;
            current_time += 1;
            
            if (angular_distance &gt;= Math.PI * 2) {
                System.out.println(&quot;Total days : &quot; + current_time / 60 / 60 / 24);
                stop = true;
            }
           }
</code></pre>
<p>Functions used in the above program-</p>
<pre><code>   public static double get_radius(double theeta, double eccentricty, double a) {
        // modern kepler formula
        return a * (1 - Math.pow(eccentricty, 2)) / (1 + eccentricty * Math.cos(theeta));
    }

    public static double get_speed(double r, double semiMajorAxis) {
        double mu = 1.32712440042 * Math.pow(10, 20);
        // vis-viva equation
        return (Math.sqrt(mu * ((2 / r) - (1 / semiMajorAxis))) / 1000);
    }
</code></pre>
<p><strong>What am I trying to do?</strong> - I am trying to calculate time taken for each planet to complete its revolution and also collect daily data for it, in this case its the Mars. To do
this I am using 1 second step to calculate angles displaced by using <a href="https://en.wikipedia.org/wiki/Vis-viva_equation" rel="nofollow noreferrer">vis-viva</a> to calculate speed and kepler formula to get radius.</p>
<p>However seems like there is a problem as the output I am getting for Mars is <code>685.4825000000001</code> <strong>days</strong> which is nearly 2 days behind of actual time i.e <code>687</code> <strong>days</strong> and for Jupiter I am getting <code>4334.061006944444</code> <strong>days</strong> which is nearly 2 day ahead of actual time <code>4332</code> <strong>days</strong>.</p>
<p>Can anybody help me fix the error or could guide me with some other way. It will really be helpfull.</p>

CC BY-SA 4.0
2022-01-02T16:51:28.703Z
You should go and delete the physics question to avoid cross posting
James K
CC BY-SA 4.0
2022-01-02T16:53:19.303Z
@JamesK ok done
Stars Tracker
CC BY-SA 4.0
2022-01-02T18:07:45.280Z
Check out integration techniques / methods and their sensitivity to different error types. Maybe start with something like https://www.ias.edu/sites/default/files/sns/Numerical-integration(2).pdf
planetmaker
CC BY-SA 4.0
2022-01-03T07:36:47.863Z
What you are doing is invalid. You are implicitly assuming a circular orbit in your calculation of angular velocity.
David Hammen
CC BY-SA 4.0
2022-01-03T07:57:08.097Z
To add to the above, your calculation always overestimates $\dot\theta$ due to the fact that Mars's orbit is eccentric, and I suspect that that is the key reason you are two days shy of the correct answer for Mars. Getting the right math is a bit ugly, but not overly so. For Jupiter, I suspect the key reason is that your one second time step is far too short. It's too short by a long shot. A good integrator (which you are not using) should take at most a few hundred steps per orbit.
David Hammen
CC BY-SA 4.0
2022-01-03T08:37:29.313Z
Thanks for all of the comments, I am trying to fix it by couple of answers below however @David Hammen thanks for your comment also i would like to add that i am not assuming circluar orbit
Stars Tracker
CC BY-SA 4.0
2022-01-03T12:13:17.453Z
Yes, you are assuming that, in calculating angular velocity as `last_speed = (get_speed(radius * 1000, semiMajorAxis * 1000) / radius)`. That is only valid for a circular orbit.
David Hammen
<orbit><astrophysics><software>
Incorrect Mars orbital time while calculating using Kepler formula in Java
Stars Tracker
# Incorrect Mars orbital time while calculating using Kepler formula in Java

(Asked by: username_0 on 2022-01-02)

Question:

I am trying to make a simple java program to do planetary orbital calculation like how many days it takes to complete the orbital along with each Earth day data but I end up getting a 4 days error. Can anybody explain why?

After fixing some problems after @PM2Ring's suggestion on other site.
I have this code
        double semiMajorAxis = 1.5237 * 149598000;
        double eccentricty = 0.0934;
        double current_time = 0;
        double last_speed;
        double angular_distance = 0; // unit : rad
        boolean stop = false;

        while (!stop) {
            double radius = get_radius(angular_distance, eccentricty, semiMajorAxis);
            last_speed = (get_speed(radius * 1000, semiMajorAxis * 1000) / radius);
            angular_distance += last_speed;
            distance_in_km = angular_distance * radius;
            current_time += 1;
            
            if (angular_distance >= Math.PI * 2) {
                System.out.println("Total days : " + current_time / 60 / 60 / 24);
                stop = true;
            }
           }

Functions used in the above program-
   public static double get_radius(double theeta, double eccentricty, double a) {
        // modern kepler formula
        return a * (1 - Math.pow(eccentricty, 2)) / (1 + eccentricty * Math.cos(theeta));
    }

    public static double get_speed(double r, double semiMajorAxis) {
        double mu = 1.32712440042 * Math.pow(10, 20);
        // vis-viva equation
        return (Math.sqrt(mu * ((2 / r) - (1 / semiMajorAxis))) / 1000);
    }

What am I trying to do? - I am trying to calculate time taken for each planet to complete its revolution and also collect daily data for it, in this case its the Mars. To do
this I am using 1 second step to calculate angles displaced by using vis-viva to calculate speed and kepler formula to get radius.
However seems like there is a problem as the output I am getting for Mars is 685.4825000000001 days which is nearly 2 days behind of actual time i.e 687 days and for Jupiter I am getting 4334.061006944444 days which is nearly 2 day ahead of actual time 4332 days.
Can anybody help me fix the error or could guide me with some other way. It will really be helpfull.

**Tags:** orbit, astrophysics, software
> **Comments:**
> **username_1**: You should go and delete the physics question to avoid cross posting  
> **username_0**: @JamesK ok done  
> **username_2**: Check out integration techniques / methods and their sensitivity to different error types. Maybe start with something like https://www.ias.edu/sites/default/files/sns/Numerical-integration(2).pdf  
> **username_3**: What you are doing is invalid. You are implicitly assuming a circular orbit in your calculation of angular velocity.  
> **username_3**: To add to the above, your calculation always overestimates $\dot\theta$ due to the fact that Mars's orbit is eccentric, and I suspect that that is the key reason you are two days shy of the correct answer for Mars. Getting the right math is a bit ugly, but not overly so. For Jupiter, I suspect the key reason is that your one second time step is far too short. It's too short by a long shot. A good integrator (which you are not using) should take at most a few hundred steps per orbit.  
> **username_0**: Thanks for all of the comments, I am trying to fix it by couple of answers below however @David Hammen thanks for your comment also i would like to add that i am not assuming circluar orbit  
> **username_3**: Yes, you are assuming that, in calculating angular velocity as `last_speed = (get_speed(radius * 1000, semiMajorAxis * 1000) / radius)`. That is only valid for a circular orbit.  

---

## 2 answers

### Answer by username_1 on 2022-01-02. Score: 3 | Accepted answer

If I follow the code, you are using Kepler's formulae to calculate the distance and speed of the planet, and then increasing the planet's position using speed × time. That is essentially a form of Euler's method for solving numerical differential equations, and as such it may be sensitive to initial conditions and rounding.

You've chosen a small step size of one second. Small steps can reduce one source of error (the essential error in Euler methods of approximation of a continuous process by a discrete one) but can introduce another error: If you you add a very small amount to a large amount, you can lose a lot of significant figure accuracy (eg to 10 sig fig, 12345.00000 + 0.00001234567890 = 12345.00001, you would have lost all but the first significant figure. Now consider what would happen with `n=12345; for i in range 10000 { n+= 1.23456789e-5}` compared with `n+= 1.23456789e-5 * 10000`)

I suspect that the issue here is a matter of programming, rather than with the formulae. Standard good practice is to do all solar system programming not in SI units, but in AU-SolarMass-Years units. These units can help avoid some of the issues with very small values, They make many of the constants, such as $G$ or $\mu$ more manageable (indeed in these units $G=\mu=4\pi^2$)

Moreover I wonder about the general method. You are using Kepler's laws, so why not just use Kepler's law for the orbital period? Or, if you want to investigate orbital dynamics, why not use Newton's law of Gravity?

> **Comments:**
> **username_3**: From what I've read, solar system programming used to use AU as the unit of distance, solar mass as the unit of mass, and a day (not a year) as the unit of time. This practice dates back to Gauss. The current practice uses the AU as a convenience distance whose value is defined as exactly 149597870700 meters rather than a calculated value based on the Gaussian gravitational constant.  

---

### Answer by username_4 on 2022-01-03. Score: 1

Too long for a comment, I refer you back to [@JamesK's answer](https://astronomy.stackexchange.com/a/48037/7982) for the answer.

In Space Exploration SE See:

* [Calculating the planets and moons based on Newtons's gravitational force](https://space.stackexchange.com/q/22915/12102)
* [How to calculate the planets and moons beyond Newtons's gravitational force?](https://space.stackexchange.com/q/23408/12102)
  + Use a reasonable ODE solver; *at least* [RK4](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods) (the classic Runge-Kutta method) or better. Do not use just the [Euler Method](https://en.wikipedia.org/wiki/Euler_method)
* [Two Body Kepler Propogation: Errors from using float / single precision](https://space.stackexchange.com/q/43726/12102)
* [this answer to *Differences between numerical propagators*](https://space.stackexchange.com/a/55906/12102) which links to interesting and helpful resources!

Also see:

* [Semi-implicit Euler method](https://en.wikipedia.org/wiki/Semi-implicit_Euler_method) as mentioned in [@DavidHammen's comment](https://space.stackexchange.com/questions/23408/how-to-calculate-the-planets-and-moons-beyond-newtonss-gravitational-force?noredirect=1&lq=1#comment103549_23409)


---

