Planet_model_Ot20.py 3.76 KB
 Laura Ketzer committed May 07, 2020 1 2 3 4 ``````import numpy as np import astropy.units as u from astropy import constants as const `````` Laura Ketzer committed Jun 19, 2020 5 `````` `````` Laura Ketzer committed May 07, 2020 6 ``````def calculate_mass_planet_Ot19(R_pl): `````` Laura Ketzer committed Jun 19, 2020 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. `````` Laura Ketzer committed May 07, 2020 12 `````` """ `````` Laura Ketzer committed Jun 19, 2020 13 `````` `````` Laura Ketzer committed May 07, 2020 14 15 16 `````` M_earth = const.M_earth.cgs.value R_earth = const.R_earth.cgs.value `````` Laura Ketzer committed Jun 19, 2020 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) `````` Laura Ketzer committed May 07, 2020 22 `````` if (rho_volatile >= 3.3): `````` Laura Ketzer committed Jun 24, 2020 23 24 `````` raise Exception("Planet with this radius is too small and"\ + " likely rocky; use LoFo14 models instead.") `````` Laura Ketzer committed May 07, 2020 25 26 `````` else: if (M_p_volatile >= 120): `````` Laura Ketzer committed Jun 24, 2020 27 28 `````` raise Exception("Planet too massive. M-R relation only valid"\ + " for <120 M_earth.") `````` Laura Ketzer committed May 07, 2020 29 30 31 `````` else: return M_p_volatile `````` Laura Ketzer committed Jun 19, 2020 32 `````` elif len(R_pl) > 1: # if R is array `````` Laura Ketzer committed May 07, 2020 33 34 `````` Ms = [] for i in range(len(R_pl)): `````` Laura Ketzer committed Jun 19, 2020 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) `````` Laura Ketzer committed May 07, 2020 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 `````` Laura Ketzer committed Jun 19, 2020 46 `````` `````` Laura Ketzer committed May 07, 2020 47 ``````def calculate_radius_planet_Ot19(M_pl): `````` Laura Ketzer committed Jun 19, 2020 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. `````` Laura Ketzer committed May 07, 2020 52 `````` """ `````` Laura Ketzer committed Jun 19, 2020 53 `````` `````` Laura Ketzer committed May 07, 2020 54 55 56 `````` M_earth = const.M_earth.cgs.value R_earth = const.R_earth.cgs.value `````` Laura Ketzer committed Jun 19, 2020 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 `````` Laura Ketzer committed May 07, 2020 60 `````` R_p_volatile = (M_pl/1.74)**(1./1.58) `````` Laura Ketzer committed Jun 19, 2020 61 62 `````` rho_volatile = (M_pl*M_earth) \ / (4/3*np.pi*(R_p_volatile*R_earth)**3) `````` Laura Ketzer committed May 07, 2020 63 `````` if (rho_volatile >= 3.3): `````` Laura Ketzer committed Jun 24, 2020 64 65 `````` raise Exception("Planet with this mass/radius is too small and"\ + " likely rocky; use LoFo14 models instead.") `````` Laura Ketzer committed May 07, 2020 66 67 `````` else: if (M_pl >= 120): `````` Laura Ketzer committed Jun 24, 2020 68 69 `````` raise Exception("Planet too massive. M-R relation only"\ + " valid for <120 M_earth.") `````` Laura Ketzer committed May 07, 2020 70 71 72 `````` else: return R_p_volatile `````` Laura Ketzer committed Jun 19, 2020 73 `````` elif len(M_pl) > 1: # if M is array `````` Laura Ketzer committed May 07, 2020 74 75 `````` Rs = [] for i in range(len(M_pl)): `````` Laura Ketzer committed Jun 19, 2020 76 `````` # R_p_volatile_i = 0.7*M_pl[i]**0.63 # if rho < 3.3 g/cm^3 `````` Laura Ketzer committed May 07, 2020 77 `````` R_p_volatile_i = (M_pl[i]/1.74)**(1./1.58) `````` Laura Ketzer committed Jun 19, 2020 78 79 `````` rho_volatile_i = (M_pl[i]*M_earth)\ / (4/3*np.pi*(R_p_volatile_i*R_earth)**3) `````` Laura Ketzer committed May 07, 2020 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 `````` Laura Ketzer committed Jun 19, 2020 87 88 `````` `````` Laura Ketzer committed May 07, 2020 89 ``````def density_planet(M_p, R_p): `````` Laura Ketzer committed Jun 19, 2020 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``````