stars.hpp

Go to the documentation of this file.

/** \file stars.hpp
 * \author Jared R. Males (jaredmales@gmail.com)
 * \brief Various utilities related to stars.
 * \ingroup astrofiles
 *
 */
 
#include "../math/gslInterpolator.hpp"
 
#include "units.hpp"
 
#ifndef mx_astro_stars_hpp
#define mx_astro_stars_hpp
 
#include <cmath>
 
#include "../ioutils/stringUtils.hpp"
#include "../ioutils/readColumns.hpp"
 
namespace mx
{
namespace astro
{
 
/// Parse a main sequence spectral type string into a numeric code.
/** Expects standard spectral type strings such as "O5V" or "F2.5" or "G8.5V".
 * The numeric code starts at 0 for "O0V", 10 for "A0V", through 90 for "Y0".  The subtypes are added to this value. For
 * instance, "G8.5V" becomes 48.5.
 *
 * Only works for main sequence types.  Any other spectral type, such as MK class I-IV, will return -1.
 *
 * \ingroup stars
 */
template <typename realT = float>
realT numSpType( std::string spType /**< [in] The spectral type string to parse.*/ )
{
 
    spType = ioutils::removeWhiteSpace( spType );
 
    spType = ioutils::toUpper( spType );
 
    realT num = 0;
 
    switch( spType[0] )
    {
    case 'O':
        num = 0;
        break;
    case 'B':
        num = 10;
        break;
    case 'A':
        num = 20;
        break;
    case 'F':
        num = 30;
        break;
    case 'G':
        num = 40;
        break;
    case 'K':
        num = 50;
        break;
    case 'M':
        num = 60;
        break;
    case 'L':
        num = 70;
        break;
    case 'T':
        num = 80;
        break;
    case 'Y':
        num = 90;
        break;
    default:
        return -1000;
    }
 
    if( spType.size() < 2 )
        return -1000;
    if( !isdigit( spType[1] ) )
        return -1000;
 
    int i = 1;
    while( i < spType.size() && ( isdigit( spType[i] ) || spType[i] == '.' ) )
        ++i;
 
    std::string subType = spType.substr( 1, i - 1 );
 
    num += ioutils::convertFromString<realT>( subType );
 
    if( spType.size() == i )
        return num;
    if( spType[i] == 'V' )
        return num;
 
    return -1 * num;
}
 
/// Provide various characteristics of main sequence stars according to their spectral type.
/** Interpolates on "A Modern Mean Dwarf Stellar Color and Effective Temperature Sequence"
 * by Eric Mamajek, available at  http://www.pas.rochester.edu/~emamajek/EEM_dwarf_UBVIJHK_colors_Teff.txt.
 *
 * Loads values of the sequence at construction (hard coded), and then provides interpolation in between the points.
 *
 * Version of the sequence used: 2017.03.14
 *
 * \ingroup stars
 */
template <typename realT>
struct mainSequence
{
    std::vector<double> m_numTypes;  ///< The numeric spectral type codes from the sequence
    std::vector<double> m_numTypesR; ///< The numeric spectral type codes from the sequence, in reversed order.
 
    std::vector<double> m_Teffs;  ///< The effective temperatures from the sequence
    std::vector<double> m_TeffsR; ///< The effective temperatures from the sequence, reveresed order.
 
    std::vector<double> m_logLs;  ///< The log10 luminosities from the sequence
    std::vector<double> m_rads;   ///< The radii from the sequence
    std::vector<double> m_masses; ///< The masses from the sequence
    std::vector<double> m_Mvs;    ///< The absolute V magnitudes from the sequence
    std::vector<double> m_V_Rcs;  ///< The V-Rc colors from the sequence
    std::vector<double> m_V_Ics;  ///< The V-Ic colors from the sequence
    std::vector<double> m_V_Kss;  ///< The V-Ks colors from the sequence
    std::vector<double> m_J_Hs;   ///< The J-H colors from the sequence
    std::vector<double> m_H_Kss;  ///< The H-Ks colors from the sequence
 
    math::gslInterpolator<math::gsl_interp_linear<double>> interpT; ///< The interpolator for effective Temperature
    double m_minT; ///< The minimum numeric spectral type for effective Temperature
    double m_maxT; ///< The maximum numeric spectral type for effective Temperature
 
    math::gslInterpolator<math::gsl_interp_linear<double>> interpRad; ///< The interpolator for effective Temperature
    double m_minRad; ///< The minimum numeric spectral type for effective Temperature
    double m_maxRad; ///< The maximum numeric spectral type for effective Temperature
 
    math::gslInterpolator<math::gsl_interp_linear<double>> interpL; ///< The interpolator for log luminosity
    double m_minL; ///< The minimum numeric spectral type for log luminosity
    double m_maxL; ///< The maximum numeric spectral type for log luminosity
 
    math::gslInterpolator<math::gsl_interp_linear<double>> interpMv; ///< The interpolator for absolute V magnitude
    double m_minMv; ///< The minimum numeric spectral type for absolute V magnitude
    double m_maxMv; ///< The maximum numeric spectral type for absolute V magnitude
 
    math::gslInterpolator<math::gsl_interp_linear<double>> interpVRc; ///< The interpolator for V-Rc color
    double m_minVRc; ///< The minimum numeric spectral type for V-Rc color
    double m_maxVRc; ///< The maximum numeric spectral type for V-Rc color
 
    math::gslInterpolator<math::gsl_interp_linear<double>> interpVIc; ///< The interpolator for V-Ic color
    double m_minVIc; ///< The minimum numeric spectral type for V-Ic color
    double m_maxVIc; ///< The maximum numeric spectral type for V-Ic color
 
    math::gslInterpolator<math::gsl_interp_linear<double>> interpVKs; ///< The interpolator for V-Ks color
    double m_minVKs; ///< The minimum numeric spectral type for V-Ks color
    double m_maxVKs; ///< The maximum numeric spectral type for V-Ks color
 
    math::gslInterpolator<math::gsl_interp_linear<double>> interpJH; ///< The interpolator for J-H color
    double m_minJH; ///< The minimum numeric spectral type for J-H color
    double m_maxJH; ///< The maximum numeric spectral type for J-H color
 
    math::gslInterpolator<math::gsl_interp_linear<double>> interpHKs; ///< The interpolator for H-Ks color
    double m_minHKs; ///< The minimum numeric spectral type for H-Ks color
    double m_maxHKs; ///< The maximum numeric spectral type for H-Ks color
 
    math::gslInterpolator<math::gsl_interp_linear<double>> interpSpTfmT; ///< The interpolator for effective Temperature
    double m_minSpTfmT; ///< The minimum numeric spectral type for effective Temperature
    double m_maxSpTfmT; ///< The maximum numeric spectral type for effective Temperature
 
    void findMinMax( std::vector<double> &seq, double &min, double &max, std::vector<double> &ref )
    {
        if( seq.size() != ref.size() )
        {
            std::cerr << "size mismatch\n";
            min = 0;
            max = 0;
        }
 
        min = ref[0];
        size_t n;
        for( n = 0; n < seq.size(); ++n )
        {
            if( seq[n] != -99 )
                break;
        }
 
        if( n >= seq.size() - 1 )
        {
            min = 0;
            max = 0;
            return;
        }
 
        min = ref[n];
 
        for( ; n < seq.size() - 1; ++n )
        {
            if( seq[n + 1] == -99 )
                break;
        }
 
        max = ref[n];
        ++n;
        for( ; n < seq.size(); ++n )
        {
            seq[n] = -99;
        }
    }
 
    mainSequence()
    {
#include "mamajek.inc"
 
        m_numTypesR.assign( m_numTypes.rbegin(), m_numTypes.rend() );
 
        findMinMax( m_Teffs, m_minT, m_maxT, m_numTypes );
        m_TeffsR.assign( m_Teffs.rbegin(), m_Teffs.rend() ); // get reversed vector after normalizing -99s.
 
        findMinMax( m_rads, m_minRad, m_maxRad, m_numTypes );
 
        findMinMax( m_logLs, m_minL, m_maxL, m_numTypes );
        findMinMax( m_Mvs, m_minMv, m_maxMv, m_numTypes );
        findMinMax( m_V_Rcs, m_minVRc, m_maxVRc, m_numTypes );
        findMinMax( m_V_Ics, m_minVIc, m_maxVIc, m_numTypes );
        findMinMax( m_V_Kss, m_minVKs, m_maxVKs, m_numTypes );
        findMinMax( m_J_Hs, m_minJH, m_maxJH, m_numTypes );
        findMinMax( m_H_Kss, m_minHKs, m_maxHKs, m_numTypes );
 
        findMinMax( m_numTypesR, m_minSpTfmT, m_maxSpTfmT, m_TeffsR );
 
        interpT.setup( m_numTypes, m_Teffs );
        interpRad.setup( m_numTypes, m_rads );
        interpL.setup( m_numTypes, m_logLs );
        interpMv.setup( m_numTypes, m_Mvs );
        interpVRc.setup( m_numTypes, m_V_Rcs );
        interpVIc.setup( m_numTypes, m_V_Ics );
        interpVKs.setup( m_numTypes, m_V_Kss );
        interpJH.setup( m_numTypes, m_J_Hs );
        interpHKs.setup( m_numTypes, m_H_Kss );
 
        interpSpTfmT.setup( m_TeffsR, m_numTypesR );
    }
 
    /// Get the interpolated effective temperature
    /**
     * \returns the effective temperature.
     * \returns -999 if the input type code is out of range.
     */
    realT Teff( realT numType /**< [in] the numeric spectral type code, see \ref mx::astro::numSpType() */ )
    {
        if( numType < m_minT || numType > m_maxT )
        {
            return -999;
        }
 
        return interpT( numType );
    }
 
    /// Get the interpolated effective temperature
    /**
     * \returns the effective temperature.
     * \returns -999 if the input type is out of range.
     */
    realT Teff( const std::string &spType /**< [in] the spectral type in standard format */ )
    {
        return Teff( numSpType( spType ) );
    }
 
    /// Get the interpolated radius
    /**
     * \returns the radius in Solar units.
     * \returns -999 if the input type code is out of range.
     */
    realT radius( realT numType /**< [in] the numeric spectral type code, see \ref mx::astro::numSpType() */ )
    {
        if( numType < m_minRad || numType > m_maxRad )
        {
            return -999;
        }
 
        return interpRad( numType );
    }
 
    /// Get the interpolated radius
    /**
     * \returns the radius in Solar units.
     * \returns -999 if the input type is out of range.
     */
    realT radius( const std::string &spType /**< [in] the spectral type in standard format */ )
    {
        return radius( numSpType( spType ) );
    }
 
    /// Get the interpolated log of luminosity
    /**
     * \returns the log of luminosity.
     * \returns -999 if the input type code is out of range.
     */
    realT logL( realT numType /**< [in] the numeric spectral type code, see \ref mx::astro::numSpType() */ )
    {
        if( numType < m_minL || numType > m_maxL )
        {
            return -999;
        }
 
        return interpL( numType );
    }
 
    /// Get the interpolated log of luminosity
    /**
     * \returns the log of luminosity.
     * \returns -999 if the input type is out of range.
     */
    realT logL( const std::string &spType /**< [in] the spectral type in standard format */ )
    {
        return logL( numSpType( spType ) );
    }
 
    /// Get the interpolated absolute V magnitude
    /**
     * \returns the aboslute V magnitude
     * \returns -999 if the input type code is out of range.
     */
    realT Mv( realT numType /**< [in] the numeric spectral type code, see \ref mx::astro::numSpType() */ )
    {
        if( numType < m_minMv || numType > m_maxMv )
        {
            return -999;
        }
 
        return interpMv( numType );
    }
 
    /// Get the interpolated absolute V magnitude
    /**
     * \returns the aboslute V magnitude
     * \returns -999 if the input type is out of range.
     */
    realT Mv( const std::string &spType /**< [in] the spectral type in standard format */ )
    {
        return Mv( numSpType( spType ) );
    }
 
    /// Get the interpolated absolute V-Rc color
    /**
     * \returns the aboslute V-Rc color
     * \returns -999 if the input type code is out of range.
     */
    realT V_Rc( realT numType /**< [in] the numeric spectral type code, see \ref mx::astro::numSpType() */ )
    {
        if( numType < m_minVRc || numType > m_maxVRc )
        {
            return -999;
        }
 
        return interpVRc( numType );
    }
 
    /// Get the interpolated absolute V-Rc color
    /**
     * \returns the aboslute V-Rc color
     * \returns -999 if the input type is out of range.
     */
    realT V_Rc( const std::string &spType /**< [in] the spectral type in standard format */ )
    {
        return V_Rc( numSpType( spType ) );
    }
 
    /// Get the interpolated absolute V-Ic color
    /**
     * \returns the aboslute V-Ic color
     * \returns -999 if the input type code is out of range.
     */
    realT V_Ic( realT numType /**< [in] the numeric spectral type code, see \ref mx::astro::numSpType() */ )
    {
        if( numType < m_minVIc || numType > m_maxVIc )
        {
            return -999;
        }
 
        return interpVIc( numType );
    }
 
    /// Get the interpolated absolute V-Ic color
    /**
     * \returns the aboslute V-Ic color
     * \returns -999 if the input type is out of range.
     */
    realT V_Ic( const std::string &spType /**< [in] the spectral type in standard format */ )
    {
        return V_Ic( numSpType( spType ) );
    }
 
    /// Get the interpolated absolute V-Ks color
    /**
     * \returns the aboslute V-Ks color
     * \returns -999 if the input type code is out of range.
     */
    realT V_Ks( realT numType /**< [in] the numeric spectral type code, see \ref mx::astro::numSpType() */ )
    {
        if( numType < m_minVKs || numType > m_maxVKs )
        {
            return -999;
        }
 
        return interpVKs( numType );
    }
 
    /// Get the interpolated absolute V-Ks color
    /**
     * \returns the aboslute V-Ks color
     * \returns -999 if the input type is out of range.
     */
    realT V_Ks( const std::string &spType /**< [in] the spectral type in standard format */ )
    {
        return V_Ks( numSpType( spType ) );
    }
 
    /// Get the interpolated absolute J-H color
    /**
     * \returns the aboslute J-H color
     * \returns -999 if the input type code is out of range.
     */
    realT J_H( realT numType /**< [in] the numeric spectral type code, see \ref mx::astro::numSpType() */ )
    {
        if( numType < m_minJH || numType > m_maxJH )
        {
            return -999;
        }
 
        return interpJH( numType );
    }
 
    /// Get the interpolated absolute J-H color
    /**
     * \returns the aboslute J-H color
     * \returns -999 if the input type is out of range.
     */
    realT J_H( const std::string &spType /**< [in] the spectral type in standard format */ )
    {
        return J_H( numSpType( spType ) );
    }
 
    /// Get the interpolated absolute H-Ks color
    /**
     * \returns the aboslute H-Ks color
     * \returns -999 if the input type code is out of range.
     */
    realT H_Ks( realT numType /**< [in] the numeric spectral type code, see \ref mx::astro::numSpType() */ )
    {
        if( numType < m_minHKs || numType > m_maxHKs )
        {
            return -999;
        }
 
        return interpHKs( numType );
    }
 
    /// Get the interpolated absolute H-Ks color
    /**
     * \returns the aboslute H-Ks color
     * \returns -999 if the input type is out of range.
     */
    realT H_Ks( const std::string &spType /**< [in] the spectral type in standard format */ )
    {
        return H_Ks( numSpType( spType ) ); /// Get the interpolated absolute H-Ks color
        /**
         * \returns the aboslute H-Ks color
         * \returns -999 if the input type code is out of range.
         */
    }
 
