PolynomialEvolution.h

#ifndef __POLYNOMIAL_EVOLUTION_H
#define __POLYNOMIAL_EVOLUTION_H

#include "../Core/Functions.h"
#include "../StellarEvolution/EvolvingStellarQuantity.h"
#include "../StellarEvolution/Interpolator.h"

#include "Common.h"
//#include "../StellarSystem.h"
//#include "../Star.h"
//#include "../Planet.h"
//#include "../OrbitSolver.h"
#include <valarray>
#include <gsl/gsl_roots.h>
//#include "YRECIO.h"
#include <fstream>
#include <boost/archive/text_oarchive.hpp>
#include <boost/archive/text_iarchive.hpp>

namespace StellarEvolution {

    class PolynomialEvolutionQuantity : public EvolvingStellarQuantity,
        Core::FunctionDerivatives {
    private:
        std::valarray<double> __poly_coef;

        double 
            __xmin,

            __xmax;

        mutable double __deriv_x;

        std::vector<double> __empty_vector;

    public:
        PolynomialEvolutionQuantity(
            const std::valarray<double> &coefficients,

            double range_low,

            double range_high,

            double derivative_x=NaN
        ) :
            __poly_coef(coefficients),
            __xmin(range_low),
            __xmax(range_high),
            __deriv_x(derivative_x)
            {}

        void select_interpolation_region(double) const {}

        double operator()(double x) const {__deriv_x = x; return order(0);}

        double range_low() const {return __xmin;}

        double range_high() const {return __xmax;}

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

        virtual double previous_discontinuity() const {return __xmin;}

        virtual double next_discontinuity() const {return __xmax;}

        virtual void enable_next_interpolation_region() const
        {assert(false);}

        const Core::FunctionDerivatives *deriv(double x) const
        {return new PolynomialEvolutionQuantity(__poly_coef,
                                                __xmin,
                                                __xmax,
                                                x);}

        double order(unsigned deriv_order = 1) const;

        friend std::ostream &operator<<(
            std::ostream &os,
            const PolynomialEvolutionQuantity &track
        );

        std::string kind() const {return "PolynomialEvolutionQuantity";}
    };//End PolynomialEvolutionQuantity class.

    class MockStellarEvolution : public StellarEvolution::Interpolator {
    private:
        std::valarray< std::valarray<double> >
            __R,

            __Iconv,

            __Irad,

            __Itot,

            __Rcore,

            __Mcore,

            __Menv,
            
            __Lum;

        double __core_formation_age;
    public:
        MockStellarEvolution(
            double core_formation_age = Core::NaN,

            const std::valarray< std::valarray<double> > &
                R = std::valarray< std::valarray<double> >(),

            const std::valarray< std::valarray<double> > &
                Iconv = std::valarray< std::valarray<double> >(),

            const std::valarray< std::valarray<double> > &
                Irad = std::valarray< std::valarray<double> >(),

            const std::valarray< std::valarray<double> > &
                Rcore = std::valarray< std::valarray<double> >(),

            const std::valarray< std::valarray<double> > &
                Mcore = std::valarray< std::valarray<double> >(),

            const std::valarray< std::valarray<double> > &
                Lum = std::valarray< std::valarray<double> >()
        );

        EvolvingStellarQuantity *operator()(QuantityID quantity,
                                            double mass,
                                            double feh) const;

        double core_formation_age() const {return __core_formation_age;}
    };

    MockStellarEvolution *make_no_evolution(double Rstar = 1.0,
                                            double Iconv = 1.0);

    MockStellarEvolution *make_linear_I_evolution();

}//End StellarEvolution namespace.

#if 0
    class StarData {
    private:
        unsigned __num_stars_created;
    public:
        double mass,
               feh,
               radius,
               age,
               conv_spin,
               rad_spin,
               tidal_Q,
               wind_strength,
               wind_sat_freq,
               coupling_timescale,
               Q_trans_width,
               disk_lock_w,
               disk_lock_time;

        PolynomialEvolutionTrack *Lrad_track,
                                 *Lconv_track;

        std::valarray<double> evolution_masses,
                              evolution_ages,
                              Lrad,
                              Lconv,
                              Iconv,
                              Irad,
                              Mrad_deriv,
                              Lconv_deriv,
                              Lrad_deriv,
                              Rrad,
                              all_radii;

        std::valarray< std::valarray<double> > r_coef,
                                               Iconv_coef,
                                               Irad_coef,
                                               Mrad_coef,
                                               Rcore_coef;

        StarData();

        void create_random_star(Star::EvolvingStar **star);

        ~StarData();
    };


