3.3: Two Layer Planet Structure Jupyter Notebook

Instructions

• To use this page, read the text between each of the python code windows, then press RUN to execute the code in the box.
• It is necessary to run each of the boxes in order.
• Don't click restart in the cells unless you come back to the top of the page and start over.
• If you modify the code, this modification will not remain after you leave or refresh this page.
• You must run simulations/cells that require user input to completion. Make sure you see the Message informing you a cell is complete before continuing to the next cell. Here's why you must do this:
  • These cells may not be re-run until the simulation is completed, otherwise no output will be generated.
  • All following cells will not run until the user input cell is completed.
• For Dropdown Menus: You do not need to rerun the cells with the menu to change your selection. However, you must rerun all of the following cells to implement the change.

Be patient, some times it can take 1-2 minutes for the juypter kernel to start.

An interactive model to explore the layered structure of far-away planetary interiors from the limited information available to us, assuming there are two distinct layers.

Key Questions: Consider these as you work your way through this page.

Initial Reference Model

• How are the initial density and core size related? For example, does a smaller starting density produce a large or small core radius when compared to a larger starting density?

Density Variation Example

• Find a density value that creates a non-physical (NaN: not-a-number) density for the core of your chosen planet? Why does this happen?
• Find a density value where the model shows a core that is larger than the planet itself? Why does this happen?

Radius Variation Example

• Find a radius value that creates a non-physical (NaN) density for the core of your chosen planet? Why does this happen?
• Find a radius value where the model shows a core is larger than the planet itself? Why does this happen?

Mass Variation Example

• Is there an upper bound on what the mass can be in the model? That is, is there a high value that breaks the model like the previous examples; ignore whether it is realistic or not.
• Find a mass value that creates a non-physical (NaN) density for the core of your chosen planet? Why does this happen?
• Find a mass value where the model shows a core that is larger than the planet itself? Why does this happen?

First let's call the necessary functions, and define a dictionary called planets holding the relevant parameters of different planets.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets

planets = {'Jupiter':{'MOI':0.254,'Mass':1898.19e24,'Radius':69911000,'Surface Density':0.16,'Compositional Density':600,'type':'metallic hydrogen'},
           'Neptune':{'MOI':0.26,'Mass':102.413e24,'Radius':24622000,'Surface Density':0.45,'Compositional Density':1000,'type':'ice'},
           'Earth':{'Mass':5.9723e24,'Radius':6371000,'MOI':0.3308,'type':'Olivine and Peridotite Upper Mantle','Surface Density':2850,
           'Compositional Density':4000}}

print('The packages have been imported and the dictionary defined; move to the next cell')

Reference Planet

Now to start your own model. Select a planet from the prompt once you run this cell.

style = {'description_width': 'initial'}
layout = {'width': '400px'}
planet = widgets.Dropdown(
    options=['Jupiter', 'Earth', 'Neptune'],
    value='Earth',
    description='Select a planet for your model:',
    disabled=False,
    style = style,
    layout = layout
)


display(planet)

Now, choose which density type you would like to use. The surface density can be measured, and the compositional density is the density of the assumed material that the mantle of the planet is made out of.

planet_choice = planet.value
print('The Surface Density of', planet_choice, 'is', planets[planet_choice]['Surface Density'], 'kg/m^3, and the compositional density is', planets[planet_choice]['Compositional Density'],'kg/m^3')
density = widgets.Dropdown(
    options=['Surface', 'Compositional', 'Other'],
    value='Surface',
    description='Select a density for your model:',
    disabled=False,
    style = style,
    layout = layout
)

display(density)

This cell will gather all the required parameters for your choices, and prompt you for your density value if you selected 'other' in the previous cell.

density_choice = density.value

