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

In this report we will be studying the transit between Venus, Earth, and the Sun. Venus along with Mercury are the only two planets we can observe in transit. The transit of Venus occurs when its revolutions place it in front of the sun between Earth which to us appears as a dot moving across the Sun. To be able to find the differences in planet positions, we use vectors to model planet positions. If we imagine a line connecting two planets in different positions to the sun, the vectors will form an angle. The dot product of these vectors is related to the angle between the vectors. If the cosine of the angle between the planets is $1$, then it means that the angle between the planets is $0$ radians, and the planets are in alignment. The time from when the dot product first equals $1$ to the next occurrence of the dot product equaling $1$ is the transit time. We will first calculate this time for an alignment of Venus and Earth, then expand to incorporate more than two planets.

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$.

# Importing calculation software
import numpy as np
import pandas as pd
import pylab as pl

#Position of a planet in radians after a certain time t
def compute_angle_position(M, m, t, r):
    G = 6.67e-11
    T = (2*np.pi)*np.sqrt((r**3)/(G*(M+m)))
    omega = (2*np.pi)/T
    return(omega*t)

#This cell simply checks if the formula above is correct. For data of earth's orbit around the sun, after 365 days written in seconds we should get a displacement of close to 2pi radians.
theta = compute_angle_position(1.99e30, 5.97e24, 365*24*3600, 1.5e11)
print(theta)

6.2540280526809795

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.

#This will compute the angular positions of Earth and Venus after one Earth year. 
theta_venus = []
theta_earth = []
times = []
for t in range(0, 365):
    times.append(t)
    theta_venus.append(compute_angle_position(1.99e30, 4.88e24, t*24*3600, 1.08e11))
    theta_earth.append(compute_angle_position(1.99e30, 5.97e24, t*24*3600, 1.5e11))

#This cell will create a graph relating the positions of Earth and Venus after one Earth year.
pl.plot(times, theta_venus, label="Venus")
pl.plot(times, theta_earth, label="Earth")
pl.xlabel("Time (Days)")
pl.ylabel("Angular displacement (radians)")
pl.legend()

<matplotlib.legend.Legend at 0x7f619703a7f0>

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.

"""The function below takes the masses of the two planets, the mass of the star,
    the radii of the planets, and the time since the last alignment, and calculates the 
    dot product of the unit radial vectors of the two planets connecting the star.
    
    PARAMETERS:
    ===========
    t = Elapsed time at which the angular position is being computed (days)
    m1 = Mass of the inner planet (default = mass of Venus)
    m2 = Mass of the outer planet (default = mass of Earth)
    m_star = Mass of the star (default = mass of the Sun)
    r1 = Radius of the orbit of the inner planet (default = radius of orbit of Venus)
    r2 = Radius of the orbit of the outer planet (default = radius of orbit of Earth)
    
    RETURNS:
    ========
    dot product of the radial vectors of the two planets with their star.
    """
    
def dot_product_rad_vect(t, m1=4.88e24, m2=5.97e24, m_star=1.99e30, r1=1.08e11, r2=1.5e11):
    twoPi = 2*np.pi
    theta_planet_inner = compute_angle_position(m_star, m1, t*24*3600, r1)
    theta_planet_outer = compute_angle_position(m_star, m2, t*24*3600, r2)

    dot_prod = np.cos(theta_planet_inner)*np.cos(theta_planet_outer) +\
               np.sin(theta_planet_inner)*np.sin(theta_planet_outer)
   
    
    return dot_prod

times = [] # To store the times for graphing
dot_prods = [] # To store the dot products for graphing
for t in range(0, 10*365):     # To graph for ten years, calculated in days
    times.append(t)
    dot_prods.append(dot_product_rad_vect(t))

# This cell will graph time in days vs dot product of radial position vectors
pl.rcParams.update({'font.size': 18})
pl.figure(figsize=(25,6))
pl.plot(times, dot_prods)
pl.xlabel("Time (days)")
pl.ylabel("Dot product of the radial position vectors")
pl.show()

Now we will find the peak points at which the dot product equals or is extremely close to one.

#This code will make sure that only the relative maximum points are taken for each peak.
T = []
for ii in range(1, len(times)-1):
    if (dot_prods[ii-1] < dot_prods[ii]) and (dot_prods[ii] > dot_prods[ii+1]):
        T.append(times[ii])

