Trouble Simulating Earth’s Elliptical Orbit Motion in Python

by

I’m attempting to simulate the motion of the Earth around the Sun in an elliptical orbit using Python. I’ve tried using basic trigonometry to calculate velocities and update positions, but I’m encountering issues with the Earth’s velocities not changing as expected along the orbit.

Here’s a simplified version of my code:

import math

earthposx = 0
earthposy = 0
sunposx = 10
sunposy = 10

velocity = 1

while True:
    dx = earthposx - sunposx
    dy = earthposy - sunposy
    
    angle = math.atan2(dy, dx)

    velx = velocity * math.cos(angle)
    vely = velocity * math.sin(angle)

    earthposx += velx
    earthposy += vely

    print(velx, vely)

However, the velocities ((velx) and (vely)) are not changing, resulting in the Earth not moving along an elliptical path as expected.

I’m seeking guidance on how to properly simulate the Earth’s elliptical orbit considering changing velocities at different points along the orbit. Should I be considering gravitational forces, numerical methods, or are there other approaches I should explore for a more accurate simulation?

Any insights or code examples on how to accurately simulate the Earth’s elliptical orbit motion in Python would be greatly appreciated. Thanks!

  • The angle you are calculating and then the velx from that are resulting in velx and vely being the same negative number. when you add that to the earthposx/y the angle with the sun in your coordinate system is the same angle. I suggest rechecking which functions you should be using to calculate those items.

    – 

  • Your starting velocity in the x and y direction is 1 which means you are moving away from your ‘sun’. might make more sense to have that separate value for x and y and have them be ‘different’ so its not move away or toward the sun in a straight line. also, you might need to do something like: velx+=velx * math.cos(angle) so that velx changes over each iteration as a function of the cosine and its current value. just some food for thought. its been a LONG time since I thought seriously about trig functions 🙂

    – 

  • @LhasaDad good guesses but the OP’s initial angle is just as equally valid as any other angle, and an initial velocity of 1 is just as valid as any other initial velocity (though this velocity is merely a magnitude and says nothing of direction). Aside from the OP’s velocity being the radius of their circle for some reason, the real issue is that nothing changes between loop iterations. The argument in the trig functions must change to get “motion”. See my answer below for a solution (no velocities needed). If you want velocity, simply take the derivative of position with respect to time.

    – 




Two major issues:

  1. If you want motion, you need to propagate time in your equations
    (i.e., the argument inside of sine and cosine needs to change)
  2. If you want elliptical motion, you need an equation for an
    ellipse not a circle.

For issue 1:

It looks like you are expecting your while loop to be the time propagation, but you need some changing variable within the sine and cosine functions. As you have it now, nothing is changing within these functions and so time is not propagating (nothing is moving). Simply add a time variable that increments every loop iteration and put it in the argument of these functions.

For issue 2:

The equation of a circle (in terms of sine and cosine) is really just a special case of an ellipse (and an ellipse is really just a special case of a conic section). An ellipse (in terms of sine and cosine) is:

x = a * cos(t)
y = b * sin(t)

Where a != b, and 0 < (1 - b**2/a**2)**0.5 < 1. A circle is when a = b and is the radius of the circle. Your code shows the velocity as the circle’s radius for some reason.

The solution:

import matplotlib.pyplot as plt
import math


earth_eccentricity = 0.017

a = 1 # astronomical unit
b = (1 - earth_eccentricity**2)**0.5 # astronomical unit
t = 0 # year

x_track = []
y_track = []

while t <= 1:
    x = a * math.cos(2*math.pi * t)
    y = b * math.sin(2*math.pi * t)
    
    t += 0.001
    
    x_track.append(x)
    y_track.append(y)
    
plt.plot(x_track, y_track)
plt.show()

Which is for Earth’s present value of eccentricity. This produces the following plot:

enter image description here

Also, note that planetary eccentricity is usually very small (often much smaller than the human eye can perceive). When plotted with an aspect ratio of 1, you get the following:

enter image description here

Leave a Comment