# SOLAR SYSTEM

## In this tutorial we take a closer look at the solar system using Python and Matplotlib, to create a scale model of the solar system. In this tutorial, we will show the size, position and orbits of each planet.

### Art and Web Design

Contact me for information about rates and availability.

### The Solar System

Since the dawn of mankind we have look towards the stars with awe and wonder. Over the centuries, astronomy-priests have identified and tracked various celestial bodies that surround our world in a better attempt to understand the universe. Exploring our solar system gives us insight into the formation of our own planet and a possible view into the possible future of our planet. Everything from the Sun, asteroids and comets, to dwarf planets, terrestrial planets and gas giants, have been explored and charted.

Before we begin coding our script, we need to collect the necessary planetary data. Most of the following data was obtained using the NASA Planetary Fact Sheet.

In this tutorial, all calculations will be done using the metric system. In some cases data is available in feet or miles (e.g. the Imperial system) and must be converted to the meters or kilometers (the metric system). If you are using the NASA Planetary Fact Sheet, don't forget to convert from diameter to radius for each celestial body by using the following code:

>>> # Calculate the radius by dividing the diameter by 2
>>> radius_planet = diameter_planet / 2

### Universal Constants

The following universal constants are used to perform various scientific calculations :

Speed of light (c) = 3 x 108 m/s
Planck's constant (h) = 6.626 x 10-34 J - s
Gravitational constant (G) = 6.67408 x 10-11 m3 kg3 s-2
Gravitational parameter - Earth (μEarth) = 3.986 x 105 km3/s2

These universal constants can easily be coded in Python:

>>> # Universal Constants
>>> speed_light = 3 * 10**8 # m/s
>>> planck_constant = 6.626 * 10**-34 # J - s
>>> gravitational_constant = 6.67408 * 10**-11 # m^3/kg^3 s^-2
>>> gravitational_parameter_earth = 3.986 * 10**5 # km^3/s^2

### Planetary Data

Planetary statistics can be found in any planetary science book or good physics book, but the internet is also a good place to search for planetary data. NASA offers planetary facts in metric and standard formats. We use the metric system for all scientific computing.

Planetary Fact Sheet
Planetary Body Mass (kg) Radius (km) Orbital Period (days) Rotation Period (hrs) Inclination (deg) Eccentricity Velocity (km/s) Mean Surface Temp (C) Distance from Sun (km)
Mercury 3.30 x 1023 2439 87.98 58.65 0 0.205 47.4 430 580 x 106
Venus 4.87 x 1024 6051 224.7 5832.24 117.4 0.007 35 480 108 x 106
Earth 5.98 x 1024 6378 365.36 24 23.5 0.017 29.8 22 150 x 106
Mars 0.64 x 1023 3393 686.98 24.56 23.98 0.094 24 -23 228 x 106
Jupiter 1.90 x 1027 71492 11.86 9.89 3.1 0.049 13.1 -150 778 x 106
Saturn 5.69 x 1026 60268 29.41 10.24 27 0.057 9.7 -80 1432 x 106
Uranus 8.68 x 1025 25559 84.04 17.25 97.9 0.046 6.8 -215 2871 x 106
Neptune 1.02 x 1026 24764 164.8 16.10 29.6 0.011 5.4 -220 4498 x 106
Sun 1.99 x 1030 696000 *** 24.47 - - 220 5600 -
Moon 7.35 x 1022 1738 27 648 6.7 0.055 1 127 (-173) **

*** Based on a distance of 30,000 light years, and travelling at a speed of 220 km/s, the Sun takes approximately 225 million years to orbit around the Milky Way.
** The Moon's average distance to Earth is 384,400 km (238,900 mi).
- Missing or incomplete data.

We will save the above planetary data as series of lists, multi-dimensional arrays or dictionaries. For our purposes, we will create lists of intergers, strings and lists of variables. These lists can then be used to by our script, printed to the screen, or plotted in bar or scatter charts. We will learn how this data can be visualized using the popular Matplotlib library. To install Matplotlib on Debian/Unix systems, open a terminal window and enter the following command:

