3
$\begingroup$

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.

$\endgroup$
7
  • 2
    $\begingroup$ You should go and delete the physics question to avoid cross posting $\endgroup$
    – James K
    Jan 2, 2022 at 16:51
  • 1
    $\begingroup$ @JamesK ok done $\endgroup$ Jan 2, 2022 at 16:53
  • 1
    $\begingroup$ Check out integration techniques / methods and their sensitivity to different error types. Maybe start with something like ias.edu/sites/default/files/sns/Numerical-integration(2).pdf $\endgroup$ Jan 2, 2022 at 18:07
  • $\begingroup$ What you are doing is invalid. You are implicitly assuming a circular orbit in your calculation of angular velocity. $\endgroup$ Jan 3, 2022 at 7:36
  • $\begingroup$ 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. $\endgroup$ Jan 3, 2022 at 7:57

2 Answers 2

3
$\begingroup$

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?

$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$ Jan 3, 2022 at 6:59
1
$\begingroup$

Too long for a comment, I refer you back to @JamesK's answer for the answer.

In Space Exploration SE See:

Also see:

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .