Gravitation 3 Project (Group 7)

Katerina Petroff

petroffk1@montclair.edu
Group Members: Katerina Petroff, Sydney Winter, Victoria Villanueva, Delvis Burgos, Jett Coronado

Objective

The purpose of this lab is to model simplified orbits of planets using code and to calculate transit times for these planets. We will use data from our solar system, but data from other celestial bodies can be used.

Theory/Introduction

Model of Orbit

To perform calculations for this report, we need to establish the model with which we will predict the movement of the planet or planets. We will assume that the orbit of each planet is perfectly circular and horizontal, with a massive star at the center of the system and constant angular velocity for each planet. For calculation purposes, we will assume that the system of the planets and the star is isolated. We will complete our calculations by modeling the planets and the center star as particles of a certain large mass, with the planets in orbit around the star.

Kepler's Laws

In the early 1600s Johannes Kepler, a German astronomer, published what are now kown as Kepler's laws of planetary motion. These laws describe the movement of orbiting planets in our solar system, and hold to this day. The three laws go as such:
1. The orbit of a planet follows an ellipse pattern with the Sun at one focus
2. A line joining the Sun and an orbiting planet will cover equal areas with equal intervals of time
3. The square of a planet's orbital period is proportional to the cube of the length of the orbit's semimajor  axis
We are assuming that the planets follow a circular orbit, or an ellipse with two foci in the same place and equal semi-major and semi-minor axis lengths. Even in a circular orbit, if the planet is orbiting at a constant speed the line defined in law 2 will still cover equal areas with equal time intervals. With the two first laws upheld, the third law is of most interest to this project. This law will allow us to calculate a planet's period, and therefore, its angular position after a certain time. Assuming a circular orbit, to derive the equation we can set centripetal force equal to gravitational force from Newton's law of gravitation. The result is this equation:

$mr\omega\ = \frac{GmM}{r^2}$

In circular motion, angular velocity in radians per second is designated as

$\omega\ = \frac{2\pi}{T}$

By substituting the second equation into the first, rearranging the variables leads us to:

$T = 2\pi\sqrt{\frac{r^3}{GM}}$

This equation only represents a planet with a circular orbit, and it does not take into account the mass of the planet. By replacing $M$ with $M + m$, we get the final equation for T, the orbital period in seconds.

$T = 2\pi\sqrt{\frac{r^3}{G(M + m)}}$

The variables are:
1. T = planet orbital period, or time for one full revolution around the star, in seconds
2. r = radius of the planet's orbit from the star (center of mass from center of mass)
3. G = gravitational constant, 6.67*10^-11
4. M = mass of the center star being orbited
5. m = mass of the planet orbiting around the center star
For this project, we will model our solar system using the sun as the center star and various planets in orbit. To find the angular position of a planet after a certain time in seconds $t$, we must find the angular velocity.

$\omega = \frac{2\pi}{T}$

$\theta = \omega t$

Transits

A transit occurs when an inner planet's position lines up with the orbit of an outer planet. When a planet is in transit it appears to the observer that the mass is moving across the face of a planet. There are two types of transits, one where the mass is smaller than the planet or when the mass is bigger than the planet, better known as occultations. From the point of view of the outer planet, the inner planet will be seen as a small dark spot slowly traveling across the sun. At the exact moment of transit when a straight line can be drawn throught the sun, the inner planet, and the outer planet, the angular positions of the inner planet and outer planet should be exactly equal. As the planets continue on their orbits, the faster moving planet will draw away from the slower planet, until eventually they line up in transit again. This time period between transits is our target. For more than two planets, the same condition of lining up holds. Since we are modeling the orbits as flat, circular orbits, our calculated transit times will be different from real transit times. This is because in the real solar system, orbits are somewhat elliptical and tilted at angles to each other. It would take much longer for a real transit to occur due to the fact that the planets must not only align in angular position, but they must also match each other at orbital nodes and intersect paths on a tilted axis from the sun. That being said, we can calculate the alignment time for a simplified model of the solar system using code, although it might be somewhat off from the true time between transits.
The main condition for alignment that we will use is that the angular positions of the planets must be equal. Although the equation for angular position above will work for a certain time, it will give the cumulative distance the planet has orbited, not the position. Therefore, faster planets will never match up again with slower planets because they will have orbited more, and a transit will never occur according to the computer. To fix this issue, we initially computed the net position of each planet, removing the revolutions the planet has already completed. We did this by subtracting the integer of the revolutions made for each planet from its total number of revolutions. While this gives us the net difference in position between planets, the function ends up jumping due to the jump from $2/pi$ radians to $0$ for each revolution. To solve this issue, we instead used a separate method for calculating difference in position using cyclic sine and cosine functions. The method is explained below. This will allow us to find the true difference in angular position between planets, despite one planet having covered more distance than the others.