\$ sudo apt-get install python-matplotlib

To install Matplotlib on a Windows operating system, open the command prompt and enter the following commands:

\$ python -m pip install -U pip setuptools
\$ python -m pip install matplotlib

Matplotlib will be used to plot various planetary features, such as the average distance of each planet from the Sun, orbital periods and more. In order to use Matplotlib in our script, we must first import it. Begin by opening a python shell. To enter the Python shell in Windows command prompt or Unix terminal, enter the following command:

\$ python
>>>

You can test the Python shell by entering a simple equation. For example, you could enter the following:

>>> 12*12
144

Begin by creating a list of strings, named planets to hold the names of each planet (to include the moon and Sun). A list is an ordered collection of values. In Python, a list is similar to what other languages call an array, but with some added functionality. Lists can contain any type of data, from strings, intergers and floating point numbers. The items that make of a list are simply called items, or elements.

>>> # List of known planets orbiting the Sun
>>> planets = ['Mercury', 'Venus', 'Earth', 'Mars', 'Jupiter', 'Saturn', 'Uranus', 'Neptune', 'Pluto', 'Sun', 'Moon']

We can reference each planetary body by referencing it's position in the list.

>>> print("The planet", planets[0], " is closest to the Sun")
The planet Mercury is closest to the Sun

First, we will identify the planetary variables for each individual planet (such as the mass, radius, density and volume), and then use these individual variables to create lists of the collected variables. Notice that some of these parameters are very small or very large and therefore are presented using the scientific notation. In Python, the ** operator is used to raise the number on the left to the power of the exponent on the right. For example, to raise 5 to the power of 3 (53), that is to multiply 5 by itself 3 times, we use the following expression:

>>> 5**3
125

When listing the following planetary variables, it is important to list each variable in the exact order as in the list of planets above. This allows us to programatically call upon each element in the list by accessing it's order in the list. Recall that elements in a list begin at 0, rather than 1. For example, in order to access the planet Mars from the above list you could use:

>>> planets[3]
'Mars'

We will begin by including the planetary mass of each planet, to include the moon and Sun. Planetary mass is a measure of the mass of each planet and is measured in kilograms (kg). The planetary masses are measured in septilion (1 followed by 24 zeros) kilograms. As you will notice, the planetary masses are listed in scientific notation and must be calculated to get the actual numbers for each planetary mass.

>>> # The mass for each Celestial body (10^24 kg)
>>> mass_mercury = 0.33 * 10**24
>>> mass_venus = 4.87 * 10**24
>>> mass_earth = 5.98 * 10**24
>>> mass_mars = 0.64 * 10**24
>>> mass_jupiter = 1898.19 * 10**24
>>> mass_saturn = 568.34 * 10**24
>>> mass_uranus = 86.8 * 10**24
>>> mass_neptune = 102.41 * 10**24
>>> mass_pluto = 0.01303 * 10**24
>>> mass_sun = 1988500 * 10**24
>>> mass_moon = 0.07 * 10**24

Once these parameters are calculated, they can be added to a list of planetary masses (more specifically, we will create a 2-dimensional array, where each element of the array is a variable pointing to the actual value), so that we can then calculate the gravitational attraction and escape velocity for each planetary body.

>>> # Mass of each planetary body as a list of variables
>>> mass_list = [mass_mercury, mass_veenus, mass_earth, mass_mars, mass_jupiter, mass_saturn, mass_uranus, mass_neptune, mass_pluto, mass_sun, mass_moon]

We could just use the above code, but we prefer to work with clean data. So we can simply print the list of planetary masses, then copy the results and replace the variables in the above list (mass_list) with the actual values.

>>> # Mass of each planetary body as a list of values
>>> mass_list
[3.3e+23, 4.87e+24, 5.98e+24, 6.4e+23, 1.89819e+27, 5.6834e+26, 8.68e+25, 1.0241e+26, 1.303e+22, 1988500000000000000000000000000, 7e+22]