if density_choice.lower() == 'other':
    print('Please enter your desired density in kg/m^3, for now stay between the Surface and Compositional Densities: ')
    rho = float(input()) 
    while rho <= planets[planet_choice]['Surface Density'] or rho >= planets[planet_choice]['Compositional Density']:
        print('Please enter a density between the Surface and Compositonal Densities')
        rho = float(input()) 

print('Here are your selected parameters:')
print('Planet:',planet_choice)
if density_choice.lower() == 'surface':
    print('Surface Density of',planet_choice,'is',planets[planet_choice]['Surface Density'],'kg/m^3')
    rho = planets[planet_choice]['Surface Density']
    if planet_choice == 'Earth':
        print('The Surface Density of Earth is taken to be the Crustal Density')
elif density_choice.lower() == 'other':
    print('Your Entered Density of',planet_choice,'is',rho,'kg/m^3')
else:
    print('Compositional Density of',planet_choice,'is',planets[planet_choice]['Compositional Density'],'kg/m^3 from',planets[planet_choice]['type'])
    rho = planets[planet_choice]['Compositional Density']
    
M = planets[planet_choice]['Mass']
R = planets[planet_choice]['Radius']
I = planets[planet_choice]['MOI']

print('The Mass of',planet_choice,'is',M,'kg')
print('The Radius of',planet_choice,'is',R/1000,'km')
print('The Moment of Inertia Factor of',planet_choice,'is',I)

print('Move on to the next cell.')

This cell defines the parameters beta, a, and b as defined in 3.2: Layered Structure of a Planet. beta depends on the planet's mass and radius. a and b depend on beta, alpha (which is 5/2 times the Moment of Inertia), and the density rho. a is the fraction between the core and total radius, and b is the proportionality constant to calculate the core density from the original density. Therefore, once you calculate \(a\), you can get the core radius by multiplying \(a\) by the total radius (however, for simplicity this model only plots a and not a * R). The density of the core is found by multiplying b by the original density, which is done in the plotting function below.

def beta(M,R):
    denominator = (4 * np.pi * R**3)
    return ((3*M)/denominator)

def a_and_b(beta,rho,alpha):
    w = beta - rho
    q = alpha * beta
    a = np.sqrt((q-rho)/w)
    b = (w / ((a**3) * rho)) + 1
    return a,b

def plotting(M,R,rho,I):
    beta1 = beta(M,R)
    a,b = a_and_b(beta1,rho,I*(5/2))
    fig, ax = plt.subplots(figsize = (7,7))
    ax.add_patch(Circle((0,0),R/R,alpha = 0.1,ec='black',fc='blue'))
    ax.add_patch(Circle((0,0),R/R,fill=False,ec='black'))
    ax.add_patch(Circle((0,0),a,fill=False,ec='black'))
    ax.add_patch(Circle((0,0),a,fc='red',color = 'red',alpha=0.2,ec='black'))
    ax.annotate('$p_1$ = '+'{:.2f}'.format(rho)+' $ kg/m^3 $',(0.05,0.97-((1-a)/2)))
    ax.annotate('$p_{core}$ = '+'{:.0f}'.format(b*rho)+' $ kg/m^3 $',(0.05,0.05))
    ax.set_ylabel('Fraction of total radius: '+ str(R/1000) +' km')
    ax.set_xlabel('Fraction of total radius: '+ str(R/1000) +' km')
    plt.show()
    
print('Here is your Reference Model of',planet_choice,'using its',density_choice,'density:')
plotting(M,R,rho,I)
print('The model is complete, move to the next cell to modify it')

Density Variation Example

Take note of the appearance of the above model. This is generally what the model looks like when it is physically possible for the input parameters to yield a solution. That's why you were limited in your choice for the range for densities if you chose the other option above. In the next series of examples, we are going to consider how the initial mass, radius, and density affect the model, and if certain values of these parameters will create an impossible planet. First, let's change the initial density of the planet. See how changing the density changes the model, and why certain density choices are impossible. You can tell if a choice is impossible if the graph does not look nice and neat like the previous one; you may also encounter non-physical values in the form of nan values (nan/NaN stands for not-a-number) or the radius of the core could end up being larger than the radius of the planet itself.