class PlanetData {
public:
    Star* star;
    StarData* sdata;
    double mass, radius;
    std::valarray<double> ages;
    std::valarray<double> semis;
    PlanetData() : star(NULL) {
        sdata = new StarData();
        sdata->create_random_star(&star);
    };
    void create_random_planet(Planet** planet);
    double get_semi(double age);

    ~PlanetData();
};

class SystemData {
public:
    Star* star;
    Planet* planet;
    StarData* sdata;
    PlanetData* pdata;
    SystemData() {
        //exit(-1);
        pdata = new PlanetData();
        pdata->create_random_planet(&planet);
        star = pdata->star;
        sdata = pdata->sdata;
    }
    void create_random_system(StellarSystem** system) {
/*      double min_age = sdata->evolution_ages[0];
        double precision = 0.01;
        double max_age = star->get_lifetime();
        OrbitSolver orb(min_age, max_age, precision,
                stellar_system_diff_eq, stellar_system_jacobian);*/
        *system = new StellarSystem(*star, *planet, "");
    }
};
#endif

namespace StellarEvolution {

    PolynomialEvolutionQuantity *exact_track(
        const std::valarray< std::valarray<double> > &poly_coef, 

        double mass,

        double low_mass_age_scaling = 0,

        double high_mass_age_scaling = 0,

        double scale_mass = NaN
    );

} //End StellarEvolution namespace.

#if 0
class PolynomialStellarEvolution : public StellarEvolution {
public:
    PolynomialStellarEvolution(
            const std::valarray<double> &masses,

            const std::valarray<double> &ages,

            const std::valarray< std::valarray<double> > &r_coef,

            const std::valarray< std::valarray<double> > &Iconv_coef,

            const std::valarray< std::valarray<double> > &Itot_coef,

            const std::valarray< std::valarray<double> > &Mrad_coef,

            const std::valarray< std::valarray<double> > &Rcore_coef,

            double low_mass_age_scaling,
            
            double high_mass_age_scaling);
};
#endif

class ExponentialPlusFunc : public Core::OneArgumentDiffFunction,
                                   Core::FunctionDerivatives
{
private:
    double __scale, __rate, __deriv_x;
    Core::OneArgumentDiffFunction *__offset;
public:
    ExponentialPlusFunc(Core::OneArgumentDiffFunction *offset,
                        double scale,
                        double rate,
                        double deriv_x = Core::NaN) :
        __scale(scale),
        __rate(rate),
        __deriv_x(deriv_x),
        __offset(offset)
    {}

    double operator()(double x) const
    {return (*__offset)(x) + __scale * std::exp(__rate * x);}

    const Core::FunctionDerivatives *deriv(double x) const
    {return new ExponentialPlusFunc(__offset, __scale, __rate, x);}

    double order(unsigned deriv_order = 1) const;

    double range_high() const {return __offset->range_high();}

    double range_low() const {return __offset->range_low();}

    Core::InterpSolutionIterator crossings(double) const
    {return Core::InterpSolutionIterator();}
};

class FuncPlusFunc : public Core::OneArgumentDiffFunction,
                            Core::FunctionDerivatives {
private:
    const OneArgumentDiffFunction *__f1, *__f2;
    double __deriv_x;
public:
    FuncPlusFunc(const OneArgumentDiffFunction *f1,
                 const OneArgumentDiffFunction *f2,
                 double deriv_x = Core::NaN) :
        __f1(f1), __f2(f2), __deriv_x(deriv_x)
    {}

    double operator()(double x) const {return (*__f1)(x) + (*__f2)(x);}

    const FunctionDerivatives *deriv(double x) const
    {return new FuncPlusFunc(__f1, __f2, x);}

    double order(unsigned deriv_order = 1) const;

    double range_high() const
    {return std::min(__f1->range_high(), __f2->range_high());}

    double range_low() const
    {return std::max(__f1->range_low(), __f2->range_low());}

    Core::InterpSolutionIterator crossings(double) const
    {return Core::InterpSolutionIterator();}
};

class PiecewiseFunction : public Core::OneArgumentDiffFunction,
                                 Core::FunctionDerivatives{
private:
    double __deriv_x, __range_low, __range_high;
    std::list<const OneArgumentDiffFunction *> __pieces;
public:
    PiecewiseFunction(const std::list<const OneArgumentDiffFunction *>
                      &pieces = std::list<const OneArgumentDiffFunction *>(),
                      double deriv_x = Core::NaN);

    void add_piece(const OneArgumentDiffFunction *piece);

    double operator()(double x) const;

    const FunctionDerivatives *deriv(double x) const
    {return new PiecewiseFunction(__pieces, x);}

    double order(unsigned deriv_order=1) const;

    double range_high() const {return __range_high;}

    double range_low() const {return __range_low;}

    Core::InterpSolutionIterator crossings(double) const
    {return Core::InterpSolutionIterator();}
};

