BinarySystem.h

Go to the documentation of this file.

#ifndef __BINARY_SYSTEM_H
#define __BINARY_SYSTEM_H

#include "../Core/SharedLibraryExportMacros.h"
#include "DissipatingBody.h"
#include "CombinedStoppingCondition.h"
#include "SecondaryDeathCondition.h"
#include "../Core/AstronomicalConstants.h"
#include "../Core/Common.h"
#include "../Core/OrbitalExpressions.h"
#include "../Core/Error.h"
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_siman.h>
#include <string>
#include <limits>
#include <iostream>

namespace Evolve {

    class LIB_PUBLIC BinarySystem {
    private:
        std::string __name;

        std::list<double>
            __semimajor_evolution,

            __eccentricity_evolution,

            //currently recorded step.
            __semimajor_rate_evolution,

            //currently recorded step.
            __eccentricity_rate_evolution;

        double
            __age,

            __semimajor,

            __eccentricity,

            __orbital_energy,

            __orbital_angmom,

            __orbit_power,

            __orbit_power_expansion_error,

            __orbit_angmom_gain,

            __orbit_angmom_gain_expansion_error;

        mutable double
            __semimajor_rate,

            __eccentricity_rate;

        Eigen::Vector3d
            __orbit_torque,

            __orbit_torque_expansion_error;

        Core::EvolModeType __evolution_mode;

        std::list<unsigned> __locked_zones;

        std::valarray<Eigen::VectorXd> __above_lock_fractions,

            __above_lock_fractions_inclination_deriv,

            __above_lock_fractions_periapsis_deriv,

            __above_lock_fractions_inertia_deriv,

            __above_lock_fractions_angmom_deriv;

        Eigen::VectorXd __above_lock_fractions_body2_radius_deriv;

        Eigen::ColPivHouseholderQR<Eigen::MatrixXd>
            __above_lock_fractions_decomp;

        DissipatingBody
            &__body1,

            &__body2;

        void find_locked_zones();

        int locked_surface_differential_equations(
            double *evolution_rates,

            bool expansion_error
        ) const;

        void locked_surface_jacobian(
            double *param_derivs,

            double *age_derivs
        ) const;

        int single_body_differential_equations(
            double *evolution_rates,

            bool expansion_error
        ) const;

        void fill_single_body_jacobian(
            double *inclination_param_derivs,

            double *periapsis_param_derivs,

            double *angmom_param_derivs,

            double *inclination_age_derivs,

            double *periapsis_age_derivs,

            double *angmom_age_derivs
        ) const;

        void single_body_jacobian(
            double *param_derivs,

            double *age_derivs
        ) const;

        double semimajor_evolution(
            double orbit_power,

            double orbit_power_deriv=Core::NaN
        ) const;

        double semimajor_evolution_expansion_error() const
        {return semimajor_evolution(__orbit_power_expansion_error);}

        double eccentricity_evolution(
            double orbit_power,

            double orbit_angmom_gain,

            //eccentricity is returned instead of the rate itself. In this
            double orbit_power_deriv=Core::NaN,

            double orbit_angmom_gain_deriv=Core::NaN,

            bool semimajor_deriv=true
        ) const;

        double eccentricity_evolution_expansion_error() const;

        void above_lock_problem_deriv_correction(
            Dissipation::QuantityEntry entry,

            bool body1_deriv,

            Eigen::MatrixXd &matrix,

            Eigen::VectorXd &rhs
        ) const;

        void calculate_above_lock_fractions(
            Eigen::VectorXd &fractions,

            Dissipation::QuantityEntry entry=Dissipation::NO_DERIV,

            bool body1_deriv=true
        );

        Eigen::VectorXd above_lock_fractions_deriv(
            Dissipation::QuantityEntry entry,

            DissipatingBody &body,

            unsigned zone_index
        );

        void fill_above_lock_fractions_deriv();

        void update_above_lock_fractions();

        void fill_orbit_torque_and_power();


        int binary_differential_equations(
            double *differential_equations,

            bool expansion_error
        ) const;


        template<typename VALUE_TYPE>
        void add_body_rate_deriv(
            const DissipatingBody &body,

            VALUE_TYPE (DissipatingBody::*func)(Dissipation::QuantityEntry,
                                                unsigned,
                                                const Eigen::VectorXd &) const,

            std::valarray<VALUE_TYPE> &orbit_rate_deriv,

            unsigned offset
        ) const;

        void fill_orbit_power_deriv(
            std::valarray<double> &orbit_power_deriv
        ) const;

        void fill_orbit_angmom_gain_deriv(
            std::valarray<double> &orbit_angmom_gain_deriv
        ) const;

        void semimajor_jacobian(
            const std::valarray<double> &orbit_power_deriv,

            bool a6p5,

            double *param_derivs,

            double &age_deriv
        ) const;

        void eccentricity_jacobian(
            const std::valarray<double> &orbit_power_deriv,

            const std::valarray<double> &orbit_angmom_gain_deriv,

            bool a6p5,

            double *param_derivs,

            double &age_deriv
        ) const;

        void angle_evolution_age_deriv(
            DissipatingBody &body,

            unsigned zone_ind,

            double sin_inc,

            double cos_inc,

            unsigned locked_zone_ind,

            double &inclination,

            double &periapsis
        ) const;

        void angle_evolution_orbit_deriv(
            Dissipation::QuantityEntry entry,

            double angmom_deriv,

            DissipatingBody &body,

            unsigned zone_ind,

            double sin_inc,

            double cos_inc,

            unsigned locked_zone_ind,

            double &inclination,

            double &periapsis
        ) const;

