Chris@16: /*
Chris@16:  [auto_generated]
Chris@16:  boost/numeric/odeint/integrate/detail/integrate_times.hpp
Chris@16: 
Chris@16:  [begin_description]
Chris@16:  Default integrate times implementation.
Chris@16:  [end_description]
Chris@16: 
Chris@16:  Copyright 2009-2012 Karsten Ahnert
Chris@16:  Copyright 2009-2012 Mario Mulansky
Chris@16: 
Chris@16:  Distributed under the Boost Software License, Version 1.0.
Chris@16:  (See accompanying file LICENSE_1_0.txt or
Chris@16:  copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16:  */
Chris@16: 
Chris@16: 
Chris@16: #ifndef BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_TIMES_HPP_INCLUDED
Chris@16: #define BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_TIMES_HPP_INCLUDED
Chris@16: 
Chris@16: #include <stdexcept>
Chris@16: 
Chris@16: #include <boost/config.hpp>
Chris@16: #include <boost/numeric/odeint/util/unwrap_reference.hpp>
Chris@16: #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
Chris@16: #include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
Chris@16: 
Chris@16: 
Chris@16: namespace boost {
Chris@16: namespace numeric {
Chris@16: namespace odeint {
Chris@16: namespace detail {
Chris@16: 
Chris@16: 
Chris@16: 
Chris@16: /*
Chris@16:  * integrate_times for simple stepper
Chris@16:  */
Chris@16: template< class Stepper , class System , class State , class TimeIterator , class Time , class Observer >
Chris@16: size_t integrate_times(
Chris@16:         Stepper stepper , System system , State &start_state ,
Chris@16:         TimeIterator start_time , TimeIterator end_time , Time dt ,
Chris@16:         Observer observer , stepper_tag
Chris@16: )
Chris@16: {
Chris@16:     BOOST_USING_STD_MIN();
Chris@16: 
Chris@16:     typename odeint::unwrap_reference< Observer >::type &obs = observer;
Chris@16: 
Chris@16:     size_t steps = 0;
Chris@16:     Time current_dt = dt;
Chris@16:     while( true )
Chris@16:     {
Chris@16:         Time current_time = *start_time++;
Chris@16:         obs( start_state , current_time );
Chris@16:         if( start_time == end_time )
Chris@16:             break;
Chris@16:         while( less_with_sign( current_time , *start_time , current_dt ) )
Chris@16:         {
Chris@16:             current_dt = min BOOST_PREVENT_MACRO_SUBSTITUTION ( dt , *start_time - current_time );
Chris@16:             stepper.do_step( system , start_state , current_time , current_dt );
Chris@16:             current_time += current_dt;
Chris@16:             steps++;
Chris@16:         }
Chris@16:     }
Chris@16:     return steps;
Chris@16: }
Chris@16: 
Chris@16: /*
Chris@16:  * integrate_times for controlled stepper
Chris@16:  */
Chris@16: template< class Stepper , class System , class State , class TimeIterator , class Time , class Observer >
Chris@16: size_t integrate_times(
Chris@16:         Stepper stepper , System system , State &start_state ,
Chris@16:         TimeIterator start_time , TimeIterator end_time , Time dt ,
Chris@16:         Observer observer , controlled_stepper_tag
Chris@16: )
Chris@16: {
Chris@16:     BOOST_USING_STD_MIN();
Chris@16: 
Chris@16:     typename odeint::unwrap_reference< Observer >::type &obs = observer;
Chris@16: 
Chris@16:     const size_t max_attempts = 1000;
Chris@16:     const char *error_string = "Integrate adaptive : Maximal number of iterations reached. A step size could not be found.";
Chris@16:     size_t steps = 0;
Chris@16:     while( true )
Chris@16:     {
Chris@16:         size_t fail_steps = 0;
Chris@16:         Time current_time = *start_time++;
Chris@16:         obs( start_state , current_time );
Chris@16:         if( start_time == end_time )
Chris@16:             break;
Chris@16:         while( less_with_sign( current_time , *start_time , dt ) )
Chris@16:         {
Chris@16:             dt = min BOOST_PREVENT_MACRO_SUBSTITUTION ( dt , *start_time - current_time );
Chris@16:             if( stepper.try_step( system , start_state , current_time , dt ) == success )
Chris@16:             {
Chris@16:                 ++steps;
Chris@16:             }
Chris@16:             else
Chris@16:             {
Chris@16:                 ++fail_steps;
Chris@16:             }
Chris@16:             if( fail_steps == max_attempts ) throw std::overflow_error( error_string );
Chris@16:         }
Chris@16:     }
Chris@16:     return steps;
Chris@16: }
Chris@16: 
Chris@16: /*
Chris@16:  * integrate_times for dense output stepper
Chris@16:  */
Chris@16: template< class Stepper , class System , class State , class TimeIterator , class Time , class Observer >
Chris@16: size_t integrate_times(
Chris@16:         Stepper stepper , System system , State &start_state ,
Chris@16:         TimeIterator start_time , TimeIterator end_time , Time dt ,
Chris@16:         Observer observer , dense_output_stepper_tag
Chris@16: )
Chris@16: {
Chris@16:     typename odeint::unwrap_reference< Observer >::type &obs = observer;
Chris@16: 
Chris@16: 
Chris@16:     if( start_time == end_time )
Chris@16:         return 0;
Chris@16: 
Chris@16:     Time last_time_point = *(end_time-1);
Chris@16: 
Chris@16:     stepper.initialize( start_state , *start_time , dt );
Chris@16:     obs( start_state , *start_time++ );
Chris@16: 
Chris@16:     size_t count = 0;
Chris@16:     while( start_time != end_time )
Chris@16:     {
Chris@16:         while( ( start_time != end_time ) && less_eq_with_sign( *start_time , stepper.current_time() , stepper.current_time_step() ) )
Chris@16:         {
Chris@16:             stepper.calc_state( *start_time , start_state );
Chris@16:             obs( start_state , *start_time );
Chris@16:             start_time++;
Chris@16:         }
Chris@16: 
Chris@16:         // we have not reached the end, do another real step
Chris@16:         if( less_eq_with_sign( stepper.current_time() + stepper.current_time_step() ,
Chris@16:                                last_time_point ,
Chris@16:                                stepper.current_time_step() ) )
Chris@16:         {
Chris@16:             stepper.do_step( system );
Chris@16:             ++count;
Chris@16:         }
Chris@16:         else if( start_time != end_time )
Chris@16:         { // do the last step ending exactly on the end point
Chris@16:             stepper.initialize( stepper.current_state() , stepper.current_time() , last_time_point - stepper.current_time() );
Chris@16:             stepper.do_step( system );
Chris@16:             ++count;
Chris@16:         }
Chris@16:     }
Chris@16:     return count;
Chris@16: }
Chris@16: 
Chris@16: 
Chris@16: } // namespace detail
Chris@16: } // namespace odeint
Chris@16: } // namespace numeric
Chris@16: } // namespace boost
Chris@16: 
Chris@16: 
Chris@16: #endif // BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_ADAPTIVE_HPP_INCLUDED