Plot g(r) of planet with two distinct density regions

[oddjobmj]

Homework Statement

Consider a simple model for the interior of the Earth: there is a spherical iron core with constant mass density ρ0 and radius a; outside the core is "rock" with constant density ρ1. Use these values for the densities: ρ0= 8.70×10³ kg/m3 and ρ1= 4.10×10³ kg/m3. The radius of the Earth is R = 6.40×10⁶ m.

Hand in a plot g(r) from 0 to 2R.
The value at r=R must be 9.81 m/s2.

Homework Equations

g(r)=\frac{G*M(r)}{r^2}

a=4.291*10⁶ m

Total Mass=\frac{4πρ0}{3}a^3+\frac{4πρ}{3}(R³-a³)

The Attempt at a Solution

Since the mass is a function of r the function will be affected by the density at the current radius. Is there a way to write a single function that will cover this or do I have to do it with a piecewise function where I use ρ0 between 0 < r < a, then ρ between a < r < R, then constant past R?

 

[cepheid]

rho(r) is by definition piecewise in this problem.

You just need to find the function M(r) in terms of rho(r). M(r) is the mass enclosed by radius r.

Hopefully when you say "constant past R", you mean constant and zero past R.

 

[oddjobmj]

Thank you, cepheid.

When I said constant past r I meant M(r) is constant past r. Yes, the density would be zero.

 

[oddjobmj]

Actually, the function for M(r) would be zero if the density were zero. That would lead to 0/r^2 which would basically mean my g(r) wouldn't change at all beyond that point. I don't believe that should be the case. So, at that point I need to simply calculate the total mass of the Earth and plug that in instead of the function for M(r).

Edit:
Can you tell what I did wrong here?

http://imgur.com/0Bp1aaM

Edit 2:

I think the second function should have (r³-a³) instead of (r-a)^3. All the lines are connected now and it looks as if it exhibits the characteristics one would expect from the system.

 

[cepheid]

Actually, the function for M(r) would be zero if the density were zero.

No it wouldn't. Remember that M(r) is the total mass enclosed by the radius r. This is non-zero for r > R. It's just constant, as you pointed out. Hint: for r > a, both densities come into play in calculating M enclosed, not just one.

 

[oddjobmj]

I'm sorry, I was unclear. You are correct. I didn't mean to imply that M(r) should actually equal zero. I just meant that if I happened to plug in zero for the density it would spit out zero. Instead I just set the mass at r>R to the total mass of the planet.

 

[cepheid]

Cool. Just out of curiosity, what does your new plot look like?

 

[oddjobmj]

http://i.imgur.com/xUNeuXg.png

 

[cepheid]

Cool. I got some very similar looking output (attached png).

#########################################################################
# g_earth.py -- simple model of Earth with const. density iron core for 
#               r < a, and const  density rock mantle for a < r < R.
#########################################################################

import numpy as np
import matplotlib.pyplot as plt
import math


R = 6.40e6      # metres
rho0 = 8.70e3   # kg/m^3
rho1 = 4.10e3   # kg/m^3
a = 4.291e6     # metres
G = 6.67384e-11 # m^3 kg^-1 s^-1

r = np.arange(1.0, 1.0e7, 1.0e3)

M = np.zeros(len(r))
g = np.zeros(len(r))

for i in xrange(len(r)):
    if (r[i] < a):
        M[i] = rho0*(4.0/3.0)*math.pi*(r[i])**3 
        g[i] = (G*M[i])/(r[i])**2
    elif (r[i] < R):
        M[i] = rho1*(4.0/3.0)*math.pi*((r[i])**3 - a**3)\
               + rho0*(4.0/3.0)*math.pi*(a**3)
        g[i] = (G*M[i])/(r[i])**2
    else:
        M[i] = rho1*(4.0/3.0)*math.pi*(R**3 - a**3)\
               + rho0*(4.0/3.0)*math.pi*(a**3)
        g[i] = (G*M[i])/(r[i])**2

plt.figure(1, figsize=(12,8), dpi=100)
plt.plot(r,g)
plt.title("Gravitational Field Strength as a Function of Radial Distance for Simple Earth Model")
plt.ylabel("g(r) [N/kg]")
plt.xlabel("r [m]")
plt.savefig("g_earth.png")
plt.show()

 

[oddjobmj]

Ah, I didn't even think about breaking out IDLE... Also, I'll sleep better tonight knowing you came up with similar results. Thank you!

 

[cepheid]

Since g ~ M/r^2, we can predict what g(r) should look like in the three regimes.

For r < a, M(r) ~ r^3, therefore g ~ r^3/r^2. In other words, g ~ r so the linear relationship is totally expected. Here I use ~ to mean "scales as" or "is proportional to."

For r > R, M(r) = const., therefore g ~ 1/r^2, which looks about right on the plot

The middle phase is tricker.