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@101:  Copyright 2011-2012 Mario Mulansky
Chris@101:  Copyright 2012 Karsten Ahnert
Chris@101:  Copyright 2012 Christoph Koke
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@101: #include <boost/throw_exception.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:     typename odeint::unwrap_reference< Observer >::type &obs = observer;
Chris@101:     typename odeint::unwrap_reference< Stepper >::type &st = stepper;
Chris@101:     typedef typename unit_value_type<Time>::type time_type;
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@101:         while( less_with_sign( current_time , static_cast<time_type>(*start_time) , current_dt ) )
Chris@16:         {
Chris@101:             current_dt = min_abs( dt , *start_time - current_time );
Chris@101:             st.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:     typename odeint::unwrap_reference< Observer >::type &obs = observer;
Chris@101:     typename odeint::unwrap_reference< Stepper >::type &st = stepper;
Chris@101:     typedef typename unit_value_type<Time>::type time_type;
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@101:         while( less_with_sign( current_time , static_cast<time_type>(*start_time) , dt ) )
Chris@16:         {
Chris@101:             // adjust stepsize to end up exactly at the observation point
Chris@101:             Time current_dt = min_abs( dt , *start_time - current_time );
Chris@101:             if( st.try_step( system , start_state , current_time , current_dt ) == success )
Chris@16:             {
Chris@16:                 ++steps;
Chris@101:                 // continue with the original step size if dt was reduced due to observation
Chris@101:                 dt = max_abs( dt , current_dt );
Chris@16:             }
Chris@16:             else
Chris@16:             {
Chris@16:                 ++fail_steps;
Chris@101:                 dt = current_dt;
Chris@16:             }
Chris@101:             if( fail_steps == max_attempts ) BOOST_THROW_EXCEPTION( 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@101:     typename odeint::unwrap_reference< Stepper >::type &st = stepper;
Chris@101:     typedef typename unit_value_type<Time>::type time_type;
Chris@16: 
Chris@16:     if( start_time == end_time )
Chris@16:         return 0;
Chris@16: 
Chris@101:     Time last_time_point = static_cast<time_type>(*(end_time-1));
Chris@16: 
Chris@101:     st.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@101:         while( ( start_time != end_time ) && less_eq_with_sign( static_cast<time_type>(*start_time) , st.current_time() , st.current_time_step() ) )
Chris@16:         {
Chris@101:             st.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@101:         if( less_eq_with_sign( st.current_time() + st.current_time_step() ,
Chris@16:                                last_time_point ,
Chris@101:                                st.current_time_step() ) )
Chris@16:         {
Chris@101:             st.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@101:             st.initialize( st.current_state() , st.current_time() , last_time_point - st.current_time() );
Chris@101:             st.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