print('Enter your new surface density value, it must be greater than 0 and less than 5600 in kg/m^3:')
rho = float(input())
while rho <= 0 or rho > 5600:
    print('Invalid Density: it must be greater than 0 or less than or equal to 5600 kg/m^3')
    rho = float(input())

handle = display(None, display_id=True)
fig, ax = plt.subplots(figsize = (7,7))
ax.add_patch(Circle((0,0),R/R,alpha = 0.1,ec='black',fc='blue'))
ax.add_patch(Circle((0,0),R/R,fill=False,ec='black'))
while rho != 0:
    
    print('Here is your new model of',planet_choice,'using a density of',rho,'kg/m^3')
    beta1 = beta(M,R)
    a,b = a_and_b(beta1,rho,I*(5/2))
    
    
    
    patch1 = ax.add_patch(Circle((0,0),a,fill=False,ec='black'))
    patch2 = ax.add_patch(Circle((0,0),a,fc='red',color = 'red',alpha=0.2,ec='black'))
    annot = ax.annotate('$p_1$ = '+'{:.2f}'.format(rho)+' $ kg/m^3 $',(0.05,0.97-((1-a)/2)))
    annot2 = ax.annotate('$p_{core}$ = '+'{:.0f}'.format(b*rho)+' $ kg/m^3 $',(0.05,0.05))
    ax.set_ylabel('Fraction of total radius: '+ str(R/1000) +' km')
    ax.set_xlabel('Fraction of total radius: '+ str(R/1000) +' km')
    handle.update(fig)
    print('Enter 0 to exit the simulation when ready')
    print('Enter your new surface density value, it must be greater than 0 or less than 5600 kg/m^3:')
    rho = float(input())
    while rho < 0 or rho > 5600:
        print('Invalid Density: it must be greater than 0 or less than 5600 kg/m^3')
        rho = float(input())
    patch1.remove()
    patch2.remove()
    annot.remove()
    annot2.remove()
plt.close()
print('The simulation is complete. Consider the key questions.')

Radius Variation Example

Now let's set the density to a constant value, we will use the compositional density of the planet you selected. This time, let's change the radius and see what happens. Like before, take note of when and why the model gives a nan or a core radius larger than the planet itself. Note that the given range of radii are given to accommodate all three planets, it may be easier to find the nan and large cores near the radius of the initial planet.

rho = planets[planet_choice]['Compositional Density'] # Going to hold rho constant at the compositional density for this example
print('The original radius of your planet is',planets[planet_choice]['Radius']/1000,'km')
print('Enter your new radius value in km, it must be greater than 0 or less than 100,000 km:')
R = float(input())
while R <= 0 or R>100000:
    print('Invalid Radius: it must be greater than 0 or less than 100,000 km')
    R = float(input())
R *= 1000 #convert to m
handle = display(None, display_id=True)
fig, ax = plt.subplots(figsize = (7,7))
ax.add_patch(Circle((0,0),R/R,alpha = 0.1,ec='black',fc='blue'))
ax.add_patch(Circle((0,0),R/R,fill=False,ec='black'))
while R != 0:
    
    print('Here is your new model of',planet_choice,'using a radius of',R,'m')
    beta1 = beta(M,R)
    a,b = a_and_b(beta1,rho,I*(5/2))
    
    
    
    patch1 = ax.add_patch(Circle((0,0),a,fill=False,ec='black'))
    patch2 = ax.add_patch(Circle((0,0),a,fc='red',color = 'red',alpha=0.2,ec='black'))
    annot = ax.annotate('$p_1$ = '+'{:.2f}'.format(rho)+' $ kg/m^3 $',(0.05,0.97-((1-a)/2)))
    annot2 = ax.annotate('$p_{core}$ = '+'{:.0f}'.format(b*rho)+' $ kg/m^3 $',(0.05,0.05))
    ax.set_ylabel('Fraction of total radius: '+ str(R/1000) +' km')
    ax.set_xlabel('Fraction of total radius: '+ str(R/1000) +' km')
    handle.update(fig)
    print('Enter 0 to exit the simulation when ready')
    print('Enter your new radius value in km, it must be greater than 0 or less than 100,000 km:')
    R = float(input())
    while R < 0 or R>100000:
        print('Invalid Radius: it must be greater than 0 or less than 100,000 km')
        R = float(input())
    R *= 1000 #convert to m
    patch1.remove()
    patch2.remove()
    annot.remove()
    annot2.remove()
