1 /*
2  [auto_generated]
3  boost/numeric/odeint/integrate/detail/integrate_const.hpp
4 
5  [begin_description]
6  integrate const implementation
7  [end_description]
8 
9  Copyright 2012-2015 Mario Mulansky
10  Copyright 2012 Christoph Koke
11  Copyright 2012 Karsten Ahnert
12 
13  Distributed under the Boost Software License, Version 1.0.
14  (See accompanying file LICENSE_1_0.txt or
15  copy at http://www.boost.org/LICENSE_1_0.txt)
16  */
17 
18 #ifndef BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_CONST_HPP_INCLUDED
19 #define BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_CONST_HPP_INCLUDED
20 
21 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
22 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
23 #include <boost/numeric/odeint/util/unit_helper.hpp>
24 #include <boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp>
25 
26 #include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
27 
28 namespace boost {
29 namespace numeric {
30 namespace odeint {
31 namespace detail {
32 
33 // forward declaration
34 template< class Stepper , class System , class State , class Time , class Observer >
35 size_t integrate_adaptive(
36         Stepper stepper , System system , State &start_state ,
37         Time &start_time , Time end_time , Time &dt ,
38         Observer observer , controlled_stepper_tag
39 );
40 
41 
42 template< class Stepper , class System , class State , class Time , class Observer >
integrate_const(Stepper stepper,System system,State & start_state,Time start_time,Time end_time,Time dt,Observer observer,stepper_tag)43 size_t integrate_const(
44         Stepper stepper , System system , State &start_state ,
45         Time start_time , Time end_time , Time dt ,
46         Observer observer , stepper_tag
47 )
48 {
49 
50     typename odeint::unwrap_reference< Observer >::type &obs = observer;
51     typename odeint::unwrap_reference< Stepper >::type &st = stepper;
52 
53     Time time = start_time;
54     int step = 0;
55     // cast time+dt explicitely in case of expression templates (e.g. multiprecision)
56     while( less_eq_with_sign( static_cast<Time>(time+dt) , end_time , dt ) )
57     {
58         obs( start_state , time );
59         st.do_step( system , start_state , time , dt );
60         // direct computation of the time avoids error propagation happening when using time += dt
61         // we need clumsy type analysis to get boost units working here
62         ++step;
63         time = start_time + static_cast< typename unit_value_type<Time>::type >(step) * dt;
64     }
65     obs( start_state , time );
66 
67     return step;
68 }
69 
70 
71 
72 template< class Stepper , class System , class State , class Time , class Observer >
integrate_const(Stepper stepper,System system,State & start_state,Time start_time,Time end_time,Time dt,Observer observer,controlled_stepper_tag)73 size_t integrate_const(
74         Stepper stepper , System system , State &start_state ,
75         Time start_time , Time end_time , Time dt ,
76         Observer observer , controlled_stepper_tag
77 )
78 {
79     typename odeint::unwrap_reference< Observer >::type &obs = observer;
80 
81     Time time = start_time;
82     const Time time_step = dt;
83     int real_steps = 0;
84     int step = 0;
85 
86     while( less_eq_with_sign( static_cast<Time>(time+time_step) , end_time , dt ) )
87     {
88         obs( start_state , time );
89         // integrate_adaptive_checked uses the given checker to throw if an overflow occurs
90         real_steps += detail::integrate_adaptive(stepper, system, start_state, time,
91                                                  static_cast<Time>(time + time_step), dt,
92                                                  null_observer(), controlled_stepper_tag());
93         // direct computation of the time avoids error propagation happening when using time += dt
94         // we need clumsy type analysis to get boost units working here
95         step++;
96         time = start_time + static_cast< typename unit_value_type<Time>::type >(step) * time_step;
97     }
98     obs( start_state , time );
99 
100     return real_steps;
101 }
102 
103 
104 template< class Stepper , class System , class State , class Time , class Observer >
integrate_const(Stepper stepper,System system,State & start_state,Time start_time,Time end_time,Time dt,Observer observer,dense_output_stepper_tag)105 size_t integrate_const(
106         Stepper stepper , System system , State &start_state ,
107         Time start_time , Time end_time , Time dt ,
108         Observer observer , dense_output_stepper_tag
109 )
110 {
111     typename odeint::unwrap_reference< Observer >::type &obs = observer;
112     typename odeint::unwrap_reference< Stepper >::type &st = stepper;
113 
114     Time time = start_time;
115 
116     st.initialize( start_state , time , dt );
117     obs( start_state , time );
118     time += dt;
119 
120     int obs_step( 1 );
121     int real_step( 0 );
122 
123     while( less_eq_with_sign( static_cast<Time>(time+dt) , end_time , dt ) )
124     {
125         while( less_eq_with_sign( time , st.current_time() , dt ) )
126         {
127             st.calc_state( time , start_state );
128             obs( start_state , time );
129             ++obs_step;
130             // direct computation of the time avoids error propagation happening when using time += dt
131             // we need clumsy type analysis to get boost units working here
132             time = start_time + static_cast< typename unit_value_type<Time>::type >(obs_step) * dt;
133         }
134         // we have not reached the end, do another real step
135         if( less_with_sign( static_cast<Time>(st.current_time()+st.current_time_step()) ,
136                             end_time ,
137                             st.current_time_step() ) )
138         {
139             while( less_eq_with_sign( st.current_time() , time , dt ) )
140             {
141                 st.do_step( system );
142                 ++real_step;
143             }
144         }
145         else if( less_with_sign( st.current_time() , end_time , st.current_time_step() ) )
146         { // do the last step ending exactly on the end point
147             st.initialize( st.current_state() , st.current_time() , end_time - st.current_time() );
148             st.do_step( system );
149             ++real_step;
150         }
151 
152     }
153     // last observation, if we are still in observation interval
154     // might happen due to finite precision problems when computing the the time
155     if( less_eq_with_sign( time , end_time , dt ) )
156     {
157         st.calc_state( time , start_state );
158         obs( start_state , time );
159     }
160 
161     return real_step;
162 }
163 
164 
165 } } } }
166 
167 #endif
168