File Integrator.hpp

File List > astrea > astro > astro > propagation > numerical > Integrator.hpp

Go to the documentation of this file

#pragma once

#include <vector>

#include <units/units.hpp>

#include <astro/astro.fwd.hpp>
#include <astro/propagation/event_detection/EventDetector.hpp>
#include <astro/state/State.hpp>
#include <astro/time/Interval.hpp>
#include <astro/types/typedefs.hpp>

namespace astrea {
namespace astro {

class Integrator {

  public:
    enum class StepMethod : EnumType {
        RK45,  
        RKF45, 
        RKF78, 
        DOP45, 
        DOP78, 
    };

    static inline Interval defaultInterval{ 0.0 * astrea::detail::time_unit, 86400.0 * astrea::detail::time_unit }; 

    Integrator() = default;

    ~Integrator() = default;

    StateHistory
        propagate(const State& state0, const Time& propTime, const EquationsOfMotion& eom, Vehicle vehicle, bool store = false, std::vector<Event> events = {});

    StateHistory
        propagate(const State& state0, const Date& endEpoch, const EquationsOfMotion& eom, Vehicle vehicle, bool store = false, std::vector<Event> events = {});

    void set_abs_tol(const Unitless& absTol);

    void set_rel_tol(const Unitless& relTol);

    void set_max_iter(const int& itMax);

    void switch_print(const bool& onOff);

    void switch_timer(const bool& onOff);

    void set_step_method(const StepMethod& stepMethod);

    void set_initial_timestep(const Time& dt0);

    void switch_fixed_timestep(const bool& onOff);

    void switch_fixed_timestep(const bool& onOff, const Time& fixedTimeStep);

    void set_timestep(const Time& fixedTimeStep);

    int n_func_evals() { return _functionEvaluations; }

  private:
    // Integrator constants
    const Unitless _EPSILON               = 0.8;    
    const Unitless _MIN_ERROR_TO_CATCH    = 2.0e-4; 
    const Unitless _MIN_ERROR_STEP_FACTOR = 5.0;    
    const Unitless _MIN_REL_STEP_SIZE     = 0.2;    

    // Iteration variables
    unsigned long _iteration          = 0;   
    unsigned _variableStepIteration   = 0;   
    unsigned long _MAX_ITER           = 1e8; 
    const unsigned _MAX_VAR_STEP_ITER = 1e3; 

    // Function evals
    int _functionEvaluations = 0; 

    // Time variables
    Time _timeStepPrevious; 
    Date _epoch0;           

    // Error variables
    Unitless _maxErrorPrevious; 

    // Butcher Tableau
    std::size_t _nStages{};                    
    static const std::size_t _MAX_STAGES = 13; 
    std::array<std::array<Unitless, _MAX_STAGES>, _MAX_STAGES> _a = {}; 
    std::array<Unitless, _MAX_STAGES> _b                          = {}; 
    std::array<Unitless, _MAX_STAGES> _bhat                       = {}; 
    std::array<Unitless, _MAX_STAGES> _db = {}; 
    std::array<Unitless, _MAX_STAGES> _c  = {}; 

    // ith order steps
    std::array<State, _MAX_STAGES> _kMatrix = {}; 
    State _statePlusKi;                           
    StatePartial _YFinalPrevious;                 

    // Clock variables
    clock_t _startClock{}; 
    clock_t _endClock{};   

    // Tolerances
    Unitless _ABS_TOL = 1.0e-13; 
    Unitless _REL_TOL = 1.0e-13; 

    // Initial step size
    Time _timeStepInitial = 60.0 * astrea::detail::time_unit; 

    // Run options
    bool _printOn = false; 
    bool _timerOn = false; 

    StepMethod _stepMethod = StepMethod::DOP45; 

    // Fake fixed step
    bool _useFixedStep  = false;                           
    Time _fixedTimeStep = 1.0 * astrea::detail::time_unit; 

    // Events
    EventDetector _eventDetector;

    StatePartial find_state_derivative(const Time& time, const State& state, const EquationsOfMotion& eom, Vehicle& vehicle);

    void setup(const std::vector<Event>& events);

    void teardown();

    void setup_butcher_tableau();

    bool try_step(Time& time, Time& timeStep, State& state, const EquationsOfMotion& eom, Vehicle& vehicle);

    Unitless find_max_error(const State& stateNew, const State& stateError) const;

    void take_fixed_step(Time& time, Time& timeStep, State& state, const EquationsOfMotion& eom, Vehicle& vehicle);

    std::pair<State, State>
        take_step(const Time& time, const Time& timeStep, const State& state, const EquationsOfMotion& eom, Vehicle& vehicle);

    bool check_error(const Unitless& maxError, const State& stateNew, const State& stateError, Time& time, Time& timeStep, State& state);

    Unitless get_relative_step_size(const Unitless& maxError) const;

    void store_final_func_eval(const Time& timeStep);

    void print_iteration(const Time& time, const State& state, const Time& timeFinal, const State& stateInitial);

    void print_performance() const;

    void startTimer();

    void endTimer();

    bool check_event(const Time& time, State& state, Vehicle& vehicle);

    bool validate_state_and_time(const Time& time, const State& state) const;
};

} // namespace astro
} // namespace astrea