Planet_model_Ot20.py 3.76 KB
Newer Older
1
2
3
4
import numpy as np
import astropy.units as u
from astropy import constants as const

5

6
def calculate_mass_planet_Ot19(R_pl):
7
8
9
10
11
    """ Estimate a planet mass (in Earth masses) based on the mass-radius
    relation derived for the volatile rich regime in Otegi et al. 2020
    (radius planet larger than 2.115 Earth radii, or density < 3.3 g/cc).
    Mass has to be smaller than 120 Earth masses. Radius should be given
    in Earth radii.
12
    """
13

14
15
16
    M_earth = const.M_earth.cgs.value
    R_earth = const.R_earth.cgs.value

17
18
19
20
21
    if (type(R_pl) == int) or (type(R_pl) == float) \
       or (type(R_pl) == np.float64):  # if R is single value
        M_p_volatile = 1.74*(R_pl)**1.58  # if rho < 3.3 g/cm^3
        rho_volatile = (M_p_volatile*M_earth) \
                        / (4/3*np.pi*(R_pl*R_earth)**3)
22
        if (rho_volatile >= 3.3):
23
24
            raise Exception("Planet with this radius is too small and"\
				 			+ " likely rocky; use LoFo14 models instead.")
25
26
        else:
            if (M_p_volatile >= 120):
27
28
                raise Exception("Planet too massive. M-R relation only valid"\
                                + " for <120 M_earth.")
29
30
31
            else:
                return M_p_volatile

32
    elif len(R_pl) > 1:  # if R is array
33
34
        Ms = []
        for i in range(len(R_pl)):
35
36
37
            M_p_volatile_i = 1.74*(R_pl[i])**1.58  # if rho < 3.3 g/cm^3
            rho_volatile_i = (M_p_volatile_i * M_earth) \
                             / (4/3*np.pi*(R_pl[i]*R_earth)** 3)
38
39
40
41
42
43
44
45
            if (rho_volatile_i >= 3.3) or (M_p_volatile_i >= 120):
                M_i = np.nan
            else:
                M_i = M_p_volatile_i
            Ms.append(M_i)
        Ms = np.array(Ms)
        return Ms

46

47
def calculate_radius_planet_Ot19(M_pl):
48
49
50
51
    """ Estimate a planet radius (in Earth radii) based on the mass-radius
    relation derived for the volatile rich regime in Otegi et al. 2020.
    Mass needs to be between ~5.7 and 120 Earth masses and should be given 
    in Earth masses.
52
    """
53

54
55
56
    M_earth = const.M_earth.cgs.value
    R_earth = const.R_earth.cgs.value

57
58
59
    if (type(M_pl) == int) or (type(M_pl) == float) \
            or (type(M_pl) == np.float64):  # if M is single value
        # R_p_volatile = 0.7*M_pl**0.63 # if rho < 3.3 g/cm^3
60
        R_p_volatile = (M_pl/1.74)**(1./1.58)
61
62
        rho_volatile = (M_pl*M_earth) \
                        / (4/3*np.pi*(R_p_volatile*R_earth)**3)
63
        if (rho_volatile >= 3.3):
64
65
            raise Exception("Planet with this mass/radius is too small and"\
                            + " likely rocky; use LoFo14 models instead.")
66
67
        else:
            if (M_pl >= 120):
68
69
                raise Exception("Planet too massive. M-R relation only"\
                                + " valid for <120 M_earth.")
70
71
72
            else:
                return R_p_volatile

73
    elif len(M_pl) > 1:  # if M is array
74
75
        Rs = []
        for i in range(len(M_pl)):
76
            # R_p_volatile_i = 0.7*M_pl[i]**0.63 # if rho < 3.3 g/cm^3
77
            R_p_volatile_i = (M_pl[i]/1.74)**(1./1.58)
78
79
            rho_volatile_i = (M_pl[i]*M_earth)\
                            / (4/3*np.pi*(R_p_volatile_i*R_earth)**3)
80
81
82
83
84
85
86
            if (rho_volatile_i >= 3.3) or (M_pl[i] >= 120):
                R_i = np.nan
            else:
                R_i = R_p_volatile_i
            Rs.append(R_i)
        Rs = np.array(Rs)
        return Rs
87
88


89
def density_planet(M_p, R_p):
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
    """ Calculate the mean density in cgs units using the radius and
    mass of the planet.

    Parameters:
    -----------
    M_p (float): mass of the planet (in Earth masses)
    R_p (float): radius of the planet (in Earth radii)

    Returns:
    --------
    rho (float): density in cgs units (g/ccm)
    """

    rho = (M_p*const.M_earth.cgs) \
        / (4./3*np.pi*(R_p*const.R_earth.cgs)**3)
    return rho.value