Then copy the results and paste it in place of the list of variables used above.

>>> mass_list = [3.3e+23, 4.87e+24, 5.98e+24, 6.4e+23, 1.89819e+27, 5.6834e+26, 8.68e+25, 1.0241e+26, 1.303e+22, 1988500000000000000000000000000, 7e+22]

Next, we will calculate the radius for each planet, based on it's diameter (NASA Planetary Fact Sheet provides the diameter of each planetary body in kilometers). For our purposes, we are only interested in the radius of each planet, as this will be used when determining the gravitational forces acting on each planet. Also, we will want to convert these numbers from kilometers to meters before using them.

If you are using the Planetary Fact Sheet displayed at the top of the page, the data has been converted for you, however, if you are using the NASA Planetary Fact Sheet, then you may want to use the following code to convert from diameter to radius and from kilometers to meters:

>>> # The radius of each Celestial body (in meters)
>>> radius_mercury = 2439.7 * 10**3
>>> radius_venus = 6051.8 * 10**3
>>> radius_earth = 6378.8 * 10**3
>>> radius_mars = 3396.2 * 10**3
>>> radius_jupiter = 71492 * 10**3
>>> radius_saturn = 60268 * 10**3
>>> radius_uranus = 25559 * 10**3
>>> radius_neptune = 24764 * 10**3
>>> radius_pluto = 1187 * 10**3
>>> radius_sun = 695700 * 10**3
>>> radius_moon = 1738 * 10**3

Again, we want to create a list containing each planetary radii and then create a new list, using the values rather than the variables.

>>> # Radius of each Celestial body as a list
[2439700, 6051800, 6378800, 3396200, 71492000, 60268000, 25559000, 24764000, 1187000, 695700000, 1738000]
>>> radius_list = [2439700, 6051800, 6378800, 3396200, 71492000, 60268000, 25559000, 24764000, 1187000, 695700000, 1738000]

We will also track each planet's polar radius, which is often much smaller than the planet's equatorial radius. We will create the list of polar radii following the same processes described above.

>>> # Polar radius of each Celestial body (in meters)
>>> polar_radius_mercury = 2439.7 * 10**3
>>> polar_radius_venus = 6051.8 * 10**3
>>> polar_radius_earth = 6356.8 * 10**3
>>> polar_radius_mars = 3376.2 * 10**3
>>> polar_radius_jupiter = 66854 * 10**3
>>> polar_radius_saturn = 54364 * 10**3
>>> polar_radius_uranus = 24973 * 10**3
>>> polar_radius_neptune = 24341 * 10**3
>>> polar_radius_pluto = 1187 * 10**3
>>> polar_radius_sun = 695700 * 10**3
>>> polar_radius_moon = 1736 * 10**3
>>> # Polar radius of each Celestial body as a list of variables
>>> # Polar radius of each Celestial body as a list of values
>>> polar_radius_list = [2439700, 6051800, 6356800, 3376200, 66854000, 54364000, 24973000, 24341000, 1187000, 695700000, 1736000]

Now we can easily calculate the gravitational acceleration of each planetary body.

>>> # The gravitational acceleration vector of each Celestial body (in Newtons)
>>> gravity_mercury = G * mass_mercury / radius_mercury**2
>>> gravity_venus = G * mass_venus / radius_venus**2
>>> gravity_earth = G * mass_earth / radius_earth**2
>>> gravity_mars = G * mass_mars / radius_mars**2
>>> gravity_jupiter = G * mass_jupiter / radius_jupiter**2
>>> gravity_saturn = G * mass_saturn / radius_saturn**2
>>> gravity_uranus = G * mass_uranus / radius_uranus**2
>>> gravity_neptune = G * mass_neptune / radius_neptune**2
>>> gravity_pluto = G * mass_pluto / radius_pluto**2
>>> gravity_sun = G * mass_sun / radius_sun**2
>>> gravity_moon = G * mass_moon / radius_moon**2
>>> # Surface gravity of each planetary body as a list
>>> gravity_planets = [3.70, 8.87, 9.81, 3.70, 24.79, 10.4, 8.87, 11.15, 0.62, 274.20, 1.55]