# This code will add lines at the graph peaks, or points of alignment
pl.rcParams.update({'font.size': 18})
pl.figure(figsize=(25,6))
pl.plot(times, dot_prods)
for thisT in T: # Plotting vertical lines at the points of alignment
    pl.plot([thisT, thisT], [-1, 1], 'r--')
pl.xlabel("Time (days)")
pl.ylabel("Dot product of the radial vectors")
pl.ylim([-1.1, 1.1])
pl.show()

# Display of alignment times
print("Consecutive alignment times for the next ten years after initial alignment:")
for tt in T:
    print("{} years".format(round(tt/365, 2)))

Consecutive alignment times for the next ten years after initial alignment:
1.58 years
3.16 years
4.73 years
6.31 years
7.89 years
9.47 years

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.

"""This function takes given masses and orbital radii of the planets as 
    an array, as well as the mass of the star and the time since the last alignment, and calculates 
    the sum of the dot products of the unit radial vectors connecting the star and one of
    the planets with the other planets.
    
    PARAMETERS:
    ===========
    t = Elapsed time at which the angular position is being computed (days)
    m_array = Array of the masses of all the planets.
    r_array = Array of the radii of the orbits of all the planets
    m_star = Mass of the star (default = mass of the Sun)
    
    RETURNS:
    ========
    Sum of the dot product of the radial vectors.
    """
def multiple_planet_alignment(t, m_array, r_array, m_star=1.99e30):
    twoPi = 2*np.pi
    theta_planet_first = compute_angle_position(m_star, m_array[0], t*24*3600, r_array[0])
    dot_prod = 0
    for ii in range(1, len(m_array)):
        theta_planet = compute_angle_position(m_star, m_array[ii], t*24*3600, r_array[ii])
        dot_prod += np.cos(theta_planet_first)*np.cos(theta_planet) +\
                    np.sin(theta_planet_first)*np.sin(theta_planet)
   
    
    return dot_prod

Before continuing, we will check to see if this function works for the Venus and Earth alignment we calculated above.

m_array = [4.88e24, 5.97e24]
r_array = [1.08e11, 1.5e11]
times_new = [] # To store the time steps for plotting
dot_prods_new = [] # To store the dot product for plotting
for t in range(0, 10*365):
    times_new.append(t)
    dot_prods_new.append(multiple_planet_alignment(t, m_array, r_array))
pl.rcParams.update({'font.size': 18})
pl.figure(figsize=(25,6))
pl.plot(times_new, dot_prods_new)
pl.xlabel("Time (days)")
pl.ylabel("Dot product of the radial vectors")
pl.show()

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 cell will graph the time in days versus sum of the dot products of the radial vectors
m_array = [3.29e23, 4.88e24, 5.97e24, 6.39e23]
r_array = [5.79e10, 1.08e11, 1.5e11, 2.28e11, 2.27e11]
times_triple = []
dot_prods_triple = [] 
for t in range(0, 100*365):
    times_triple.append(t)
    dot_prods_triple.append(multiple_planet_alignment(t, m_array, r_array))
pl.rcParams.update({'font.size': 18})
pl.figure(figsize=(25,6))
pl.plot(times_triple, dot_prods_triple)
pl.xlabel("Time (days)")
pl.ylabel("Sum of dot products of the radial vectors")
pl.show()

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.

# Finding the peaks of the graph above
peaks_triple = []
peak_values = []
for ii in range(1, len(times_triple)-1):
    if (dot_prods_triple[ii-1] < dot_prods_triple[ii]) and (dot_prods_triple[ii] > dot_prods_triple[ii+1]):
        peak_values.append(dot_prods_triple[ii])
        peaks_triple.append(times_triple[ii])

# Making sure that only the peaks are listed in the graph
peaks_triple = np.array(peaks_triple)
peak_values = np.array(peak_values)

# Graphing the times and relative peaks of the function
pl.rcParams.update({'font.size': 18})
pl.figure(figsize=(25,6))
pl.plot(peaks_triple, peak_values)
pl.xlabel("The times at which peaks are seen in days")
pl.ylabel("The peak values ")
pl.show()

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.

cases = peak_values > 2.9 # How many cases of the sum of dot products > 2.9
print(np.sum(cases))

8

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.

cases = peak_values > 2.99 # How many cases of the sum of dot products > 2.99
print(np.sum(cases))

1

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.

# Printing the alignment time(s) in Earth years
peaks_triple[cases]/365
print("Years until next alignment: {} years".format(peaks_triple[cases]/365))

Years until next alignment: [64.65479452] years

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.