Writing code for a solar system

A Solar System in Pieces

To implement the model of orbits I started in a previous post I'll get the most basic parts down and then try to improve the model by making the code work a little better.

Before I describe the solution, it's important to remember that reading is sequential (this post), but building models and programming is most definitely not. It wasn't until I had finished the first coding pass that I came back around and distilled the implementation into three main components:

1. Real world data. Have the model use some approximately accurate data about the orbit of planets.
2. Angular velocity. Obtain or calculate the angular velocity of Mars and Earth.
3. x,y coordinates. Use the angular velocity and trigonometry to determine x and y coordinates of each planet on an (x,y) grid.

Why Angular Velocity?

Not only did it take me some time to formalize the implementation of the conceptual model into three components, it wasn't until after I had programmed the basic model, written some blog posts and made up the images that I actually figured out why an explicit calculation of angular velocity is important in this model.

Even though trigonometry is important to the implementation, Angular Velocity is really the core of this approach to modelling orbits.

It allows us to achieve two important things: Determine relative location as a function of time and closely related, pseudo-realistically model the speed of a planet in its orbit.

Real world data

I'll walk through the basic implementation of the model, found here:

Distance from the Sun

The average distance of Earth and Mars from the Sun becomes the radius in each planet's circular orbit. Earth's distance from the Sun is about 149600000 km, while Mars is about 227900000 km away. Mars is thus 1.5 times further from the Sun in this model and, I assume, real life.

Orbital period

  # Earth
  # Using averages of distance from the sun. 
  earthRadius = 149600000 # km
  earthOrbitalPeriod = 365 * 86400 # seconds in an Earth year

  # Mars
  marsRadius = 227900000 # km
  marsOrbitalPeriod = 687 * 86400 # seconds in a Mars year
                

The orbital period - the time for Earth and Mars to go around the sun - is required for the calculation of angular velocity, and I expressed it in seconds. So, 365 days * seconds in a day == 86400. I suppose any unit of time could be used, and the choice depends on purpose and context.

Choosing the units

Here the choice is mostly about convenience. I guessed seconds would give enough points on the orbit to work with, and could be converted to other units of time pretty easily. I didn't want to work with a unit like years as I'd end up working in 1/365 year-units, or some other fraction. You know, I honestly didn't give it much thought.

Units, details
This is the kind of decision that usually becomes less arbitrary as the importance of the project increases. In the case of units, either the tools (code libraries, math) you work with will have well established conventions, or the conventions will be apparent from the context you're working in. Deviating from convention would require a good reason, not least because it will be more work for you. Take the path of least resistance.

Other data points

  firstQuarterEarthYearInDays = 365 / 4
  halfEarthYearInDays = 365 / 2
  threeQuarterEarthYearInDays = 365 / 4 * 3
  fullEarthYearInDays = 365
                

Pedestrian figures can come in handy and it's often worth writing them down explicitly in the code as variables or constants. In this example, Earth days in a complete orbit (365), and the same for Mars (687) are named with variables.

There are a few other figures here, like the number of days in half or a quarter of an orbit. These look silly at first, but I used them to print out the coordinates of Earth at points in the "orbit" I could verify as correct as a very basic form of testing.

I haven't worried about code quality at this point, but naming "magical" numbers like this can enhance the readability of the code. If something like 365/4 wasn't named, it might not be obvious when read later what that value was about. Giving it a variable name is a form of documentation. I wouldn't take this as a great example of this idea in practice, I'm just getting it all down.

Angular velocity

  # Earth
  # Using averages of distance from the sun. 
  earthOrbitalPeriod = 365 * 86400 # s
  earthAngularVelocity = 2 * math.pi / earthOrbitalPeriod # rad/s

  # Mars
  marsOrbitalPeriod = 687 * 86400
  marsAngularVelocity = 2 * math.pi / marsOrbitalPeriod # rad/s
                

The next step is to determine the angular velocity of Earth and Mars.

The angular velocity of a planet in orbit around a central point is the 360 degrees or 2PI radians of the entire orbit, divided by the time for a full orbit. Consistent with the choice of seconds as the standard unit of time, I expressed the angular velocity per-second.