Next we will calculate the volume of each planetary body. A planet's volume tells us how much space it occupies in space. To calulate the volume for each planetary body, we will use the following equation:

V = 4/3 π r3

Before we begin calculating the volume for each planetary body, we need to convert the radius_list from meters to kilometers. We will do this by dividing each element in the radius_list by 103.

[2439.7, 6051.8, 6378.8, 3396.2, 71492.0, 60268.0, 25559.0, 24764.0, 1187.0, 695700.0, 1738.0]

Now we can use this new list to quickly calculate the volume for each planetary body. These measurements closely match the data calculated by NASA for each planet (Planetary Fact Sheets).

>>> volume_list = [4/3*3.14*x**3 for x in radius_list_km]
>>> # Scale the volumes down by dividing each element by 10**10
>>> volume_list_scaled = [x / 10**10 for x in volume_list]
>>> # Volume of each planetary body - scaled (10^10 km^3)
>>> volume_list_scaled
[6.079637193992961, 92.79446788788701, 108.66391334405307, 16.400162740269586, 152982.1375259083, 91649.20951548433, 6990.366701453341, 6358.1500180048215, 0.7001974769893333, 140972497.787736, 2.197949455210667]
>>> volume_list = volume_list_scaled

The average density of each planetary body is measured in kilograms per cubic meter. Dividing mass and volume we gives us the density of each planetary body using the following equation:

ρ = m / V

Before proceeding, we must again scale down the mass of each planetary body. Once the masses are in the correct scale, we can calculate the density of each planet. Again, these calculations are compared to the NASA Planetary Fact Sheets linked above.

