11979115482

Committed 22 Nov 2024 07:45PM UTC coverage: 65.499% (-0.2%) from 65.737%

Build # 11979115482

Build Type

Pull #403

github

Committed by

web-flow

Commit Message

Merge 0d3dd536f into 86cc76cd9

Pull Request Pull Request #403: Improved version tracking and deprecated SVN

Run Details

36 of 54 new or added lines in 2 files covered. (66.67%)

43 existing lines in 5 files now uncovered.

9534 of 14556 relevant lines covered (65.5%)

0.65 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

94.74

/EXOSIMS/PlanetPhysicalModel/FortneyMarleyCahoyMix1.py

from EXOSIMS.Prototypes.PlanetPhysicalModel import PlanetPhysicalModel
import astropy.units as u
import numpy as np
import scipy.interpolate as interpolate
import os
import inspect
import pickle
from scipy.io import loadmat


class FortneyMarleyCahoyMix1(PlanetPhysicalModel):
    """Planet density models based on Fortney & Marley, albedo models based on
    Cahoy.  Intended for use with the Kepler-like planet population modules.

    Args:
        **specs:
            user specified values

    Attributes:


    Notes:
    1. The calculation of albedo is based solely on the semi-major axis
    and uses a uniform distribution of metallicities to interpolate albedo from
    the grid in Cahoy et al. 2010.

    """

    def __init__(self, **specs):

        PlanetPhysicalModel.__init__(self, **specs)

        # define albedo interpolant:
        smas = [0.8, 2, 5, 10]
        fes = [1, 3, 10, 30]
        ps = np.array(
            [
                [0.322, 0.241, 0.209, 0.142],
                [0.742, 0.766, 0.728, 0.674],
                [0.567, 0.506, 0.326, 0.303],
                [0.386, 0.260, 0.295, 0.279],
            ]
        )

        grid_a, grid_fe = np.meshgrid(smas, fes)
        self.albedo_pts = np.vstack((grid_a.flatten(), grid_fe.flatten())).T
        self.albedo_vals = ps.T.flatten()

        # define conversion functions for icy/rocky/metal planets
        # ice/rock fraction - Fortney et al., 2007 Eq. 7 (as corrected in paper erratum)
        # units are in Earth masses and radii, frac = 1 is pure ice
        self.R_ir = (
            lambda frac, M: (0.0912 * frac + 0.1603) * np.log10(M) ** 2.0
            + (0.333 * frac + 0.7387) * np.log10(M)
            + (0.4639 * frac + 1.1193)
        )
        self.M_ir = lambda frac, R: 10.0 ** (
            (
                -(0.333 * frac + 0.7387)
                + np.sqrt(
                    (0.333 * frac + 0.7387) ** 2.0
                    - 4.0 * (0.0912 * frac + 0.1603) * (0.4639 * frac + 1.1193 - R)
                )
            )
            / (2.0 * (0.0912 * frac + 0.1603))
        )

        # rock/iron fraction - Fortney et al., 2007 Eq. 8 (as corrected in paper
        # erratum). units are in Earth masses and radii, frac = 1 is pure rock
        self.R_ri = (
            lambda frac, M: (0.0592 * frac + 0.0975) * np.log10(M) ** 2.0
            + (0.2337 * frac + 0.4938) * np.log10(M)
            + (0.3102 * frac + 0.7932)
        )
        self.M_ri = lambda frac, R: 10.0 ** (
            (
                -(0.2337 * frac + 0.4938)
                + np.sqrt(
                    (0.2337 * frac + 0.4938) ** 2.0
                    - 4.0 * (0.0592 * frac + 0.0975) * (0.3102 * frac + 0.7932 - R)
                )
            )
            / (2.0 * (0.0592 * frac + 0.0975))
        )

        # find and load the Fortney et al. Table 4 data (gas giant densities)
        # data is al in Jupiter radii, Earth masses and AU
        filename = "Fortney_etal_2007_table4.p"
        datapath = os.path.join(self.cachedir, filename)
        if os.path.exists(datapath):
            try:
                with open(datapath, "rb") as f:
                    self.ggdat = pickle.load(f)
            except UnicodeDecodeError:
                with open(datapath, "rb") as f:
                    self.ggdat = pickle.load(f, encoding="latin1")
        else:
            classpath = os.path.split(inspect.getfile(self.__class__))[0]
            matpath = os.path.join(classpath, "fortney_table4.mat")
            self.ggdat = loadmat(matpath, squeeze_me=True, struct_as_record=False)
            with open(datapath, "wb") as f:
                pickle.dump(self.ggdat, f)

        self.ggdat["dist"] = self.ggdat["dist"] * u.AU
        self.ggdat["planet_mass"] = self.ggdat["planet_mass"] * u.earthMass
        self.ggdat["radii"] = self.ggdat["radii"] * u.jupiterRad
        # convert to Earth radius
        Rtmp = self.ggdat["radii"].to("earthRad").value
        # remove NANs
        Rtmp[np.isnan(Rtmp)] = 0.0
        self.giant_pts = np.array(
            [
                self.ggdat["x1"].flatten().astype(float),
                self.ggdat["x3"].flatten().astype(float),
                Rtmp.flatten().astype(float),
            ]
        ).T
        self.giant_vals = self.ggdat["x2"].flatten()

        self.giant_pts2 = np.array(
            [
                self.ggdat["x1"].flatten().astype(float),
                self.ggdat["x3"].flatten().astype(float),
                self.ggdat["x2"].flatten().astype(float),
            ]
        ).T
        self.giant_vals2 = Rtmp.flatten().astype(float)

    def calc_albedo_from_sma(self, a, prange=None):
        """Helper function for calculating albedo.

        We assume a uniform distribution of metallicities, and then interpolate the
        grid from Cahoy et al. 2010.

        Args:
            a (astropy Quanitity array):
               Semi-major axis values

        Args:
            a (~astropy.units.Quantity(~numpy.ndarray(float))):
               Semi-major axis values
            prange (arrayLike):
                [Min, Max] geometric albedo. Defaults to [0.367, 0.367]

        Returns:
            ~numpy.ndarray(float):
                Albedo values

        .. warning::
            The prange input is ignored.

        """

        # grab the sma values and constrain to grid
        atmp = np.clip(a.to("AU").value, 0.8, 10.0)

        # generate uniform fe grid:
        fetmp = np.random.uniform(size=atmp.size, low=1, high=30)

        p = interpolate.griddata(
            self.albedo_pts, self.albedo_vals, (atmp, fetmp), method="cubic"
        )

        return p

    def calc_mass_from_radius(self, Rp):
        """Helper function for calculating mass given the radius.

        The calculation is done in two steps, first covering all things that can
        only be ice/rock/iron, and then things that can be giants.

        Args:
            Rp (astropy Quantity array):
                Planet radius in units of Earth radius

        Returns:
            Mp (astropy Quantity array):
                Planet mass in units of Earth mass

        """

        Mp = np.zeros(Rp.shape)

        # first, the things up to the min giant radius but greater
        # than 2 Earth radii (assumed to be icy)
        inds = (Rp <= np.nanmin(self.ggdat["radii"])) & (Rp > 2 * u.earthRad)
        Rtmp = Rp[inds]
        fracs = np.random.uniform(size=Rtmp.size, low=0.5, high=1.0)
        Mp[inds] = self.M_ir(fracs, Rtmp.to("earthRad").value)

        # everything under 2 Earth radii can by ice/rock/iron
        inds = Rp <= 2 * u.earthRad
        Rtmp = Rp[inds]
        Mtmp = np.zeros(Rtmp.shape)
        fracs = np.random.uniform(size=Rtmp.size, low=-1.0, high=1.0)

        # ice/rock and rock/iron
        icerock = fracs < 0
        Mtmp[icerock] = self.M_ir(
            np.abs(fracs[icerock]), Rtmp[icerock].to("earthRad").value
        )
        rockiron = fracs >= 0
        Mtmp[rockiron] = self.M_ri(
            np.abs(fracs[rockiron]), Rtmp[rockiron].to("earthRad").value
        )
        Mp[inds] = Mtmp

        # everything else is a giant.  those above the table limit
        # are inflated close-in things that are undetectable
        inds = Rp > np.nanmax(self.ggdat["radii"])
        Mp[inds] = np.max(self.ggdat["planet_mass"]).to("earthMass").value

        inds = (Rp > np.nanmin(self.ggdat["radii"])) & (
            Rp <= np.nanmax(self.ggdat["radii"])
        )
        Rtmp = Rp[inds]
        Mtmp = interpolate.griddata(
            self.giant_pts,
            self.giant_vals,
            (
                np.random.uniform(low=0, high=100, size=Rtmp.size),
                np.exp(
                    np.log(0.02)
                    + (np.log(9.5) - np.log(0.02)) * np.random.uniform(size=Rtmp.size)
                ),
                Rtmp.to("earthRad").value,
            ),
        )
        if np.any(np.isnan(Mtmp)):
            inds2 = np.isnan(Mtmp)
            Mtmp[inds2] = (
                (1.33 * u.g / u.cm**3.0 * 4 * np.pi * Rtmp[inds2] ** 3.0 / 3.0)
                .to("earthMass")
                .value
            )
        Mp[inds] = Mtmp

        Mp = Mp * u.earthMass

        return Mp

    def calc_radius_from_mass(self, Mp):
        """Helper function for calculating radius given the mass.

        The calculation is done in two steps, first covering all things that can
        only be ice/rock/iron, and then things that can be giants.

        Args:
            Mp (astropy Quantity array):
                Planet mass in units of Earth mass

        Returns:
            Rp (astropy Quantity array):
                Planet radius in units of Earth radius

        """

        Rp = np.zeros(Mp.shape)

        # Everything below the tabulated mass values is treated as
        # ice/rock/iron
        inds = Mp <= np.nanmin(self.ggdat["planet_mass"])
        Mtmp = Mp[inds]
        Rtmp = np.zeros(Mtmp.shape)
        fracs = np.random.uniform(size=Mtmp.size, low=-1.0, high=1.0)

        # ice/rock and rock/iron
        icerock = fracs < 0
        Rtmp[icerock] = self.R_ir(
            np.abs(fracs[icerock]), Mtmp[icerock].to("earthMass").value
        )
        rockiron = fracs >= 0
        Rtmp[rockiron] = self.R_ri(
            np.abs(fracs[rockiron]), Mtmp[rockiron].to("earthMass").value
        )
        Rp[inds] = Rtmp

        # everything else is a giant.
        inds = Mp > np.nanmin(self.ggdat["planet_mass"])
        Mtmp = Mp[inds]
        Rp[inds] = interpolate.griddata(
            self.giant_pts2,
            self.giant_vals2,
            (
                np.random.uniform(low=0, high=100, size=Mtmp.size),
                np.exp(
                    np.log(0.02)
                    + (np.log(9.5) - np.log(0.02)) * np.random.uniform(size=Mtmp.size)
                ),
                Mtmp.to("earthMass").value,
            ),
        )

        # things that failed
        inds = np.isnan(Rp) | (Rp == 0.0)
        if np.any(inds):
            rho = 1.33 * u.g / u.cm**3.0
            Rp[inds] = (
                ((3 * Mp[inds] / rho / np.pi / 4.0) ** (1 / 3.0)).to("earthRad").value
            )

        Rp = Rp * u.earthRad

        return Rp