plt.close()
print('The simulation is complete. Consider the key questions.')

Mass Variation Example

Finally, lets set the radius to its original length and change the mass of your planet. Like before, take note of when and why the model gives a nan or a core radius larger than the planet itself.

R = planets[planet_choice]['Radius'] # Going to hold radius constant for this example
print('The original mass of your planet is',planets[planet_choice]['Mass'])
print('Enter your new mass value, it must be greater than 0 and in kg:')
print('The number you enter will be multiplied by 10^24')
M = float(input())
while M <= 0:
    print('Invalid Mass: it must be greater than 0')
    M = float(input())
M *= 10**24
handle = display(None, display_id=True)
fig, ax = plt.subplots(figsize = (7,7))
ax.add_patch(Circle((0,0),R/R,alpha = 0.1,ec='black',fc='blue'))
ax.add_patch(Circle((0,0),R/R,fill=False,ec='black'))
while M != 0:
    
    print('Here is your new model of',planet_choice,'using a mass of','{:.2e}'.format(M),'kg')
    beta1 = beta(M,R)
    a,b = a_and_b(beta1,rho,I*(5/2))
    
    
    
    patch1 = ax.add_patch(Circle((0,0),a,fill=False,ec='black'))
    patch2 = ax.add_patch(Circle((0,0),a,fc='red',color = 'red',alpha=0.2,ec='black'))
    annot = ax.annotate('$p_1$ = '+'{:.2f}'.format(rho)+' $ kg/m^3 $',(0.05,0.97-((1-a)/2)))
    annot2 = ax.annotate('$p_{core}$ = '+'{:.0f}'.format(b*rho)+' $ kg/m^3 $',(0.05,0.05))
    ax.set_ylabel('Fraction of total radius: '+ str(R/1000) +' km')
    ax.set_xlabel('Fraction of total radius: '+ str(R/1000) +' km')
    handle.update(fig)
    print('Enter 0 to exit the simulation when ready')
    print('Enter your new mass value, it must be greater than 0 and in kg:')
    M = float(input())
    while M < 0:
        print('Invalid Mass: it must be greater than 0')
        M = float(input())
    M *= 10**24
    patch1.remove()
    patch2.remove()
    annot.remove()
    annot2.remove()
plt.close()
print('The simulation is complete. Consider the key questions.')

Key Questions: Consider these as you work your way through this page.

Initial Reference Model

• How are the initial density and core size related? For example, does a smaller starting density produce a large or small core radius when compared to a larger starting density?

Density Variation Example

• Find a density value that creates a nan density for the core of your chosen planet? Why does this happen?
• Find a density value where the model shows a core is larger than the planet itself? Why does this happen?

Radius Variation Example

• Find a radius value that creates a nan density for the core of your chosen planet? Why does this happen?
• Find a radius value where the model shows a core is larger than the planet itself? Why does this happen?

Mass Variation Example

• Is there an upper bound on what the mass can be in the model? ie is there a high value that breaks the model like the previous examples; ignore whether it is realistic or not.
• Find a mass value that creates a nan density for the core of your chosen planet? Why does this happen?
• Find a mass value where the model shows a core is larger than the planet itself? Why does this happen?