Chris@16: /*
Chris@16:  [auto_generated]
Chris@16:  boost/numeric/odeint/integrate/detail/integrate_const.hpp
Chris@16: 
Chris@16:  [begin_description]
Chris@16:  integrate const implementation
Chris@16:  [end_description]
Chris@16: 
Chris@101:  Copyright 2012 Mario Mulansky
Chris@101:  Copyright 2012 Christoph Koke
Chris@101:  Copyright 2012 Karsten Ahnert
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: #ifndef BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_CONST_HPP_INCLUDED
Chris@16: #define BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_CONST_HPP_INCLUDED
Chris@16: 
Chris@16: #include <boost/numeric/odeint/util/unwrap_reference.hpp>
Chris@16: #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
Chris@16: #include <boost/numeric/odeint/util/unit_helper.hpp>
Chris@16: #include <boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp>
Chris@16: 
Chris@16: #include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
Chris@16: 
Chris@16: namespace boost {
Chris@16: namespace numeric {
Chris@16: namespace odeint {
Chris@16: namespace detail {
Chris@16: 
Chris@16: // forward declaration
Chris@16: template< class Stepper , class System , class State , class Time , class Observer >
Chris@16: size_t integrate_adaptive(
Chris@16:         Stepper stepper , System system , State &start_state ,
Chris@16:         Time &start_time , Time end_time , Time &dt ,
Chris@16:         Observer observer , controlled_stepper_tag
Chris@16: );
Chris@16: 
Chris@16: 
Chris@16: template< class Stepper , class System , class State , class Time , class Observer >
Chris@16: size_t integrate_const(
Chris@16:         Stepper stepper , System system , State &start_state ,
Chris@16:         Time start_time , Time 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@16:     
Chris@16:     Time time = start_time;
Chris@16:     int step = 0;
Chris@101:     // cast time+dt explicitely in case of expression templates (e.g. multiprecision)
Chris@101:     while( less_eq_with_sign( static_cast<Time>(time+dt) , end_time , dt ) )
Chris@16:     {
Chris@16:         obs( start_state , time );
Chris@101:         st.do_step( system , start_state , time , dt );
Chris@16:         // direct computation of the time avoids error propagation happening when using time += dt
Chris@16:         // we need clumsy type analysis to get boost units working here
Chris@16:         ++step;
Chris@16:         time = start_time + static_cast< typename unit_value_type<Time>::type >(step) * dt;
Chris@16:     }
Chris@16:     obs( start_state , time );
Chris@16: 
Chris@16:     return step;
Chris@16: }
Chris@16: 
Chris@16: 
Chris@16: 
Chris@16: template< class Stepper , class System , class State , class Time , class Observer >
Chris@16: size_t integrate_const(
Chris@16:         Stepper stepper , System system , State &start_state ,
Chris@16:         Time start_time , Time 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@16:     
Chris@16:     Time time = start_time;
Chris@16:     const Time time_step = dt;
Chris@101:     int real_steps = 0;
Chris@16:     int step = 0;
Chris@16:     
Chris@101:     while( less_eq_with_sign( static_cast<Time>(time+time_step) , end_time , dt ) )
Chris@16:     {
Chris@16:         obs( start_state , time );
Chris@101:         real_steps += detail::integrate_adaptive( stepper , system , start_state , time , time+time_step , dt ,
Chris@16:                                     null_observer() , controlled_stepper_tag() );
Chris@16:         // direct computation of the time avoids error propagation happening when using time += dt
Chris@16:         // we need clumsy type analysis to get boost units working here
Chris@101:         step++;
Chris@16:         time = start_time + static_cast< typename unit_value_type<Time>::type >(step) * time_step;
Chris@16:     }
Chris@16:     obs( start_state , time );
Chris@101: 
Chris@101:     return real_steps;
Chris@16: }
Chris@16: 
Chris@16: 
Chris@16: template< class Stepper , class System , class State , class Time , class Observer >
Chris@16: size_t integrate_const(
Chris@16:         Stepper stepper , System system , State &start_state ,
Chris@16:         Time start_time , Time 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@16:     
Chris@16:     Time time = start_time;
Chris@16:     
Chris@101:     st.initialize( start_state , time , dt );
Chris@16:     obs( start_state , time );
Chris@16:     time += dt;
Chris@16: 
Chris@16:     int obs_step( 1 );
Chris@16:     int real_step( 0 );
Chris@16:     
Chris@101:     while( less_with_sign( static_cast<Time>(time+dt) , end_time , dt ) )
Chris@16:     {
Chris@101:         while( less_eq_with_sign( time , st.current_time() , dt ) )
Chris@16:         {
Chris@101:             st.calc_state( time , start_state );
Chris@16:             obs( start_state , time );
Chris@16:             ++obs_step;
Chris@16:             // direct computation of the time avoids error propagation happening when using time += dt
Chris@16:             // we need clumsy type analysis to get boost units working here
Chris@16:             time = start_time + static_cast< typename unit_value_type<Time>::type >(obs_step) * dt;
Chris@16:         }
Chris@16:         // we have not reached the end, do another real step
Chris@101:         if( less_with_sign( static_cast<Time>(st.current_time()+st.current_time_step()) ,
Chris@16:                             end_time ,
Chris@101:                             st.current_time_step() ) )
Chris@16:         {
Chris@101:             while( less_eq_with_sign( st.current_time() , time , dt ) )
Chris@16:             {
Chris@101:                 st.do_step( system );
Chris@16:                 ++real_step;
Chris@16:             }
Chris@16:         }
Chris@101:         else if( less_with_sign( st.current_time() , end_time , st.current_time_step() ) )
Chris@16:         { // do the last step ending exactly on the end point
Chris@101:             st.initialize( st.current_state() , st.current_time() , end_time - st.current_time() );
Chris@101:             st.do_step( system );
Chris@16:             ++real_step;
Chris@16:         }
Chris@16:         
Chris@16:     }
Chris@16:     // last observation, if we are still in observation interval
Chris@16:     if( less_eq_with_sign( time , end_time , dt ) )
Chris@16:     {
Chris@101:         st.calc_state( time , start_state );
Chris@16:         obs( start_state , time );
Chris@16:     }
Chris@16:     
Chris@16:     return real_step;
Chris@16: }
Chris@16: 
Chris@16: 
Chris@16: } } } }
Chris@16: 
Chris@16: #endif