Problem

The goal of this project is to model planetary orbits through code using Kepler's third law. We will develop a function to find the time between transits or alignments of four planets in our solar system, Mercury, Venus, Earth, and Mars. We will first develop a function to calculate the alignment times of Venus and Earth, then expand the function to include more than two planets. The function will make use of radial unit vectors extending from the sun to the respective planets. Technically, if the planets are aligned in one position during transit, the angle between each planet's vector should be zero degrees. Using this, we know that the cosine of the angle between each should be one, and therefore the sum of the radial unit vectors should be a whole number. The number depends on the number of planets, as will be shown below. The algorithms we will use will include graphing tools, vector calculations, and defining functions. There will be notes in the cells of code to explain the exact use of each cell.

Code Conditions

The explanation above was for two planets. Each time a new planet is added, the transit time must take into account another body. Theoretically, during a perfect transit the dot products of the position vectors of the planets should all equal $1$. Therefore, if the separately calculated angles between the inner planet and all of the outer planets of the system all equal $1$, then the angles between the planets are all $0$ radians and the planits are in complete alignment or transit. This means that the sum of the dot products of these vectors (one between each outside planet and the inside planed) should equal the number of outside planets, which would be three in a four planet system. This value is what we will be searching for. When the function hits a maximum at or very near to three, the time at which this event occurs will be the time of transit. We will calculate this in days from the original transit. Our function will calculate the positions of the planets each day from the original transit, so the result might be slightly lower than three, even on the correct day. We will find the times of these occurences above a certain threshold slightly smaller than three to get the alignment times. The time between transits for four planets will be larger than three or two planets, as there are more planets with difference angular velocities to account for.

Base Formulas

These cells will define the base formulas and constants we will be using such as orbital period $T$, angular momentum $\omega$, and angular position $\theta$.

Planet Misalignment

If we do not take into account the full revolutions the planets undergo after intial alignment, the planets will never line up due to their different angular velocities, as demonstrated below. We are going to begin with calculating the transit of Venus from Earth.
While this shows the difference in angular distance traveled between the planets, it does not show the real angle between the planets at a certain time.

Alignment of Venus and Earth

As explained in the code section of the theory, we will find the angle between the position vectors of the planets. The vector functions when written out look like this:

$\hat{r}_v = \cos\theta_{\rm inner} \hat{i} + \sin\theta_{\rm inner} \hat{j}$

$\hat{r}_e = \cos\theta_{\rm outer} \hat{i} + \sin\theta_{\rm outer} \hat{j}$

$\hat{r}_{\rm inner} . \hat{r}_{\rm outer} = \cos\theta_{\rm inner} \cos\theta_{\rm outer} + \sin\theta_{\rm inner} \sin\theta_{\rm outer}$

where $\hat{r}_v$ is the position vector of Venus and $\hat{r}_e$ is the position vector of Earth. $\theta_{\rm inner}$ is the angular displacement of the inner planet Venus at a certain time and $\theta_{\rm outer}$ is the angular displacement of the outer planet, or Earth in this case.
Now we will find the peak points at which the dot product equals or is extremely close to one.

Extending code to include four planets

The code above only included two planets in orbit, Venus and Earth. We will now extend this code to take into account more than two planets. An example using four planets of this solar system, Mercury, Venus, Earth, and Mars, will be calculated. We will find the dot product of of the unit radius vector of the innermost planet with the other planets. During planetary alignment, all of the dot products should equal one.
Before continuing, we will check to see if this function works for the Venus and Earth alignment we calculated above.
This graph is exactly the same as the graph created above. So, we will now extend this system to include all four of the mentioned planets. Orbital radii and planet masses were found online.
This system has run for 100 years. There are many relative maximums, but we are only interested in those that are at or near 3. These peaks are the instances of alignment. Since our step size is one day, a peak that indicates alignment might not necessarily reach exatly 3 as it might fall between two days, so we are searching for peaks that rise above a certain threshold. These occurences will give the times of the next alignments.
The graph above shows the instances of relative maximums. Most are far below 3, so we will set a certain threshold above which we will consider the planets in alignment.
The 8 means that within this 100 year period after the first alignment, there will be eight more alignments. An alignment between two planets is rare in itself, between four planets even more so. Some relative maximums that are not true alignments might have crossed the $2.9$ threshold. Below is the same function if we raise the threshold to 2.99.
One alignment within the next 100 years after an initial alignment seems plausible for four planets orbiting at various speeds and radii from the sun. Below the code will show the time for the next alignment in years after the initial alignment.
Therefore, in reference to our model, a strong alignment between Mercury, Venus, Earth, and Mars will occur about 64.7 years after initial alignment. This time is much less than the time it would take for a real life transit due to tilted orbits and elliptical paths, but it works for our model.