1	from EXOSIMS.Prototypes.PlanetPhysicalModel import PlanetPhysicalModel	1✔
2	import astropy.units as u	1✔
3	import numpy as np	1✔
4	import scipy.interpolate as interpolate	1✔
5	import os	1✔
6	import inspect	1✔
7	import pickle	1✔
8	from scipy.io import loadmat	1✔
9
10
11	class FortneyMarleyCahoyMix1(PlanetPhysicalModel):	1✔
12	"""Planet density models based on Fortney & Marley, albedo models based on
13	Cahoy. Intended for use with the Kepler-like planet population modules.
14
15	Args:
16	**specs:
17	user specified values
18
19	Attributes:
20
21
22	Notes:
23	1. The calculation of albedo is based solely on the semi-major axis
24	and uses a uniform distribution of metallicities to interpolate albedo from
25	the grid in Cahoy et al. 2010.
26
27	"""
28
29	def __init__(self, **specs):	1✔
30
31	PlanetPhysicalModel.__init__(self, **specs)	1✔
32
33	# define albedo interpolant:
34	smas = [0.8, 2, 5, 10]	1✔
35	fes = [1, 3, 10, 30]	1✔
36	ps = np.array(	1✔
37	[
38	[0.322, 0.241, 0.209, 0.142],
39	[0.742, 0.766, 0.728, 0.674],
40	[0.567, 0.506, 0.326, 0.303],
41	[0.386, 0.260, 0.295, 0.279],
42	]
43	)
44
45	grid_a, grid_fe = np.meshgrid(smas, fes)	1✔
46	self.albedo_pts = np.vstack((grid_a.flatten(), grid_fe.flatten())).T	1✔
47	self.albedo_vals = ps.T.flatten()	1✔
48
49	# define conversion functions for icy/rocky/metal planets
50	# ice/rock fraction - Fortney et al., 2007 Eq. 7 (as corrected in paper erratum)
51	# units are in Earth masses and radii, frac = 1 is pure ice
52	self.R_ir = (	1✔
53	lambda frac, M: (0.0912 * frac + 0.1603) * np.log10(M) ** 2.0
54	+ (0.333 * frac + 0.7387) * np.log10(M)
55	+ (0.4639 * frac + 1.1193)
56	)
57	self.M_ir = lambda frac, R: 10.0 ** (	1✔
58	(
59	-(0.333 * frac + 0.7387)
60	+ np.sqrt(
61	(0.333 * frac + 0.7387) ** 2.0
62	- 4.0 * (0.0912 * frac + 0.1603) * (0.4639 * frac + 1.1193 - R)
63	)
64	)
65	/ (2.0 * (0.0912 * frac + 0.1603))
66	)
67
68	# rock/iron fraction - Fortney et al., 2007 Eq. 8 (as corrected in paper
69	# erratum). units are in Earth masses and radii, frac = 1 is pure rock
70	self.R_ri = (	1✔
71	lambda frac, M: (0.0592 * frac + 0.0975) * np.log10(M) ** 2.0
72	+ (0.2337 * frac + 0.4938) * np.log10(M)
73	+ (0.3102 * frac + 0.7932)
74	)
75	self.M_ri = lambda frac, R: 10.0 ** (	1✔
76	(
77	-(0.2337 * frac + 0.4938)
78	+ np.sqrt(
79	(0.2337 * frac + 0.4938) ** 2.0
80	- 4.0 * (0.0592 * frac + 0.0975) * (0.3102 * frac + 0.7932 - R)
81	)
82	)
83	/ (2.0 * (0.0592 * frac + 0.0975))
84	)
85
86	# find and load the Fortney et al. Table 4 data (gas giant densities)
87	# data is al in Jupiter radii, Earth masses and AU
88	filename = "Fortney_etal_2007_table4.p"	1✔
89	datapath = os.path.join(self.cachedir, filename)	1✔
90	if os.path.exists(datapath):	1✔
91	try:	1✔
92	with open(datapath, "rb") as f:	1✔
93	self.ggdat = pickle.load(f)	1✔
94	except UnicodeDecodeError:	×
95	with open(datapath, "rb") as f:	×
96	self.ggdat = pickle.load(f, encoding="latin1")	×
97	else:
98	classpath = os.path.split(inspect.getfile(self.__class__))[0]	1✔
99	matpath = os.path.join(classpath, "fortney_table4.mat")	1✔
100	self.ggdat = loadmat(matpath, squeeze_me=True, struct_as_record=False)	1✔
101	with open(datapath, "wb") as f:	1✔
102	pickle.dump(self.ggdat, f)	1✔
103
104	self.ggdat["dist"] = self.ggdat["dist"] * u.AU	1✔
105	self.ggdat["planet_mass"] = self.ggdat["planet_mass"] * u.earthMass	1✔
106	self.ggdat["radii"] = self.ggdat["radii"] * u.jupiterRad	1✔
107	# convert to Earth radius
108	Rtmp = self.ggdat["radii"].to("earthRad").value	1✔
109	# remove NANs
110	Rtmp[np.isnan(Rtmp)] = 0.0	1✔
111	self.giant_pts = np.array(	1✔
112	[
113	self.ggdat["x1"].flatten().astype(float),
114	self.ggdat["x3"].flatten().astype(float),
115	Rtmp.flatten().astype(float),
116	]
117	).T
118	self.giant_vals = self.ggdat["x2"].flatten()	1✔
119
120	self.giant_pts2 = np.array(	1✔
121	[
122	self.ggdat["x1"].flatten().astype(float),
123	self.ggdat["x3"].flatten().astype(float),
124	self.ggdat["x2"].flatten().astype(float),
125	]
126	).T
127	self.giant_vals2 = Rtmp.flatten().astype(float)	1✔
128
129	def calc_albedo_from_sma(self, a, prange=None):	1✔
130	"""Helper function for calculating albedo.
131
132	We assume a uniform distribution of metallicities, and then interpolate the
133	grid from Cahoy et al. 2010.
134
135	Args:
136	a (astropy Quanitity array):
137	Semi-major axis values
138
139	Args:
140	a (~astropy.units.Quantity(~numpy.ndarray(float))):
141	Semi-major axis values
142	prange (arrayLike):
143	[Min, Max] geometric albedo. Defaults to [0.367, 0.367]
144
145	Returns:
146	~numpy.ndarray(float):
147	Albedo values
148
149	.. warning::
150	The prange input is ignored.
151
152	"""
153
154	# grab the sma values and constrain to grid
155	atmp = np.clip(a.to("AU").value, 0.8, 10.0)	1✔
156
157	# generate uniform fe grid:
158	fetmp = np.random.uniform(size=atmp.size, low=1, high=30)	1✔
159
160	p = interpolate.griddata(	1✔
161	self.albedo_pts, self.albedo_vals, (atmp, fetmp), method="cubic"
162	)
163
164	return p	1✔
165
166	def calc_mass_from_radius(self, Rp):	1✔
167	"""Helper function for calculating mass given the radius.
168
169	The calculation is done in two steps, first covering all things that can
170	only be ice/rock/iron, and then things that can be giants.
171
172	Args:
173	Rp (astropy Quantity array):
174	Planet radius in units of Earth radius
175
176	Returns:
177	Mp (astropy Quantity array):
178	Planet mass in units of Earth mass
179
180	"""
181
182	Mp = np.zeros(Rp.shape)	1✔
183
184	# first, the things up to the min giant radius but greater
185	# than 2 Earth radii (assumed to be icy)
186	inds = (Rp <= np.nanmin(self.ggdat["radii"])) & (Rp > 2 * u.earthRad)	1✔
187	Rtmp = Rp[inds]	1✔
188	fracs = np.random.uniform(size=Rtmp.size, low=0.5, high=1.0)	1✔
189	Mp[inds] = self.M_ir(fracs, Rtmp.to("earthRad").value)	1✔
190
191	# everything under 2 Earth radii can by ice/rock/iron
192	inds = Rp <= 2 * u.earthRad	1✔
193	Rtmp = Rp[inds]	1✔
194	Mtmp = np.zeros(Rtmp.shape)	1✔
195	fracs = np.random.uniform(size=Rtmp.size, low=-1.0, high=1.0)	1✔
196
197	# ice/rock and rock/iron
198	icerock = fracs < 0	1✔
199	Mtmp[icerock] = self.M_ir(	1✔
200	np.abs(fracs[icerock]), Rtmp[icerock].to("earthRad").value
201	)
202	rockiron = fracs >= 0	1✔
203	Mtmp[rockiron] = self.M_ri(	1✔
204	np.abs(fracs[rockiron]), Rtmp[rockiron].to("earthRad").value
205	)
206	Mp[inds] = Mtmp	1✔
207
208	# everything else is a giant. those above the table limit
209	# are inflated close-in things that are undetectable
210	inds = Rp > np.nanmax(self.ggdat["radii"])	1✔
211	Mp[inds] = np.max(self.ggdat["planet_mass"]).to("earthMass").value	1✔
212
213	inds = (Rp > np.nanmin(self.ggdat["radii"])) & (	1✔
214	Rp <= np.nanmax(self.ggdat["radii"])
215	)
216	Rtmp = Rp[inds]	1✔
217	Mtmp = interpolate.griddata(	1✔
218	self.giant_pts,
219	self.giant_vals,
220	(
221	np.random.uniform(low=0, high=100, size=Rtmp.size),
222	np.exp(
223	np.log(0.02)
224	+ (np.log(9.5) - np.log(0.02)) * np.random.uniform(size=Rtmp.size)
225	),
226	Rtmp.to("earthRad").value,
227	),
228	)
229	if np.any(np.isnan(Mtmp)):	1✔
UNCOV 230	inds2 = np.isnan(Mtmp)	×
UNCOV 231	Mtmp[inds2] = (	×
232	(1.33 * u.g / u.cm*3.0 4 * np.pi * Rtmp[inds2] ** 3.0 / 3.0)
233	.to("earthMass")
234	.value
235	)
236	Mp[inds] = Mtmp	1✔
237
238	Mp = Mp * u.earthMass	1✔
239
240	return Mp	1✔
241
242	def calc_radius_from_mass(self, Mp):	1✔
243	"""Helper function for calculating radius given the mass.
244
245	The calculation is done in two steps, first covering all things that can
246	only be ice/rock/iron, and then things that can be giants.
247
248	Args:
249	Mp (astropy Quantity array):
250	Planet mass in units of Earth mass
251
252	Returns:
253	Rp (astropy Quantity array):
254	Planet radius in units of Earth radius
255
256	"""
257
258	Rp = np.zeros(Mp.shape)	1✔
259
260	# Everything below the tabulated mass values is treated as
261	# ice/rock/iron
262	inds = Mp <= np.nanmin(self.ggdat["planet_mass"])	1✔
263	Mtmp = Mp[inds]	1✔
264	Rtmp = np.zeros(Mtmp.shape)	1✔
265	fracs = np.random.uniform(size=Mtmp.size, low=-1.0, high=1.0)	1✔
266
267	# ice/rock and rock/iron
268	icerock = fracs < 0	1✔
269	Rtmp[icerock] = self.R_ir(	1✔
270	np.abs(fracs[icerock]), Mtmp[icerock].to("earthMass").value
271	)
272	rockiron = fracs >= 0	1✔
273	Rtmp[rockiron] = self.R_ri(	1✔
274	np.abs(fracs[rockiron]), Mtmp[rockiron].to("earthMass").value
275	)
276	Rp[inds] = Rtmp	1✔
277
278	# everything else is a giant.
279	inds = Mp > np.nanmin(self.ggdat["planet_mass"])	1✔
280	Mtmp = Mp[inds]	1✔
281	Rp[inds] = interpolate.griddata(	1✔
282	self.giant_pts2,
283	self.giant_vals2,
284	(
285	np.random.uniform(low=0, high=100, size=Mtmp.size),
286	np.exp(
287	np.log(0.02)
288	+ (np.log(9.5) - np.log(0.02)) * np.random.uniform(size=Mtmp.size)
289	),
290	Mtmp.to("earthMass").value,
291	),
292	)
293
294	# things that failed
295	inds = np.isnan(Rp) \| (Rp == 0.0)	1✔
296	if np.any(inds):	1✔
297	rho = 1.33 * u.g / u.cm**3.0	1✔
298	Rp[inds] = (	1✔
299	((3 * Mp[inds] / rho / np.pi / 4.0) ** (1 / 3.0)).to("earthRad").value
300	)
301
302	Rp = Rp * u.earthRad	1✔
303
304	return Rp	1✔

dsavransky / EXOSIMS / 11979115482

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous