OrbitSolver.h

Go to the documentation of this file.

#ifndef __ORBIT_SOLVER_H
#define __ORBIT_SOLVER_H

#include "../Core/SharedLibraryExportMacros.h"
#include "../Core/AstronomicalConstants.h"
#include "../Core/Common.h"
#include "../Core/OrbitalExpressions.h"
#include "../Core/Common.h"
#include "BinarySystem.h"
#include "CombinedStoppingCondition.h"
#include "ExternalStoppingConditions.h"
#include "StopInformation.h"
#include "StopHistoryInterval.h"
#include <math.h>
#include <list>
#include <vector>
#include <stdlib.h>
#include <fstream>
#include <iostream>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_poly.h>
#include <sstream>
#include <iostream>
#include <limits>
#include <ctime>

namespace Evolve {

    typedef int (*GSL_ODE_TYPE)(double, const double*, double*, void*);
    typedef int (*GSL_JAC_TYPE)(double, const double*, double*, double*, void*);
    typedef bool (*STOP_EVOL_TYPE)(double, const double*, void*);

    LIB_LOCAL int stellar_system_diff_eq(
            double age,

            const double *parameters,

            double *derivatives,

            void *system);

    LIB_LOCAL int stellar_system_jacobian(
            double age,

            const double *parameters,

            double *param_derivs,

            double *age_derivs,

            void *system_mode);

    class LIB_LOCAL ExtremumInformation {
    private:
        double __x,

               __y;
    public:
        ExtremumInformation(
                double x=Core::Inf,

                double y=Core::NaN) : __x(x), __y(y) {}

        double x() const {return __x;}
        double &x() {return __x;}

        double y() const {return __y;}
        double &y() {return __y;}
    };//End ExtremumInformation class.

    class LIB_PUBLIC OrbitSolver {
    private:
        double
            __end_age,

            __precision,

            __e_order_downgrade_threshold,

            __last_e_order_upgrade_age;

        std::list<double> __tabulated_ages;

        std::list<Core::EvolModeType> __tabulated_evolution_modes;

        std::valarray<size_t> __skip_history_zerocrossing;

        std::valarray<double> __skip_history_extremum;

        std::list<double> __stop_history_ages,

            __discarded_stop_ages;

        std::list< std::valarray<double> >
            __orbit_history,
            __orbit_deriv_history, 
            __stop_cond_history,

            __stop_deriv_history,

            __stop_cond_discarded,

            __stop_deriv_discarded;

        StoppingCondition *__stopping_conditions;

        std::vector<StopInformation> __stop_info;

        bool __print_progress;

        time_t __evolution_start_time;

        double __runtime_limit;

#ifndef NDEBUG
        void output_history_and_discarded(std::ostream &os);
#endif

        void clear_discarded();

        void insert_discarded(double age,
                const std::valarray<double> &current_stop_cond,
                const std::valarray<double> &current_stop_deriv);

        void add_to_evolution(
                double age,

                Core::EvolModeType evolution_mode,

                BinarySystem &system);

        double go_back(double max_age,
                       BinarySystem &system,
                       std::valarray<double> &orbit);


        void clear_history();

        StopHistoryInterval select_stop_condition_interval(bool crossing,
                size_t cond_ind, size_t max_points) const;

        ExtremumInformation extremum_from_history_no_deriv(
                size_t condition_index) const;

        ExtremumInformation extremum_from_history(size_t condition_index) const;

        double crossing_from_history_no_deriv(
                size_t condition_index) const;

        double crossing_from_history(
                size_t condition_index) const;

        void initialize_skip_history(const StoppingCondition &stop_cond,
                StoppingConditionType stop_reason);

        bool at_exact_condition(double previous_age,
                                const StopInformation &stop_info);

    /*  void update_skip_history(
                const std::valarray<double> &current_stop_cond,
                const StopInformation &stop_info);*/

        bool acceptable_step(double current_age,
                             double previous_age,
                             const StopInformation &stop_info);

        //(e-order should be increased).
        int check_expansion_error(
            const std::valarray<double> &derivatives,

            const std::valarray<double> &expansion_errors
        );

        StopInformation update_stop_condition_history(
            double age,

            const std::valarray<double> &orbit,

            const std::valarray<double> &derivatives,

            const std::valarray<double> &expansion_errors,

            Core::EvolModeType evolution_mode,

            StoppingConditionType stop_reason=NO_STOP,

            bool ignore_e_order_decrease=false
        );

        void reject_step(
            double &age,

            StopInformation &stop,

            BinarySystem &system,

            std::valarray<double> &orbit,

            double &max_next_t,

            double &step_size
#ifndef NDEBUG
            , std::string reason
#endif
        );

        StopInformation evolve_until(
            BinarySystem &system,

            double &max_age,

            std::valarray<double> &orbit,

            StoppingConditionType &stop_reason,

            double max_step,

            Core::EvolModeType evolution_mode
        );

        CombinedStoppingCondition *get_stopping_condition(
            BinarySystem &system
        );

        double stopping_age(
                double age,

                const BinarySystem &system,

                const std::list<double> &required_ages);

        void reached_stopping_condition(
            double stop_age,

            StoppingConditionType stop_reason
        );

        void adjust_eccentricity_order(
            BinarySystem &system,

            const std::valarray<double> &orbit,

            Core::EvolModeType evolution_mode,

            bool must_increase = false
        );


        void reset(BinarySystem &system);
    public:
        OrbitSolver(
            double max_age,

            double required_precision,

            bool print_progress=false
        );

        void operator()(
            BinarySystem &system,

            double max_step=Core::Inf,

            const std::list<double> &required_ages=std::list<double>(),

            double max_runtime=0
        );

        const std::list<double> &evolution_ages() const
        {return __tabulated_ages;}

        const std::list<Core::EvolModeType> &mode_evolution() const
        {return __tabulated_evolution_modes;}

        ~OrbitSolver()
        {if(__stopping_conditions) delete __stopping_conditions;}

    }; //End OrbitSolver class.

}//End Evolve namespace.

#endif