phoenixSpectrum.hpp

/** \file phoenixSpectra.hpp
 * \author Jared R. Males
 * \brief Utilities for working with spectra from phoenix code.
 * \ingroup astrophot
 *
 */
 
#ifndef phoenixSpectrum_hpp
#define phoenixSpectrum_hpp
 
#include "../math/vectorUtils.hpp"
#include "../ioutils/fileUtils.hpp"
 
namespace mx
{
namespace astro
{
 
/// A spectrum from the Phoenix model, for use with ioutils::astro::astroSpectrum.
/** \ingroup astrophot_spectra
 */
template <typename _units>
struct phoenixSpectrum
{
    typedef _units units;
    typedef typename units::realT realT;
 
    static const bool freq = false;
 
    /// Convert from A to SI m
    static constexpr realT wavelengthUnits = static_cast<realT>( 1e10 );
 
    /// Convert from erg s-1 cm-2 A-1 to SI W m-3
    static constexpr realT fluxUnits =
        static_cast<realT>( 1e7 ) / ( static_cast<realT>( 1e4 ) * static_cast<realT>( 1e10 ) );
 
    static constexpr const char *dataDirEnvVar = "PHOENIX_DATADIR";
 
    typedef std::string paramsT; ///< The parameter is a string name.
 
    static std::string fileName( const std::string &name )
    {
        return name;
    }
 
    static int
    readSpectrum( std::vector<realT> &rawLambda,
                  std::vector<realT> &rawSpectrum,
                  const std::string &path,
                  const paramsT &params ///< [in] the parameters are passed in case needed to construct the spectrum
    )
    {
        if( ioutils::readColumns( path, rawSpectrum ) < 0 )
            return -1;
 
        if( ioutils::readColumns( ioutils::parentPath( path ) + "/wavelength.dat", rawLambda ) < 0 )
            return -1;
 
        return 0;
    }
 
    static void scaleSpectrum( std::vector<realT> &spectrum, realT radius, realT distance )
    {
        for( int i = 0; i < spectrum.size(); ++i )
        {
            spectrum[i] *=
                pow( radius / distance * ( constants::radJupiter<units>() / constants::parsec<units>() ), 2 );
        }
    }
};
 
/// Read in, crop, scale, and re-write a Phoenix spectrum data file.
/** For working with spectra obtained from: http://perso.ens-lyon.fr/france.allard/
 *
 * The Phoenix spectra contain many columns of line data which are not often used, are
 * not necessarily sorted by wavelength, and are sometimes formated so that points with leading
 * minus signs are not in a separate column.  We also have to change the fortan D to e.
 *
 * We also want to apply the dilution factor and take the power of 10.
 *
 * This function deals with these issues, producing a two-column space-delimited
 * file in a specified wavelength range.
 *
 * References:
 * - https://phoenix.ens-lyon.fr/Grids/BT-Settl/README
 * - https://phoenix.ens-lyon.fr/Grids/FORMAT
 *
 * NOTE: This overwrites the input file!
 *
 *
 * \tparam floatT the floating point type in which to work
 *
 * \ingroup astrophot_spectra
 */
template <typename floatT>
void rewritePhoenixSpectrum(
    const std::string &filename, ///< [in.out] complete name of the file to rewrite
    floatT lmin,                 ///< [in]  minimum wavelength to rewrite [microns]
    floatT lmax,                 ///< [in] maximum wavelemngth to rewrite [microns]
    int sepWavelength =
        0, ///< [in] [optional] controls how wavelength is handled.  0=> wavelength is included in file as first column.
           ///< 1=> wavelength is written to a separate 'wavelength.dat' file. -1=> means don't write wavelength.
    floatT DF = -8.0 ///< [in] [optional] the dilution factor.  See references.
)
{
    float lambda, flambda;
 
    std::ifstream fin;
 
    fin.open( filename );
 
    if( !fin.good() )
    {
        std::cerr << "Error opening file: " << filename << "\n";
        return;
    }
 
    std::vector<floatT> lambdas;
    std::vector<floatT> flambdas;
 
    int lineSize = 1024;
    char line[lineSize];
    std::string sline;
 
    fin.getline( line, lineSize );
    sline = line;
 
    if( sline.length() < 25 )
    {
        std::cerr << "Error reading file: " << filename << "\n";
    }
 
    lambda = ioutils::convertFromString<floatT>( sline.substr( 1, 12 ) );
 
    // First have to diagnose if this has the column problem
    int nst;
    if( sline[13] == '-' )
        nst = 13;
    else
        nst = 14;
 
    // convert the D to e
    size_t dpos = sline.find( 'D', nst );
    sline[dpos] = 'e';
 
    flambda = ioutils::convertFromString<floatT>( sline.substr( nst, 12 ) );
 
    if( lambda >= lmin * 1e4 && lambda <= lmax * 1e4 )
    {
        lambdas.push_back( lambda );
        flambdas.push_back( flambda );
    }
 
    fin.getline( line, lineSize );
    sline = line;
 
    while( fin.good() )
    {
        if( sline.length() < 25 )
            continue;
 
        dpos = sline.find( 'D', nst );
        sline[dpos] = 'e';
 
        lambda = ioutils::convertFromString<floatT>( sline.substr( 1, 12 ) );
        flambda = ioutils::convertFromString<floatT>( sline.substr( nst, 12 ) );
 
        if( lambda >= lmin * 1e4 && lambda <= lmax * 1e4 )
        {
            lambdas.push_back( lambda );
            flambdas.push_back( flambda );
        }
 
        fin.getline( line, lineSize );
        sline = line;
    }
 
    fin.close();
 
    std::vector<size_t> idx = math::vectorSortOrder( lambdas );
 
    std::ofstream fout;
 
    std::string fname = filename; // + ".t";
    fout.open( filename );
    fout.precision( 10 );
 
    for( int i = 0; i < lambdas.size(); ++i )
    {
        if( sepWavelength == 0 )
        {
            fout << lambdas[idx[i]] << "    ";
        }
        fout << pow( 10, flambdas[idx[i]] + DF ) << "\n";
    }
    fout.close();
 
    if( sepWavelength == 1 )
    {
        fname = ioutils::parentPath( filename ) + "/wavelength.dat";
 
        fout.open( fname );
        fout.precision( 10 );
        for( int i = 0; i < lambdas.size(); ++i )
        {
            fout << lambdas[idx[i]] << "\n";
        }
        fout.close();
    }
 
} // rewritePhoenixSpectrum
 
/// Call rewritePhoenixSpectrum for all files in a directory.
/** \ingroup astropho_spectra
 */
template <typename floatT>
void rewritePhoenixSpectrumBatch( const std::string &dir, floatT lmin, floatT lmax, floatT DF = -8.0 )
{
    std::vector<std::string> flist = ioutils::getFileNames( dir, "lte", "", ".7" );
 
    int sepWavelength = 1;
    for( int i = 0; i < flist.size(); ++i )
    {
        rewritePhoenixSpectrum( flist[i], lmin, lmax, sepWavelength, DF );
        if( i == 0 )
            sepWavelength = -1;
    }
}
 
} // namespace astro
} // namespace mx
 
#endif // phoenixSpectrum_hpp