Before we can explore unknown planets and have all kinds of fun in space, let's have a closer look at our solar system. We know it's pretty similar to Earth's solar system, with 8 planets orbiting a star about 3.45 times as massive as Earth's sun. It's safe to assume that the planets orbit the sun, but how can we predict their orbts? Just because a planet is in one place right now, doesn't mean it's in the same place a second from now. Matter of fact something would be really strange if it was.
It's time we explored how planets move and figure out how we can predict their movement.
Philosophers of the past had a firm belief that the circle was the Universe's perfect shape, and therefore, all planetary orbits must be circular. You might already know that they were a little off with this assumption, as planets move in elliptical, not circular, orbits.
Sometime around the early 1600s a man named Tycho Brahe made detailed observations of the planets, only for Johannes Kepler to take all the glory when he used these observations to deduce three laws concerning planetary motion:
- Every planet's orbit is elliptical with the sun at a focus.
- A line joining the sun and a planet sweeps out equal areas over equal time periods.
- The square of a planet's orbital period, P, (in earth years) is proportional to the cube of the semi-major axis, a, of its orbit (in AU). In more mathematical terms, \(P^2=a^3\) .
Here, AU stands for Astronomical Unit, which is the average distance between the sun and the earch, approximately 150 million km (more on this unit later).
Keplers first law simply describes the shape of the orbit, but it was revolutionary at the time. It stems from observation that the planets orbital velocity varies, which ties in with Kepler's second law.
Kepler's second law states that a line segment joining a planet and the Sun sweeps out equal areas during equal time intervals. From this we get an idea of the planets' velocity, as they need higher velocity closer to the sun in order to sweep out the same area as when its further away.
Keplers third law gives us a relation of the size of the orbit compared to the orbital period. In short, the further away a planet is from its star, the more oval the elliptical orbit becomes, or \(p^2 = a^3\). If you like, you can derive these laws with us, if not you can simply believe with all your heart that we are telling the truth.
We can use these laws when we explore how the planets in our system move! We know that the planets move in an elliptical orbit and therefore we can describe this orbit using the expression for an ellipse \(r = \frac{a(1-e^2)}{1+e\cos{f}}\), where, \(e = \sqrt{1+\frac{b^2}{a^2}}\), and \(f = \theta-\omega\). We can plot our orbits using this expression for the position as a function of the angle. Keep in mind that this expression is given in polar coordinates (that depends on the distance \(r\) and the angle \(\theta\)), and as we are working in a coordinate system with an x-axis and a y-axis (Cartesian system), we need to transform from polar coordinates to Cartesian coordinates. We set:
\(x = r \cos \theta \\ y = r \sin \theta\)
We now have everything we need, and can plot our orbits for angles in the interval \([0,2\pi]\).
From this we can see the orbits of our planets around the star, and the size of our solar system. It's a whole 120 AU across! In comparison, the Kuiper belt extends to about 50 AU from the sun, meaning if we measured from one side of the Kuiper belt to another, we'd get 100 AU. You might look at this plot and exclaim "Hey, these are circles you morons," but if we look closely, they do not all have their center in the origin, which means they aren't symmetrican and thus cannot be circles. If we look at our planets eccentricities we see that these values are very small, which will produce almost circular orbits.
This holds for most larger bodies (like planets), while smaller bodies such as asteroids will have more elliptical orbits.
Unfortunately, this general expression for the shape of an ellipse does not let us know where the planets are at any given time. Solar systems are really really really big, so there is some advantage in knowing where planets should be. You know how it is when you look away for a second and suddenly you have no idea where that mischievous planet went. Or even worse, you turn around and suddenly a giant gas planet has caught you in its gravitational orbit and you've no way of escaping! We really need an expression for the position as a function of time.
We need to know the force on each planet if we want to know their acceleration (\(a = \frac{F}{m}\)). If we're realistic here, the planets exert a force on each other, as well as a force on the sun, but we're ignoring this as we don't have the resources needed to store all the extra information this would create. In our simulation, the only force we use is the gravitational force from the sun on the planets.
\(\vec{F} = \gamma\frac{m_sm_p}{r^3}\vec{r}\)
Thus, we have the acceleration, and we already know the velocity and position at our start time. We have our star mass, all of our planets masses, and the distance they are from the star at our start time. Thus, we also know their positions at the start time \(t=0\).. The force is pointed along the vector between the planet and star, and should 'pull' the planets towards the star.
In order to find the position as a function of time, we're going to be using the relation between force, acceleration, velocity, and position. We'll be using a similar method as the last time, only this time we need a method that's even more precise. Introducing: The leap frog method! We mentioned Euler-Cromer in our last post, and it's time to step up the game. Let's dive right into it, study this closely and see what you can make of it before we explain further:
\(t_{i+1}= t_{i} + \Delta t \\ a_{i+1} = - \frac{G m }{|r_i|^3}r_i \\ r_{i+1} = r_i + v_i \Delta t + \frac{1}{2}a_i \Delta t^2 \\ a_{i+1} = - \frac{G m }{|r_{i+1}|^3}r_{i+1} \\ v_{i+1} = v_i + \frac{1}{2}(a_i + a_{i+1})\Delta t \)
As with Euler-Cromer, for every calculation we need to take a small step, \(\Delta t\), forward in time, and every calculation is dependent on the previous one. The integration begins with the acceleration, which is known from the forces acting on the planets. As the velocity is the derivative of position, and acceleration is the second derivative of position, we can integrate to find both of them.
With all numerical calculations, we will always lose a small amount of intormation, as our computers can only handle so many numbers at a time. Unfortunately, that "small" amount can easily turn into quite a large amount when we do enough calculations.
Using the leap frog method, our Frog of Integration is constantly jumping back and forth in our calculations. We notice that the acceleration is involved in every step our integration, which means that we lose less information for every loop, as the same information is being used through every step. Think of it like a whisper game: every time the sentence passes through another person, the information will get slightly tweaked. With the leap frog method, whenever someone whispers something into your ear, the person two seats away from you will also whisper you their version of the sentence. This way you have more information to go on, and will hopefully be able to whisper something sensible to the person after you!
We have many methods for solving second order differential equations, and while leap frog isn't the most precise (although better than Euler-Cromer), the energy is conserved. Without energy conservation, we will either die a very cold or a very hot death, depending on wether our kinetic energy will outweigh the potential (we're flung into space, which is cold), or the other way around (we will fall into the sun, but will at least be nice and warm for a brief moment before we're incinerated). As we want our planets to stay in orbit after thousands of calculations, we will be going forth with our leaping frog <3
You probably noticed we've been measuring distances in Astronomical Units, or AU. When working with astronomical bodies, we often use other units than those we are familiar with in our daily lives. Describing the position of bodies millions of meters apart in meters will quickly lead to some problems. It's hard to keep track of numbers that are 6-7 ciphers long, and missing one ciper in your results can quickly lead you to believe that something is ten (or evern a hundred) times as close as it really is.
That's why we work in astronomical units. 1 AU equals 149 597 870 700 meters, or the average distance between the Earth and its sun. We also use AU/year for our velocities, because astronomical bodies tend to move pretty fast in comparison to what we observe in our daily lives (you are currently orbiting the sun at about 107.000 km/h!).
Now that we have gotten all the practicalities out of the way, let's see what happens when we press 'play' on our simulation!
And here we have our final orbits. "Orbits...?" Yes, the outer orbits in our solar system aren't completed, as we have only simulated over 20 years. Still, they are moving pretty fast, as the outer planet seems to be able to complete its orbit in about 40 years. To put it into perspective, the outermost planet in our solar system, Neptune (sorry, Pluto!), needs 165 years to complete one revolution around the star. This may be because the period of the planet is dependent on the mass of the star, which in our case is about 3.5 times as large as that of that the Earth's sun. A planet in our simulated solar system, of the same mass and at the same distance, would then have a period roughly 1/3 as long as if the planet was in the Earth's solar system.
Most of the orbits fit very well with what we found earlier, but some seem a little off. This could be due to numerical mistakes, but more likely it's because we haven't rotated all the orbits. The angle, f, in our analytical expression is defined as the angle with respect to perhiel (when the planet and star are closest), and so we could rotate some orbits to take this into account. Unfortunately, management is breathing down our necks to stop with the simulations and launch this rocket, and so we couldn't find the time.
What about the planets closer to the sun? We choose to take a closer look at the four innermost planets.
They all have fully drawn circles. As ll of these planets have a mass and radius similar to or smaller than the Earth, it makes sense that they manage to get all the way around their star in 20 years.
We also notice that the circles look eerily perfect. After all, we are doing this on a computer, which usually leads to small errors, and therefore circles that don't completely overlap. So how come that our circles look so good, even though they have been drawn 20 times, or more, on top of each other? The secret is energy conservation! Because of our integration method, the planets don't lose their energy, and therefore stay in orbit.
How can we know for sure that these orbits are correct?
Are our computers lying to us? We decide to check wether they satisfy Kepler's second and third law!
Before we start we make a few simplifications that still should produce satisfuing results.
Kepler's second law describes the area swept over by our planet in a given time.
\(\frac{dA}{dt} = \frac{1}{2}r^2d\theta dt\)
Once again we have an equation we'd 
prefer to not solve analytically, but this time we don't bother with any frogs; we approximate! If we calculate the area for very small pieces of the area and add them together, we can get close to the actual area. You may be familiar with this technique from integrals: we divide the area into small boxes and add them together. In this ellipse, we will be working with triangles instead of boxes.
Let's zoom in on a part of our ellipse and see what happens if we divide our elliptical pie into small slices.
If we let the angle, \(\theta\), be small enough, we can treat the rounded shape of our slice as a straight line. Thus, we have a triangle, and can use trigonometry to find \(\theta\). We consider the small segment \(r d \theta\) to be the height of our triangle.
Now we can easily calculate the area and distance traveled by the planet!
First, we know the area of a triangle is \(A = \frac{1}{2}gh\). An arc length between two points is given 
by \(\frac{2\pi r}{P} \Delta t\), where \(r\) is the radius and \(P\) is the period, or the time of one revolution. This distance will naturally be curved, but we approximate it as straight, and considet it our 'height'. The base of our triangle will be the radius in \(t_0\), or the length of the position vector between the star and our planet. Inserting this into A, we are left with the expresion \(\frac{\pi r^2}{P} \Delta t\). We use Kepler's third law \(P = \sqrt{\frac{4 \pi^2}{G(m_s + m_p)}a^3}\), improved by Newton, to find the period \(P\).
After chopping up the area in fine pieces and adding them together, we check if they match Kepler's second law by calculating two areas, and then subtracting one from the other. If there is no difference, they are the same. Seeing as we have done quite a few simplifications, we allow a small error margin by checking if our calculations satisfy the inequality \(\mid A_2 - A_1 \mid < \epsilon\), where \(\epsilon\) is a small number. If the absolute value of the difference is less than our chosen tolerance, we consider this a satisfying result.
After doing our calculations, we find that our method holds for \(\epsilon = 0.01\), but not for \(\epsilon = 0.001\). On an astronomical scale, we are not too satisfied by this difference in areas, as we are working with large numbers. As we will see shortly, the distance travelled along the orbit during this period of time is 1.28 AU, or 191.485.274 kilometers. That is also the approximated height in our triangle, so an error of 1% will be quite large. One way to fix this would be to slice our orbital pie into even finer pieces, but we leave that up to people with more powerful computers than we have here at our basement space agency.
While we are on the subject of hammering a circle into a triangle, let's see how we can calculate the distance traveled by the planet between two points in time (even if we spoiled the result, the method is half the fun!) We do this by using the height in our triangle. We add all the 'heights' in our triangle between \(t_0\) and \(t_1\), and we have our distance!
We can also calculate the mean velocity during this period of time by summing up every velocity our planet has during this period, and dividing by the number of velocities we have:
\(\langle v \rangle = \frac{1}{N} \sum_{n=1}^{N} v_i\)
We choose to look at our home planet over a period of 0.2 years to find both mean velocity and distance, and we get the following results:
This gives us an orbital speed of 87.890 km/h, compared to the Earths orbital speed of 107.153 km/h. As we are not several magnitudes off, we can conclude that this is a likely mean velocity for our planet in this period. The velocity will vary only slightly if we calculate it for a longer period in time, as our orbits are almost circular. This is because our orbits are close to being perfectly circular, in which case the mean velocity would be constant.
We have confirmed that our planets behave according to Kepler's second law, but what about Kepler's third law:
\(P^2 = a^3\)
Where P is the orbital period of a planet around the sun and a is the semi-major axis. Kepler didn't get it quite right, though. Newton put the last finish on the equation, introducing the gravitational constant, G, and the masses of the sun and the planet:
\(P^2 = \frac{4 \pi^2}{G(m_s + m_p)}a^3\)
This may seem like a dramatic change, but does it really matter…? G is a small number, sure, but compared to the sum of the masses of two large bodies like the sun and a planet, the denominator will still be large. This means that the constant before \(a^3\) will sometimes be small enough to be insignificant. Let's see if our planets behave according to Kepler's third law, alpha- and beta version! The way we choose to do this is by shuffling around the two laws so we get an expression for \(P\). By setting both of them equal, we get
\(\sqrt{\frac{4 \pi^2}{G(m_s + m_p)}a^3} = \sqrt{a^3} \Longleftrightarrow \sqrt{\frac{4 \pi^2}{G(m_s + m_p)}a^3} - \sqrt{a^3} = 0\)
By subtracting one from the other, we can get the difference between them. If that difference is 0, they are equal. As we are doing this numerically, we add a small tolerance to account for numerical errors. We are also taking the absolute value of the difference, as a negative difference will always be less than our tolerance of 0.01.
As we can see, none of our planets hold up. Planet 3 gets pretty close, but 0.7 years off is quite a lot for a planet with a period of 1.5 years— *record scratch* Look at Planet 0's period. We concluded that our home planet and star are pretty close to the solar system's Earth and sun, so how come our period is 38 years? This is probably due to a) a bug fluttering around in our system or b) a sheer lack of competence in two second-year physics students turned rocket scientists. Looking at the data Kepler used vs modern data, we can see that he was closer than what our program predicts. His calculation of the Earths period was off by 0.0064 days, so we have got to give him some credit!
Despite this error, we believe that our orbits hold up, and the mistake lies in this test- the orbits look good and follow Kepler's second law after all.
Our planets move around in what looks to be nice orbits, but we did make a lot of approximations and simplifications in our calculations- maybe you'd like to have a look at a more realistic model? Tune in next time, when we for one post only (!!!), actually take into account that the star, too, is orbiting something :0