EvolvingStellarQuantity.h

Go to the documentation of this file.

#ifndef __EVOLVING_STELLAR_QUANTITY_H
#define __EVOLVING_STELLAR_QUANTITY_H

#include "../Core/SharedLibraryExportMacros.h"
#include "AllowedGridGrowth.h"
#include "RemoveLogDeriv.h"
#include "SumDerivatives.h"
#include "InterpolatedDerivatives.h"
#include "mass_feh_interp.h"
#include "../Core/Functions.h"
#include "../Core/InterpSolutionIterator.h"
#include "../Core/InterpolatingFunctionALGLIB.h"
#include "../Core/Error.h"
#include <valarray>
#include <list>
#include <string>
#include <iostream>
#include <fstream>
#include <sstream>

namespace StellarEvolution {

    using Core::OneArgumentDiffFunction;
    using Core::FunctionDerivatives;
    using Core::InterpSolutionIterator;

    class LIB_PUBLIC EvolvingStellarQuantity : public OneArgumentDiffFunction {
    private:
        double 
            __mass,

            __feh;


        double 
            __min_age,

            __max_age; 

        bool __log_age;

        bool __log_quantity;

        bool __initially_zero;

        std::valarray<double> 
            __track_masses,

            __track_feh;

        std::valarray<double> 
            __min_interp_ages,

            __max_interp_ages;

        std::vector<double> __interp_grid_change_ages;

        mutable std::vector<double>::const_iterator __next_grid_change_age;

        mutable size_t 
            __mass_index_above,

            __mass_index_below,

            __feh_index_above,

            __feh_index_below;

        std::vector<const OneArgumentDiffFunction *> __evolution_tracks;

        mutable alglib::real_1d_array
            __interp_masses,

            __interp_feh;

        mutable size_t 
            __min_interp_mass_index,

            __max_interp_mass_index,

            __min_interp_feh_index,

            __max_interp_feh_index;

        inline size_t track_index(
            size_t mass_index,

            size_t feh_index
        ) const
        {return feh_index * __track_masses.size() + mass_index;}

        inline bool track_in_range(
            size_t track_i,

            double age) const
        {
            assert(track_i < __min_interp_ages.size());
            assert(track_i < __max_interp_ages.size());
            return (__min_interp_ages[track_i] <= age
                    &&
                    __max_interp_ages[track_i] >= age);
        }

        inline bool track_in_range(
            size_t mass_i,

            size_t feh_i,

            double age) const
        {return track_in_range(track_index(mass_i, feh_i), age);}

        void check_grid_range() const;

        void find_cell(
            const std::valarray<double> &boundaries,

            double value,

            size_t &below_index,

            size_t &above_index
        ) const;

        void set_interp_age_ranges();

        double evaluate_track(
            double age,

            const OneArgumentDiffFunction &track,

            const FunctionDerivatives **derivatives) const;

        void check_grid_expansion_directions(
            AllowedGridGrowth &grow,

            double age
        ) const;

        bool expand_grid(
            const AllowedGridGrowth &grow,

            double age
        ) const;

        void update_interpolation_grid() const;

        double interpolate(
            double age, 
            const FunctionDerivatives **derivatives=NULL
        ) const;

        inline double age_to_interp_param(double age) const
        {return age_to_interp_param(age, __mass, __feh);}

        inline double interp_param_to_age(double interp_param) const
        {return interp_param_to_age(interp_param, __mass, __feh);}

    protected:
        virtual double age_to_interp_param(
            double age,

            double mass,

            double feh
        ) const;

        virtual double interp_param_to_age(
            double interp_param,

            double mass,

            double feh
        ) const;

    public:
        EvolvingStellarQuantity() {}

        EvolvingStellarQuantity(
            double mass, 

            double feh, 

            const std::valarray<double> &track_masses,

            const std::valarray<double> &track_feh,

            //mass index varies faster.
            const std::vector<const OneArgumentDiffFunction *> 
            &evolution_tracks,

            bool log_age=true,

            bool log_quantity=true,

            bool starts_zero=false
        );

        virtual void select_interpolation_region(double age) const;

        virtual double operator()(double age) const
        {return interpolate(age);}

        virtual const FunctionDerivatives *deriv(double age) const;

        virtual double range_high() const {return __max_age;}

        virtual double range_low() const 
        {return (__initially_zero ? -Core::Inf : __min_age);}

        virtual const std::vector<double> &discontinuities() const
        {return __interp_grid_change_ages;}

        virtual double next_discontinuity() const
        {return *__next_grid_change_age;}

        virtual double previous_discontinuity() const;

        virtual void enable_next_interpolation_region() const;

        InterpSolutionIterator crossings(double =0) const
        {
            throw Core::Error::Runtime(
                "Called EvolvingStellarQuantity::crossings, which are ill"
                "defined."
            );
        }
    }; //End EvolvingStellarQuantity class declaration.

} //End StellarEvolution namespace

#endif