        void fill_orbit_torque_deriv(
            Dissipation::QuantityEntry entry,

            DissipatingBody &body,

            unsigned zone_ind,

            std::valarray<Eigen::Vector3d> &orbit_torque_deriv
        ) const;

        void fill_zone_torque_deriv(
            Dissipation::QuantityEntry entry,

            DissipatingBody &body,

            unsigned zone_ind,

            std::valarray<Eigen::Vector3d> &zone_torque_deriv
        ) const;

        void inclination_evolution_zone_derivs(
            Dissipation::QuantityEntry entry,

            DissipatingBody &body,

            unsigned zone_ind,

            double zone_x_torque_above,

            double zone_x_torque_below,

            const std::valarray<Eigen::Vector3d> &zone_torque_deriv,

            const Eigen::Vector3d &orbit_torque,

            const std::valarray<Eigen::Vector3d> &orbit_torque_deriv,

            const std::valarray<Eigen::VectorXd> &above_frac_deriv,

            double sin_inc,

            double cos_inc,

            unsigned locked_zone_ind,

            double *result
        ) const;

        void periapsis_evolution_zone_derivs(
            Dissipation::QuantityEntry entry,

            DissipatingBody &body,

            unsigned zone_ind,

            double zone_y_torque_above,

            double zone_y_torque_below,

            const std::valarray<Eigen::Vector3d> &zone_torque_deriv,

            double orbit_y_torque,

            const std::valarray<Eigen::Vector3d> &orbit_torque_deriv,

            const std::valarray<Eigen::VectorXd> &above_frac_deriv,

            double sin_inc,

            double cos_inc,

            unsigned locked_zone_ind,

            double *result
        ) const;

        void spin_angmom_evolution_zone_derivs(
            Dissipation::QuantityEntry entry,

            DissipatingBody &body,

            unsigned zone_ind,

            double zone_z_torque_above,

            double zone_z_torque_below,

            const std::valarray<Eigen::Vector3d> &zone_torque_deriv,

            const std::valarray<Eigen::VectorXd> &above_frac_deriv,

            unsigned locked_zone_ind,

            double *result
        ) const;

        void binary_jacobian(
            double *param_derivs,

            double *age_derivs
        ) const;

        void fill_locked_surface_orbit(std::valarray<double> &orbit) const;

        void fill_binary_orbit(std::valarray<double> &orbit) const;

        void fill_single_orbit(std::valarray<double> &orbit) const;

    public:
        BinarySystem(
            DissipatingBody &body1,

            DissipatingBody &body2,

            const std::string &system_name=""
        ) :
        __name(system_name),
        __above_lock_fractions(Dissipation::NUM_DERIVATIVES),
        __body1(body1),
        __body2(body2)
        {}

        const std::string get_name() const {return __name;}

        virtual int configure(
            bool initialize,

            double age,

            double semimajor,

            double eccentricity,

            const double *spin_angmom,

            const double *inclination,

            const double *periapsis,

            Core::EvolModeType evolution_mode
        );

        int configure(
            bool initialize,

            double age,

            const double *parameters,

            Core::EvolModeType evolution_mode
        );

        double age() const {return __age;}

        const DissipatingBody &primary() const {return __body1;}

        const DissipatingBody &secondary() const {return __body2;}

        unsigned number_zones() const
        {return (__evolution_mode == Core::BINARY
                 ? __body1.number_zones() + __body2.number_zones()
                 : __body1.number_zones());}

        unsigned number_locked_zones() const
        {return __body1.number_locked_zones() + __body2.number_locked_zones();}

        double semimajor() const {return __semimajor;}

        double eccentricity() const {return __eccentricity;}

        Core::EvolModeType fill_orbit(
            std::valarray<double> &orbit
        ) const;


        double above_lock_fraction(
            unsigned locked_zone_index,

            Dissipation::QuantityEntry entry=Dissipation::NO_DERIV,

            unsigned deriv_zone_index=0,

            bool secondary_radius=false
        );

        int differential_equations(
            double age,

            const double *parameters,

            Core::EvolModeType evolution_mode,

            double *differential_equations,

            bool expansion_error=false
        );

        int jacobian(
            double age,

            const double *parameters,

            Core::EvolModeType evolution_mode,

            double *param_derivs,

            double *age_derivs
        );

        void check_for_lock(
            int orbital_freq_mult,

            int spin_freq_mult,

            unsigned short body_index,

            unsigned zone_index,

            short direction
        );

        virtual double minimum_separation(
            bool deriv=false
        ) const;

        Core::EvolModeType evolution_mode() {return __evolution_mode;}

        virtual void secondary_died();

        virtual void release_lock(
            unsigned locked_zone_index,

            short direction
        );

        virtual void add_to_evolution();

        virtual void reset_evolution();

        virtual void rewind_evolution(
            unsigned nsteps
        );

        virtual CombinedStoppingCondition *stopping_conditions();

        virtual void reached_critical_age(double age);

        virtual double next_stop_age() const;

        const std::list<double> &semimajor_evolution() const
        {return __semimajor_evolution;}

        const std::list<double> &semimajor_evolution_rate() const
        {return __semimajor_rate_evolution;}

        const std::list<double> &eccentricity_evolution() const
        {return __eccentricity_evolution;}

        const std::list<double> &eccentricity_evolution_rate() const
        {return __eccentricity_rate_evolution;}

        virtual void change_e_order(
            unsigned new_e_order
        )
        {
            __body1.change_e_order(new_e_order, *this, true);
            __body2.change_e_order(new_e_order, *this, false);
        }

        virtual unsigned eccentricity_order() const;

        virtual ~BinarySystem() {}

    }; //End BinarySystem class.

} //End Evolve namespace.

#endif