libs/numeric/odeint/examples/harmonic_oscillator_units.cpp

   1 /*
   2  Copyright 2011-2013 Karsten Ahnert
   3  Copyright 2011-2013 Mario Mulansky
   4
   5  Distributed under the Boost Software License, Version 1.0.
   6  (See accompanying file LICENSE_1_0.txt or
   7  copy at http://www.boost.org/LICENSE_1_0.txt)
   8  */
   9
  10
  11 #include <iostream>
  12 #include <vector>
  13
  14
  15 /* WARNING: Compilation in debug mode might consume enormous memory
  16    (e.g. ~2GB on gcc 4.4 )
  17 */
  18
  19 // first increase fusion macro variables (now done by odeint itself)
  20 //#define BOOST_FUSION_INVOKE_MAX_ARITY 15
  21 //#define BOOST_RESULT_OF_NUM_ARGS 15
  22
  23 //[ units_define_basic_quantities
  24
  25 #include <boost/numeric/odeint.hpp>
  26 #include <boost/numeric/odeint/algebra/fusion_algebra.hpp>
  27 #include <boost/numeric/odeint/algebra/fusion_algebra_dispatcher.hpp>
  28
  29 #include <boost/units/systems/si/length.hpp>
  30 #include <boost/units/systems/si/time.hpp>
  31 #include <boost/units/systems/si/velocity.hpp>
  32 #include <boost/units/systems/si/acceleration.hpp>
  33 #include <boost/units/systems/si/io.hpp>
  34
  35 #include <boost/fusion/container.hpp>
  36
  37 using namespace std;
  38 using namespace boost::numeric::odeint;
  39 namespace fusion = boost::fusion;
  40 namespace units = boost::units;
  41 namespace si = boost::units::si;
  42
  43 typedef units::quantity< si::time , double > time_type;
  44 typedef units::quantity< si::length , double > length_type;
  45 typedef units::quantity< si::velocity , double > velocity_type;
  46 typedef units::quantity< si::acceleration , double > acceleration_type;
  47 typedef units::quantity< si::frequency , double > frequency_type;
  48
  49 typedef fusion::vector< length_type , velocity_type > state_type;
  50 typedef fusion::vector< velocity_type , acceleration_type > deriv_type;
  51 //]
  52
  53
  54
  55 //[ units_define_ode
  56 struct oscillator
  57 {
  58     frequency_type m_omega;
  59
  60     oscillator( const frequency_type &omega = 1.0 * si::hertz ) : m_omega( omega ) { }
  61
  62     void operator()( const state_type &x , deriv_type &dxdt , time_type t ) const
  63     {
  64         fusion::at_c< 0 >( dxdt ) = fusion::at_c< 1 >( x );
  65         fusion::at_c< 1 >( dxdt ) = - m_omega * m_omega * fusion::at_c< 0 >( x );
  66     }
  67 };
  68 //]
  69
  70
  71 //[ units_observer
  72 struct streaming_observer
  73 {
  74     std::ostream& m_out;
  75
  76     streaming_observer( std::ostream &out ) : m_out( out ) { }
  77
  78     struct write_element
  79     {
  80         std::ostream &m_out;
  81         write_element( std::ostream &out ) : m_out( out ) { };
  82
  83         template< class T >
  84         void operator()( const T &t ) const
  85         {
  86             m_out << "\t" << t;
  87         }
  88     };
  89
  90     template< class State , class Time >
  91     void operator()( const State &x , const Time &t ) const
  92     {
  93         m_out << t;
  94         fusion::for_each( x , write_element( m_out ) );
  95         m_out << "\n";
  96     }
  97 };
  98 //]
  99
 100
 101 int main( int argc , char**argv )
 102 {
 103 //    typedef dense_output_runge_kutta
 104 //    <
 105 //        controlled_runge_kutta
 106 //        <
 107 //            runge_kutta_dopri5< state_type , double , deriv_type , time_type , fusion_algebra >
 108 //        >
 109 //    > stepper_type;
 110
 111     //[ units_define_stepper
 112     typedef runge_kutta_dopri5< state_type , double , deriv_type , time_type > stepper_type;
 113
 114     state_type x( 1.0 * si::meter , 0.0 * si::meter_per_second );
 115
 116     integrate_const( make_dense_output( 1.0e-6 , 1.0e-6 , stepper_type() ) , oscillator( 2.0 * si::hertz ) ,
 117                      x , 0.0 * si::second , 100.0 * si::second , 0.1 * si::second , streaming_observer( cout ) );
 118     //]
 119
 120     return 0;
 121 }