task01.ipynb

Kernel: Python 2

Task 01

import numpy as np
import matplotlib.pylab as plt
from scipy import integrate
%pylab inline

Populating the interactive namespace from numpy and matplotlib

WARNING: pylab import has clobbered these variables: ['plt']
`%matplotlib` prevents importing * from pylab and numpy

This homework is an activity intended to qualify work class and apply the fixed-point iteration method for solving a problem related to an estimation of the radius of an exoplanet.

Due to: Abr 10

---

1) One variable equations methods (30%)

In class several methods to find roots were developed. Submit the codes you did during class for at least 2 of the methods seen in class. Show the roots you find for an specific function.

2) Estimating the Radius of an Exoplanet (70%)

When a new planet is discovered, there are different methods to estimate its physical properties. Many times is only possible to estimate either the planet mass or the planet radius and the other property has to be predicted through computer modelling.

If one has the planet mass, a very rough way to estimate its radius is to assume certain composition (mean density) and a homogeneous distribution (a very bad assumption!). For example, for the planet Gliese 832c with a mass $M= 5.40 M_{\oplus}$, if we assume an earth-like composition, i.e. $\bar \rho_{\oplus} = 5520\ kg/m^3$, we obtain:

$$R_{g832c} = \left( \frac{3 M_{g832c}}{ 4 \pi \bar\rho_{\oplus} } \right)^{1/3} \approx 1.75 R_{\oplus}$$

That would be the planet radius if the composition where exactly equal to earth's.

A more realistic approach is assuming an internal one-layer density profile like:

$$\rho(r) = \rho_0 \exp\left( -\frac{r}{L} \right)$$

where $\rho_0$ is the density at planet centre and $L$ is a characteristic lenght depending on the composition. From numerical models of planet interiors, the estimated parameters for a planet of are $M= 5.40 M_{\oplus}$ are approximately $\rho_0 = 18000\ kg/m^3$ and $L = 6500\ km$.

Integrating over the planet volume, we obtain the total mass as

$$M = 4\pi \int_0^R \rho(r)r^2dr$$

This is a function of the mass in terms of the planet radius.

Solving the equation $M(R) = M_{g832c}$ it would be possible to find a more realistic planet radius. However when using numerical models, it is not possible to approach the solution from the left side as a negative mass makes no sense.

a. Solve the previous problem and find the radius of Gliese 832c using your own version of the Fixed-point iteration algorithm.

##Solution 2)

To find the radius of Gliese, the method fixed point is used. Although the integration can be done analytically, it can also be done numerically. But in this specific case, it is recommended to use an analytical form, not only because it is more efficient but it avoids more errors associated to numerical approximations.

#Constants 
Rt = 6.4E6 #m
Mt = 5.98E24 #kg

#Parameters
L = 6.5E6/Rt # Normalized with Earth radius
Rho0 = 18000*(Rt**3/Mt) # Density normalized with earth mass
Mg = 5.4 #Mass of Gliese normalized with earth mass

# Fix point method
def fix_point( f, Niter, p0, tol = 1e-17 ):
    
    g = lambda x: x-f(x)
    pi = [p0,]
    
    n = 1
    while n<Niter:
        pi.append( g(pi[n-1]) )
        if((pi[n]-pi[n-1])/pi[n])<tol:
            return pi[n]
        n+=1    
    return "The method do not converge with Nmax iterations", Niter

#Integrating over the planet volume for the density distribution
def F1(R):
    return 4*np.pi*Rho0*(2*L**3-L*np.exp(-R/L)*(2*L**2+2*L*R+R**2))-Mg

#Numerical integration
fi = lambda x: exp(-x/L)*x**2
def F2(R):
    return 4*np.pi*Rho0*integrate.quad(fi,0,R)[0]-Mg

# Gliese 832c radius using the expression analytically integrated
fix_point(F1,2000,1.792)

1.7885838956845579

# Gliese 832c radius using the numerical integrated expression
fix_point(F2,2000,1.792)

1.788583895684572