planet.ipynb

Kernel: Unknown Kernel

import numpy as np
%pylab inline
import matplotlib.pyplot as plt

Populating the interactive namespace from numpy and matplotlib

#RK4 integrator
def RK4_step( f, y, r, h ):
    #Creating solutions
    K0 = h*f(r, y)
    K1 = h*f(r + 0.5*h, y + 0.5*K0)
    K2 = h*f(r + 0.5*h, y + 0.5*K1)
    K3 = h*f(r + h, y + K2)
    y1 = y + 1/6.0*(K0 + 2*K1 + 2*K2 + K3 )
    #Returning solution
    return y1

#Defining Bisection function
def Bisection( f, a, b, Nmax, printer=False ):
    #verifying the STEP1, a and b with different signs
    if f(a)*f(b)>0:
        print "Error, f(a) and f(b) should have opposite signs"
        return False
    #Assigning the current extreme values, STEP2
    ai = a
    bi = b
    #Iterations
    n = 1
    while n<=Nmax:
        #Bisection, STEP3
        pi = (ai+bi)/2.0
        #Evaluating function in pi, STEP4 and STEP5
        if printer:
            print "Value for %d iterations:"%(n),pi
        #Condition A
        if f(pi)*f(ai)>0:
            ai = pi
        #Condition B
        elif f(pi)*f(ai)<0:
            bi = pi
        #Condition C: repeat the cycle
        n+=1
    #Final result
    return pi

#Parameters of the system
#Cavendish constant
G = 6.67e-11
#Bulk Modulus
Ks = 2.8e11

#Surface pressure
P_surf = 1.0e5
#Surface gravity
g_surf = 9.8
#Total mass
M_p = 5.97e24
#Surface density
rho_surf = 3500.0

#Delta radius
h = -10000.0
#Number of iterations for bisection
Nbis = 20

#Dynamic function of the system
def f( r, y ):
    #Extracting values
    P = y[0]
    g = y[1]
    m = y[2]
    rho = y[3]
    #Defining derivatives
    dP = -rho*g
    dg = 4*np.pi*G*rho - 2*G*m/r**3
    dm = 4*np.pi*r**2*rho
    drho = -rho**2*g/Ks
    
    return np.array( [dP, dg, dm, drho] )

def Mr( R ):
    #Array of radius
    rarray = np.arange( R, 0, h )
    #Initializing variables
    P = [ P_surf, ]
    g = [ g_surf, ]
    m = [ M_p, ]
    rho = [ rho_surf, ]
    #Numerical integration
    i = 0
    for r in rarray:
        y = [ P[i], g[i], m[i], rho[i] ]
        Pi, gi, mi, rhoi = RK4_step( f, y, r, h )
        #Storing solutions
        P.append( Pi )
        g.append( gi )
        m.append( mi )
        rho.append( rhoi )
        i += 1
    return m[-1]/M_p

Rarray = np.linspace( 6400000.0, 7000000.0, 400 )
Mrarray = []
for R in Rarray:
    Mrarray.append( Mr(R) )
plt.figure( figsize=(8,8) )
plt.plot( Rarray, Mrarray )
plt.ylim( (-0.1,0.1) )

(-0.1, 0.1)