>>> # Create a new list of planetary masses by dividing each element in mass_list by 10**19
>>> new_mass_list = [x / 10**19 for x in mass_list]
>>> # Mean Density of each Celestial body (kg/m^3)
>>> mean_density_mercury = (new_mass_list[0] / volume_list[0] # 5428 (kg/m^3)
>>> mean_density_venus = (new_mass_list[1] / volume_list[1] # 5248 (kg/m^3)
>>> mean_density_earth = (new_mass_list[2] / volume_list[2] # 5503 (kg/m^3)
>>> mean_density_mars = (new_mass_list[3] / volume_list[3] # 3902 (kg/m^3)
>>> mean_density_jupiter = (new_mass_list[4] / volume_list[4] # 1241 (kg/m^3)
>>> mean_density_saturn = (new_mass_list[5] / volume_list[5] # 620 (kg/m^3)
>>> mean_density_uranus = (new_mass_list[6] / volume_list[6] # 1242 (kg/m^3)
>>> mean_density_neptune = (new_mass_list[7] / volume_list[7] # 1611 (kg/m^3)
>>> mean_density_pluto = (new_mass_list[8] / volume_list[8] # 1861 (kg/m^3)
>>> mean_density_sun = (new_mass_list[9] / volume_list[9] # 1411 (kg/m^3)
>>> mean_density_moon = (new_mass_list[10] / volume_list[10] # 3185 (kg/m^3)
>>> # Mean density of each Celestial body as a list (kg/m^3)
>>> mean_density_list = (mean_density_mercury, mean_density_venus, mean_density_earth, mean_density_mars, mean_density_jupiter, mean_density_saturn, mean_density_uranus, mean_density_neptune, mean_density_sun, mean_density_moon)

>>> # Planetary cores
>>> core_mercury = "Iron-nickel"
>>> core_venus = "Iron-nickel"
>>> core_earth = "Solid iron-nickel (inner core); Liquid iron-nickel (outer core)"
>>> core_mars = "Iron-iron sulfide"
>>> core_jupiter = "Rocky"
>>> core_saturn = "Rocky"
>>> core_uranus = "Rock-ice"
>>> core_neptune = "Rock-ice"
>>> core_pluto = "Rocky"
>>> core_sun = "Hydrogen"
>>> core_moon = "Solid iron"
>>> # List of planetary cores
>>> planetary_cores = [core_mercury, core_venus, core_earth, core_mars, core_jupiter, core_saturn, core_uranus, core_neptune, core_pluto, core_sun, core_moon]

### Orbital Parameters

In the next section, we will list the semimajor axis, or average distance from the Sun, of each planetary body in kilometers. Keep in mind, all planets follow an elliptical orbit, rather than a perfectly circular one. The point in an orbit when the planet is farthest from the Sun is known as aphelion, likewise, when the planet is closest to the Sun it is known as perihelion. We can also plot these distances, but for our solar system model we will only use the average distances.

>>> # Distance of each planet (Mercury - Pluto) from from Sun in kilometers (as list)
>>> distance_sun = [58000000, 108000000, 150000000, 228000000, 778000000, 1432000000, 2871000000, 4498000000, 5914000000]

At this point we could convert the distances to AU by dividing each distance by 150000000. We can automate this using the following Python code:

>>> # Astronomical unit (1 AU = 150000000 km)
>>> astronomical_unit = 150000000 # 92960000 miles
>>> distance_au = [x / 150000000 for x in distance_sun]
>>> distance_au
[0.38666666666666666, 0.72, 1.0, 1.52, 5.1866666666666665, 9.546666666666667, 19.14, 29.986666666666668, 39.42666666666667]

Other orbital parameters of interest to scientists include the orbital period of each planetary body. The orbital period is the time in Earth days for a planet to orbit the Sun from one vernal equinox to the next. This is also known as the tropical orbit period. We will use the orbital periods provided by NASA (link above).

>>> # Orbit of planet (in Earth days)
>>> orbital_periods = [87.98, 224.7, 365.36, 686.98, 11.86, 29.41, 84.04, 164.8, 248.6]

The speed at which a planetary body travels around the Sun is known as it's orbital velocity, and is measured in kilometers per second. To calculate the orbit of each planet around the Sun, all we need to know is the average distance of the planet from the Sun and the time it takes for the planet to complete one revolution around the Sun. To calculate the planet's orbital velocity we will use the following equation:

v = 2πr / T

>>> # Orbital velocity of each planetary body (in km/s)
>>> orbital_velocity_mercury = (2*3.14*distance_sun[0])/(orbital_periods[0]*24*60*60)
>>> orbital_velocity_mercury
47.91703501637578

We could use the above equation to calculate the velocity of an orbiting satellite or spacecraft, such as the International Space Station. For the International Space Station (ISS), we will collect some more data from NASA. For the radius (r) we will use the radius of the Earth (6378.8 km) plus the altitude of the ISS (408 km). The orbital period of the ISS is 92.49 minutes, which needs to be converted to seconds. Using the above equation we can easily calculate the orbital velocity of the ISS in Python.

>>> T = 92.49 * 60 # 92.46 minutes converted to 5549.4 seconds
>>> r = 6378.8 + 408 # radius of Earth plus altitude of ISS in km
>>> v = (2 * 3.14 * r)/T # velocity is in km/s
>>> v
7.680308501820017

We can confirm with the European Space Agency (ESA) that our calculations are correct. The average orbital velocity of the ISS, according to the ESA, is approximately 7.6 - 7.7 km/s.

### Statistical Analysis

We can use Python and Matplotlib to create visualizations to help us beter understand the relationships between the planets and the Sun. For example, we could show the distances of each planet from the Sun.

>>> import matplotlib.pyplot as plt
>>> plt.figure(1)
>>> plt.barh(planets, distance_sun, color=(0.1, 0.1, 0.1, 0.1), edgecolor='blue')
>>> plt.xlabel('Distance (km)')
>>> plt.ylabel('Planets')
>>> plt.title('Distance of each planet from the Sun')
>>> ax = plt.gca()   # This section turns off scientific notation from changing the xticks
>>> ax.get_xaxis().get_major_formatter().set_scientific(False)
plt.show()

Which will display the following figure.

Looking at the above bar graph allows us to easily visualize the vast distances between planets. Upon closer inspection, we begin to see a pattern to the distances between the planets. It appears that a planet is missing in the spot where the Asteroid Belt is located, between Mars and Jupiter. In the next section of this tutorial we will explore different ways to calculate the weight of a spacecraft or other object on each planetary body.

### Scale Model of the Solar System

The following Python script uses Matplotlib to draw the entire solar system, to scale. I have chosen the date of September 4, 2018 for the approximate placement of the planets in their orbits. I added the orbits for each planet and labels to help identify each planet.

``````
################################################################################
#                                                                              #
#                             Solar System Model                               #
#                        Written by Johnathan Nicolosi                         #
#                         Created September 04, 2018                           #
#                                                                              #
################################################################################

import matplotlib.pyplot as plt

radius_list = [2439700, 6051800, 6378800, 3396200, 71492000, 60268000, 25559000, 24764000, 695700000]

distance_meters = [58000000000, 108000000000, 150000000000, 228000000000, 778000000000, 1432000000000, 2871000000000,
4498000000000]

fig = plt.figure('Solar System')

plt.title('Solar System - 09/04/2018')
plt.axes()

# Sun - host star

# Mercury - orbital path
orbit_mercury = plt.Circle((0, 0), radius=distance_meters[0], fc='none', ec='black')
# Mercury - planet
# Mercury - label
mercury_label = ax.annotate('Mercury', xy=(41031000000, 41007000000), xytext=(41031000000, 41007000000),
horizontalalignment='left')

# Venus - orbital path
orbit_venus = plt.Circle((0, 0), radius=distance_meters[1], fc='none', ec='black')
# Venus - planet
# Venus - label
venus_label = ax.annotate('Venus', xy=(28820000000, -104100000000), xytext=(28820000000, -104100000000),
horizontalalignment='left')

# Earth - orbital path
orbit_earth = plt.Circle((0, 0), radius=distance_meters[2], fc='none', ec='black')
# Earth - planet
# Earth - label
earth_label = ax.annotate('Earth', xy=(140770000000, -51830000000), xytext=(140770000000, -51830000000),
horizontalalignment='left')

# Mars - orbital path
# Mars - planet
# Mars - label
mars_label = ax.annotate('Mars', xy=(165650000000, -156680000000), xytext=(165650000000, -156680000000),
horizontalalignment='left')

# Jupiter - orbital path
orbit_jupiter = plt.Circle((0, 0), radius=distance_meters[4], fc='none', ec='black')
# Jupiter - planet
# Jupiter - label
jupiter_label = ax.annotate('Jupiter', xy=(-447700000000, -636057000000), xytext=(-447700000000, -636057000000),
horizontalalignment='left')

# Saturn - orbital path
orbit_saturn = plt.Circle((0, 0), radius=distance_meters[5], fc='none', ec='black')
# Saturn - planet
# Saturn - label
saturn_label = ax.annotate('Saturn', xy=(150400000000, -1424150000000), xytext=(150400000000, -1424150000000),
horizontalalignment='left')

# Uranus - orbital path
orbit_uranus = plt.Circle((0, 0), radius=distance_meters[6], fc='none', ec='black')
# Uranus - planet
# Uranus - label
uranus_label = ax.annotate('Uranus', xy=(2488900000000, 1431300000000), xytext=(2488900000000, 1431300000000),
horizontalalignment='left')

# Neptune - orbital path
orbit_neptune = plt.Circle((0, 0), radius=distance_meters[7], fc='none', ec='black')
# Neptune - planet
# Neptune - label
neptune_label = ax.annotate('Neptune', xy=(4430200000000, -778050000000), xytext=(4430170000000, -778040000000),
horizontalalignment='right')

plt.ylabel('Distance in meters')
plt.xlabel('Distance in meters')

plt.axis('scaled')
plt.show()
``````

When you run the above code it will create the following diagram.

Zooming in reveals the locations of each of the planets within their orbit. You could zoom in on each planet, including the sun, but this will reveal very little information.

Johnathan Nicolosi - 20 Sep 2018