class FunctionRatio: public Core::OneArgumentDiffFunction,
                            Core::FunctionDerivatives {
private:
    double __deriv_x;
    const OneArgumentDiffFunction *__f1, *__f2;
public:
    FunctionRatio(const OneArgumentDiffFunction *f1,
                  const OneArgumentDiffFunction *f2,
                  double deriv_x = Core::NaN) :
        __deriv_x(deriv_x), __f1(f1), __f2(f2) {}

    double operator()(double x) const {return (*__f1)(x) / (*__f2)(x);}

    const FunctionDerivatives *deriv(double x) const
    {return new FunctionRatio(__f1, __f2, x);}

    double order(unsigned deriv_order=1) const;

    double range_high() const
    {return std::min(__f1->range_high(), __f2->range_high());}

    double range_low() const
    {return std::max(__f1->range_low(), __f2->range_low());}

    Core::InterpSolutionIterator crossings(double) const
    {return Core::InterpSolutionIterator();}
};

class FunctionToPower : public Core::OneArgumentDiffFunction,
                               Core::FunctionDerivatives {
private:
    const OneArgumentDiffFunction *__f;
    double __power, __deriv_x;
public:
    FunctionToPower(const OneArgumentDiffFunction *f,
                    double power,
                    double deriv_x = Core::NaN) :
        __f(f), __power(power), __deriv_x(deriv_x) {}

    double operator()(double x) const
    {return std::pow((*__f)(x), __power);}

    const FunctionDerivatives *deriv(double x) const
    {return new FunctionToPower(__f, __power, x);}

    double order(unsigned deriv_order=1) const;

    double range_high() const {return __f->range_high();}

    double range_low() const {return __f->range_low();}

    Core::InterpSolutionIterator crossings(double) const
    {return Core::InterpSolutionIterator();}
};

class ScaledFunction : public Core::OneArgumentDiffFunction,
                              Core::FunctionDerivatives {
private:
    const OneArgumentDiffFunction *__f;
    double __scale, __deriv_x;
public:
    ScaledFunction(const OneArgumentDiffFunction *f,
                   double scale,
                   double deriv_x = Core::NaN) :
        __f(f), __scale(scale), __deriv_x(deriv_x) {}

    double operator()(double x) const
    {return __scale*(*__f)(x);}

    const FunctionDerivatives *deriv(double x) const
    {return new ScaledFunction(__f, __scale, x);}

    double order(unsigned deriv_order=1) const;

    double range_high() const {return __f->range_high();}

    double range_low() const {return __f->range_low();}

    Core::InterpSolutionIterator crossings(double) const
    {return Core::InterpSolutionIterator();}
};

class LogFunction : public Core::OneArgumentDiffFunction,
                           Core::FunctionDerivatives {
private:
    const OneArgumentDiffFunction *__f;
    double __deriv_x;
public: 
    LogFunction(const OneArgumentDiffFunction *f,
                double deriv_x = Core::NaN) :
        __f(f), __deriv_x(deriv_x) {}

    double operator()(double x) const {return std::log((*__f)(x));}

    const FunctionDerivatives *deriv(double x) const
    {return new LogFunction(__f, x);}

    double order(unsigned deriv_order=1) const;

    double range_high() const {return __f->range_high();}

    double range_low() const {return __f->range_low();}

    Core::InterpSolutionIterator crossings(double) const
    {return Core::InterpSolutionIterator();}
};

class CosFunction : public Core::OneArgumentDiffFunction,
                           Core::FunctionDerivatives {
private:
    const OneArgumentDiffFunction *__f;
    double __deriv_x;
public: 
    CosFunction(const OneArgumentDiffFunction *f,
                double deriv_x = Core::NaN) :
        __f(f), __deriv_x(deriv_x) {}

    double operator()(double x) const {return std::cos((*__f)(x));}

    const FunctionDerivatives *deriv(double x) const
    {return new CosFunction(__f, x);}

    double order(unsigned deriv_order=1) const;

    double range_high() const {return __f->range_high();}

    double range_low() const {return __f->range_low();}

    Core::InterpSolutionIterator crossings(double) const
    {return Core::InterpSolutionIterator();}
};


#if 0
std::valarray<double> tabulate_track(PolynomialEvolutionTrack *track,
        std::valarray<double> ages, unsigned deriv_order=0);

double eval_poly(const std::valarray< std::valarray<double> > &poly_coef,
                 double mass, double age, double low_mass_age_scaling=0,
        double high_mass_age_scaling=0, double scale_mass=NaN);
#endif

double solve(double guess_x,
             double abs_precision,
             double rel_precision,
             double (*f)(double x, void *params),
             double (*df) (double x, void *params),
             void (*fdf) (double x, void *params, double *f, double *df),
             void *params);

#endif