It's also worth thinking about how averages are affecting the model here. By assuming the orbits are perfect circles with a radius equal to the average distance from the sun at all points on the circle, we don't have to worry about how irregularities are going to change the angular velocity - it will be a constant.

Better questions
The use of angular velocity is interesting in that it's a "derived value," not an observed value from the natural world - like the distance from the Sun. It's also a constant, it doesn't change at any point in the orbit. But if orbits are irregular in real life, does the angular velocity change throughout different parts of the orbit? Is angular velocity even a useful concept in a more realistic model?

The attempt to build a simple model helps us ask more interesting questions about the problem, because better questions come to you as you become more familiar with the language of the problem. In this case, some light searching suggests new avenues of knowledge. Kepler... I've seen that name around.

Angular velocity + trigonometry = x,y coordinates

Here's where we get to the trigonometry. Knowing the angular velocity of each planet and their distance to the sun (radius), we can figure out their positions on an x,y grid using trigonometry.

It took me awhile to figure out how to turn some of this into working code, and then to become comfortable with the ins and outs of how it worked. That to say, I had a solution before I understood it well. As curious as this seems, it's quite common when programming. You'll look at the output and think "that seems right," without fully comprehending what's going on in detail.

It's a regular part of the process when you're learning the problem space.

Orbits, angles and trigonometry

Having derived the angular velocity of each planet, the total angle (in radians) can be calculated for an elapsed period of time through which a planet is moving in its orbit.

To use Earth as an example, its angular velocity in this model is 1.99e-07 radians per second. Take a look at the next image, which shows one quarter of Earth's orbit. Let's figure out the angle, A1, for position 1.

First we need the time it would take for Earth to move through that much orbit. That looks to be about one third, of one quarter, of the total orbit. So, 365 days * 1/4 * 1/3 * 86400 seconds in a day = 2628000 seconds. About 30 days.

To get the angle we take the previous figure, which is the total time the Earth spends moving that far in its orbit, and multiply by the angular velocity of Earth. So 2628000 seconds * 1.99e-07 radians/second = some-angle-in-radians.

Cosine and Sine

Having established how to calculate the angle of the planet as a function of time, notice now that you can draw a right angled triangle between the Sun, the Earth and the x or y axis.

Calculating the x coordinate of a planet

Notice also that the length of the adjacent side of the triangle is the x coordinate of Earth, while the length of the opposite side from angle A is the y coordinate of Earth in its orbit. We need to get these lengths to find Earth's (x,y) coordinates at this point in time on our x,y grid!

Calculating the y coordinate of a planet

I'm not the person to teach trignometry... suffice to say, the basic insight is that the cosine function is the ratio of the adjacent side and hypoteneuse, while the sine function is the ratio of the opposite side and hypoteneuse in a right triangle.

If we know the angle and the radius (the hypoteneuse), we can use these functions to calculate the opposite and adjacent sides of the triangle. These values are our <x and y!

A model in code

  import math

  secondsInADay = 86400

  def get_planet_xy(radius, days, angularVelocity):
    
    x = radius * math.cos(angularVelocity * days * secondsInADay)
    y = radius * math.sin(angularVelocity * days * secondsInADay)
    return x, y
                

The end result of all this is one simple Python function that uses the Python math library implementation of the math.cos() and math.sin() function to return the x and y coordinate of a planet.

This is a general function, and takes as parameters the distance of the planet from the Sun (radius), the number of days elapsed from some start point, and the angular velocity of the planet. If you plug these values, for any planet, into the function it will return the x and y coordinate at that point in time.

Modelling the orbit through time

The following image shows two points on the "orbit" of a planet for which we can calculate its (x,y) coordinates using this method. The next step is to take this basic implementation and extend it so we can learn something of how the orbits of Earth and Mars affect the distance between the two over time.

Making it work

There are a number of different ways this approach could be written in code - all depending on the purpose to which the code is put. It doesn't really matter at the moment because the objective was to get something down that does the basics.

That being said, I don't like this code for a variety of reasons and there are a number of obvious ways it could be improved. The next post will discuss some of the ideas and issues you'll face writing code, and then a few examples of how to improve this implementation.