planets.hpp

Go to the documentation of this file.

/** \file planets.hpp
 * \author Jared R. Males (jaredmales@gmail.com)
 * \brief Various utilities related to planets.
 * \ingroup astrofiles
 *
 */
 
#ifndef __mx_astro_planets_hpp__
#define __mx_astro_planets_hpp__
 
#include "../math/constants.hpp"
 
#include "units.hpp"
 
namespace mx
{
namespace astro
{
 
/** \ingroup planets
 * @{
 */
 
/// An ad-hoc planetary mass-to-radius relationship (old version)
/** The goal of this function is to provide a radius given an exoplanet mass, for lightly-irradiated exoplanets.  By
  lightly-irradiated we mean (roughly) planet's at
  * Mercury's separation or further, scaled for stellar luminosity.
  * Here we make use of the transition from rocky to gas-dominated composition at \f$ 1.6  R_e \f$ identified by Rogers
  \cite rogers_2015
  * (see also Marcy et al. (2014) \cite marcy_2014).  Below this radius we assume Earth composition and so
  * \f$ R \propto M^{1/3}\f$.  Above this we scale with a power law matched to the mean radius and mass of Uranus and
  Neptune, which defines the radius between
  * \f$ 1.6^3 M_\oplus \f$ and this Uranus-Neptune mean point.  Above this point we use
  * a polynomial fit (in log(M)) to points including the Uranus/Neptune mean, Saturn, Jupiter, and above Jupiter's mass
  the average points from the 4.5 Gyr 1 AU models from
  * Fortney et al. (2007) \cite fortney_2007.  Above 3591.1 \f$ M_\oplus \f$ (\f$\sim 11 M_{jup}\f$) we scale as \f$
  M^{-1/8} \f$ based on the curve shown in Fortney et al. (2011) \cite fortney_2010.
  *
  * \f$
    \frac{R}{R_\oplus} =
    \begin{cases}
      \left(\frac{M}{M_\oplus}\right)^{1/3}, & M < 4.1 M_\oplus \\
      0.62\left(\frac{M}{M_\oplus}\right)^{0.67}, & 4.1 M_\oplus \le M < 15.84 M_\oplus \\
      14.0211 - 44.8414 \log_{10}\left(\frac{M}{M_\oplus}\right) + 53.6554 \log_{10}^2\left(\frac{M}{M_\oplus}\right)
        -25.3289\log_{10}^3\left(\frac{M}{M_\oplus}\right) + 5.4920\log_{10}^4\left(\frac{M}{M_\oplus}\right) - 0.4586
  \log_{10}^5\left(\frac{M}{M_\oplus}\right), & 15.84 \le M < 3591.1 M_\oplus \\ 32.03
  \left(\frac{M}{M_\oplus}\right)^{-1/8}, & 3591.1 M_\oplus \le M \end{cases} \f$
  *
  * \image html planet_mrdiag_2020.png "The ad-hoc mass-to-radius relationship compared to known planets.  Blue
  circles are the Solar system.  Black circles indicate expoplanets with Teq below that of a blackbody at Mercury's
  separation (413 K)."
  *
  * This function makes use of the units type system (\ref astrounits) so it can be used with Earth masses, Jupiter
  masses, kg (SI units), etc.
  *
  * \returns the estimated radius of the planet.
  *
  * \tparam units is the units-type specifying the units of mass.  See \ref astrounits.
  */
template <typename units>
typename units::realT planetMass2Radius( typename units::realT mass /**< The mass of the planet. */ )
{
    typedef typename units::realT realT;
 
    using namespace mx::astro::constants;
 
    if( mass < 4.1 * massEarth<units>() )
    {
        return pow( mass / massEarth<units>(), math::third<realT>() ) * radEarth<units>();
    }
    else if( mass < static_cast<realT>( 15.84 ) * massEarth<units>() )
    {
        return 0.62 * pow( mass / massEarth<units>(), static_cast<realT>( 0.67 ) ) * radEarth<units>();
    }
    else if( mass < 3591.1 * massEarth<units>() )
    {
        realT logM = log10( mass / massEarth<units>() );
 
        return ( static_cast<realT>( 14.0211 ) - static_cast<realT>( 44.8414 ) * logM +
                 static_cast<realT>( 53.6554 ) * pow( logM, 2 ) - static_cast<realT>( 25.3289 ) * pow( logM, 3 ) +
                 static_cast<realT>( 5.4920 ) * pow( logM, 4 ) - 0.4586 * pow( logM, 5 ) ) *
               radEarth<units>();
    }
    else
    {
        return static_cast<realT>( 32.03 ) * pow( mass / massEarth<units>(), static_cast<realT>( -0.125 ) ) *
               radEarth<units>();
    }
}
typename units::realT planetMass2Radius( typename units::realT mass /**< The mass of the planet. */ ) {…}
 
/// An ad-hoc planetary mass-to-radius relationship.
/** The goal of this function is to provide a radius given an exoplanet mass, for cool, lightly-irradiated exoplanets.
 * Here cool means \f$T_{eq} < 1000 \f$ K.
 *
 * We make use of the transition from rocky to gas-dominated composition at \f$ 1.6  R_\oplus \f$ identified by Rogers
 * \cite rogers_2015 (see also Marcy et al. (2014) \cite marcy_2014).  Below this radius we assume Earth composition and
 * so \f$ R \propto M^{1/3}\f$.
 *
 * Above \f$ 15 M_\oplus\f$ we use the empirical relationship of Thorngren et al (2019) \cite{Thorngren_2019}.
 *
 * Between these two cases we scale with a power law matched to the endpoints.
 *
 * Above  \f$ 12 M_{Jup} we scale as \f$ M^{-1/8} \f$ based on the curve shown in Fortney et al. (2011) \cite
 * fortney_2010.
 *
 * \image html planet_mrdiag_2020.png "The ad-hoc mass-to-radius relationship compared to known planets.  Blue
 * circles are the Solar system.  Points with error bars are from the NASA exoplanet catalog, selected for \f$T_{eq} <
 * 1000\f$ K
 *
 * This function makes use of the units type system (\ref astrounits) so it can be used with Earth masses, Jupiter
 * masses, kg (SI units), etc.
 *
 * \returns the estimated radius of the planet.
 *
 * \tparam units is the units-type specifying the units of mass.  See \ref astrounits.
 */
template <typename units>
typename units::realT planetMass2RadiusWThorngren( typename units::realT mass /**< The mass of the planet. */ )
{
    typedef typename units::realT realT;
 
    using namespace mx::astro::constants;
 
    if( mass < 4.1 * massEarth<units>() )
    {
        return pow( mass / massEarth<units>(), math::third<realT>() ) * radEarth<units>();
    }
    else if( mass < static_cast<realT>( 15.00 ) * massEarth<units>() )
    {
        return 0.6413 * pow( mass / massEarth<units>(), static_cast<realT>( 0.6477 ) ) * radEarth<units>();
    }
    else if( mass < 12 * massJupiter<units>() )
    {
        realT lM = log10( mass / massJupiter<units>() );
 
        return ( static_cast<realT>( 0.96 ) + static_cast<realT>( 0.21 ) * lM -
                 static_cast<realT>( 0.20 ) * pow( lM, 2 ) ) *
               radJupiter<units>();
    }
    else
    {
        return static_cast<realT>( 29.97 ) * pow( mass / massEarth<units>(), static_cast<realT>( -0.125 ) ) *
               radEarth<units>();
    }
}
 
/// Empirical mass-to-radius relationship for cool EGPs from Thorngren et al. (2019)
/** A polynomial in log(M) fit to cool (T < 1000 K) exoplanets by Thorngren et al. (2019) \cite Thorngren_2019.
 * Only valid for \f$ 15 M_\oplus < M < 12 M_{Jup} \f$
 *
 * \returns the estimated radius of the planet.
 *
 * \tparam units is the units-type specifying the units of mass.  See \ref astrounits.
 */
template <typename units>
typename units::realT planetMass2RadiusThorngren( typename units::realT mass /**< The mass of the planet. */ )
{
    typedef typename units::realT realT;
 
    using namespace mx::astro::constants;
 
    realT lM = log10( mass / massJupiter<units>() );
 
    return ( 0.96 + 0.21 * lM - 0.2 * pow( lM, 2 ) ) * radJupiter<units>();
}
 
/// The planetary mass-to-radius of Fabrycky et al. (2014)..
/** A broken power law for low mass planets from Fabrycky et al. (2014)\cite fabrycky_2014 (see also \cite
  lissauer_2011}.
  *
  * \f$
    \frac{R}{R_e} =
    \begin{cases}
      \left(\frac{M}{M_e}\right)^{1/3}, & M < 1 M_e \\
      \left(\frac{M}{M_e}\right)^{1/2.06}, & 1 M_e \le M
      \end{cases}
     \f$
  *
  * This makes use of the units type system (\ref astrounits) so it can be used with Earth masses, Jupiter masses, kg
  (SI units), etc.
  *
  * \returns the estimated radius of the planet.
  *
  * \tparam units is the units-type specifying the units of mass.  See \ref astrounits.
  */
template <typename units>
typename units::realT planetMass2RadiusFab2014( typename units::realT mass /**< The mass of the planet. */ )
{
    typedef typename units::realT realT;
 
    using namespace mx::astro::constants;
 
    if( mass < massEarth<units>() )
    {
        return pow( mass / massEarth<units>(), math::third<realT>() ) * radEarth<units>();
    }
    else
    {
        return pow( mass / massEarth<units>(), static_cast<realT>( 1 ) / static_cast<realT>( 2.06 ) ) *
               radEarth<units>();
    }
}
 
///@}
} // namespace astro
} // namespace mx
#endif