    realT spTFromTeff( realT Teff )
    {
        if( Teff < m_minSpTfmT || Teff > m_maxSpTfmT )
        {
            return -999;
        }
 
        return std::floor( 2 * interpSpTfmT( Teff ) + 0.5 ) / 2; // round to nearest 0.5 types.
    }
};
 
namespace maintenance
{
 
/// Read in the main sequence table of Mamajek, and construct the vectors for input into the mainSequence class.
/** This is used to create mamajek.inc.
 * The table should be copied to a text file, and all `...` replaced with -99, then any remaining . replaced with space.
 */
void makeMSTable( const std::string &fn )
{
    std::vector<std::string> SpT;
    std::vector<double> Teff;
    std::vector<double> logT;
    std::vector<double> logL;
    std::vector<double> Mbol;
    std::vector<double> BCv;
    std::vector<double> Mv;
    std::vector<double> B__V;
    std::vector<double> Bt__Vt;
    std::vector<double> G__V;
    std::vector<double> Bp__Rp;
    std::vector<double> G__Rp;
    std::vector<double> M_G;
    std::vector<double> b__y;
    std::vector<double> U__B;
    std::vector<double> V__Rc;
    std::vector<double> V__Ic;
    std::vector<double> V__Ks;
    std::vector<double> J__H;
    std::vector<double> H__Ks;
    std::vector<double> Ks__W1;
    std::vector<double> W1__W2;
    std::vector<double> W1__W3;
    std::vector<double> W1__W4;
    std::vector<double> M_J;
    std::vector<double> M_Ks;
    std::vector<double> i__z;
    std::vector<double> z__Y;
    std::vector<double> R_Rsun;
    std::vector<double> Msun;
 
    ioutils::readColumns( fn,
                          SpT,
                          Teff,
                          logT,
                          logL,
                          Mbol,
                          BCv,
                          Mv,
                          B__V,
                          Bt__Vt,
                          G__V,
                          Bp__Rp,
                          G__Rp,
                          M_G,
                          b__y,
                          U__B,
                          V__Rc,
                          V__Ic,
                          V__Ks,
                          J__H,
                          H__Ks,
                          Ks__W1,
                          W1__W2,
                          W1__W3,
                          W1__W4,
                          M_J,
                          M_Ks,
                          i__z,
                          z__Y,
                          R_Rsun,
                          Msun );
 
    std::cout << "//" << fn << "\n";
    std::cout << "         m_numTypes = {" << numSpType( SpT[0] );
    for( int n = 1; n < SpT.size(); ++n )
        std::cout << ',' << numSpType( SpT[n] );
    std::cout << "};\n";
 
    std::cout << "         m_Teffs = {" << Teff[0];
    for( int n = 1; n < SpT.size(); ++n )
        std::cout << ',' << Teff[n];
    std::cout << "};\n";
 
    std::cout << "         m_logLs = {" << logL[0];
    for( int n = 1; n < SpT.size(); ++n )
        std::cout << ',' << logL[n];
    std::cout << "};\n";
 
    std::cout << "         m_rads = {" << R_Rsun[0];
    for( int n = 1; n < SpT.size(); ++n )
        std::cout << ',' << R_Rsun[n];
    std::cout << "};\n";
 
    std::cout << "         m_masses = {" << Msun[0];
    for( int n = 1; n < SpT.size(); ++n )
        std::cout << ',' << Msun[n];
    std::cout << "};\n";
 
    std::cout << "         m_Mvs = {" << Mv[0];
    for( int n = 1; n < SpT.size(); ++n )
        std::cout << ',' << Mv[n];
    std::cout << "};\n";
 
    std::cout << "         m_V_Rcs = {" << V__Rc[0];
    for( int n = 1; n < SpT.size(); ++n )
        std::cout << ',' << V__Rc[n];
    std::cout << "};\n";
 
    std::cout << "         m_V_Ics = {" << V__Ic[0];
    for( int n = 1; n < SpT.size(); ++n )
        std::cout << ',' << V__Ic[n];
    std::cout << "};\n";
 
    std::cout << "         m_V_Kss = {" << V__Ks[0];
    for( int n = 1; n < SpT.size(); ++n )
        std::cout << ',' << V__Ks[n];
    std::cout << "};\n";
 
    std::cout << "         m_J_Hs = {" << J__H[0];
    for( int n = 1; n < SpT.size(); ++n )
        std::cout << ',' << J__H[n];
    std::cout << "};\n";
 
    std::cout << "         m_H_Kss = {" << H__Ks[0];
    for( int n = 1; n < SpT.size(); ++n )
        std::cout << ',' << H__Ks[n];
    std::cout << "};\n";
}
} // namespace maintenance
} // namespace astro
} // namespace mx
 
#endif