Subject: [Boost-commit] svn:boost r59988 - in sandbox/odeint: . boost/numeric boost/numeric/odeint boost/numeric/odeint/concepts boost/numeric/odeint/detail libs/numeric/odeint/doc libs/numeric/odeint/examples
From: karsten.ahnert_at_[hidden]
Date: 2010-02-28 05:22:45
Author: karsten
Date: 2010-02-28 05:22:44 EST (Sun, 28 Feb 2010)
New Revision: 59988
URL: http://svn.boost.org/trac/boost/changeset/59988
Log:
small changes and introducing the solar system example in the tutorial section
Added:
sandbox/odeint/boost/numeric/odeint/gram_schmitt.hpp (contents, props changed)
sandbox/odeint/boost/numeric/odeint/stepper_base.hpp (contents, props changed)
sandbox/odeint/libs/numeric/odeint/examples/solar_system.cpp (contents, props changed)
Text files modified:
sandbox/odeint/ToDo | 41 ++++++++++++++++++
sandbox/odeint/boost/numeric/odeint.hpp | 2
sandbox/odeint/boost/numeric/odeint/concepts/state_concept.hpp | 1
sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp | 13 ++++-
sandbox/odeint/boost/numeric/odeint/stepper_euler.hpp | 38 +++++++++++-----
sandbox/odeint/boost/numeric/odeint/stepper_rk78_fehlberg.hpp | 2
sandbox/odeint/libs/numeric/odeint/doc/Jamfile | 6 -
sandbox/odeint/libs/numeric/odeint/doc/tutorial.qbk | 87 +++++++++++++++++++++++++++++++++++++++
sandbox/odeint/libs/numeric/odeint/examples/Jamfile | 2
9 files changed, 168 insertions(+), 24 deletions(-)
Modified: sandbox/odeint/ToDo
==============================================================================
--- sandbox/odeint/ToDo (original)
+++ sandbox/odeint/ToDo 2010-02-28 05:22:44 EST (Sun, 28 Feb 2010)
@@ -3,10 +3,21 @@
Karsten:
* Resizer functionality should be called explicitly and not in every call of do_step or try_step
* Add new order function (because we need three orders for rk78)
-* deriv every stepper from stepper_base< stepper >, because do_step( system , x , t , dt ) is only implemented ones ( and calls the special do_step( system , x , dxdt , t , dt )
+* deriv every stepper from stepper_base< stepper >, because do_step( system , x , t , dt ) is only implemented once ( and calls the special do_step( system , x , dxdt , t , dt )
* In all steppers:
* const unsigned short -> short
* typedef typename traits_type::container_type container_type statt typedef Container container_type
+* in iterator algebra die funktionen welche maximas suchen durch die entsprechenden boost funktionen ersetzen
+* Concept checks, für vernünftige Fehlermeldungen
+ * auch für checks ob abs() vorhanden ist
+* controlled stepper standard noncopyable machen
+* controlled stepper bs noncopyable
+* all stepper noncopyable
+* check in all steppers the order
+* rk78 die fehler implementieren
+
+
+
Mario:
* code cleanup:
@@ -64,3 +75,31 @@
+
+
+
+
+
+
+= Wish list =
+
+* dense output
+* replacable iterator algebra, with specialized algebras for parallelized calculations / cuda ...
+* implicit methods
+* Dormand Prince methods 54 und 8
+* error methods for symplectic integrators
+
+
+
+
+
+
+
+
+
+
+
+
+= Some issues =
+
+* zip iterator is slow, hence writing the iterator algebra with zip iterators is not an option
Modified: sandbox/odeint/boost/numeric/odeint.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint.hpp 2010-02-28 05:22:44 EST (Sun, 28 Feb 2010)
@@ -41,4 +41,6 @@
#include <boost/numeric/odeint/integrator_adaptive_stepsize.hpp>
#include <boost/numeric/odeint/integrator_constant_stepsize.hpp>
+#include <boost/numeric/odeint/gram_schmitt.hpp>
+
#endif // BOOST_NUMERIC_ODEINT_HPP
Modified: sandbox/odeint/boost/numeric/odeint/concepts/state_concept.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/concepts/state_concept.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/concepts/state_concept.hpp 2010-02-28 05:22:44 EST (Sun, 28 Feb 2010)
@@ -26,6 +26,7 @@
public:
typedef typename X::iterator iterator; // requires iterator typedef
+ typedef typename X::value_type value_type; // requires value_type typedef
// requires iterator being ForwardIterator
BOOST_CONCEPT_ASSERT((ForwardIterator<iterator>));
Modified: sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp 2010-02-28 05:22:44 EST (Sun, 28 Feb 2010)
@@ -541,11 +541,13 @@
T a_x,
T a_dxdt )
{
+ using std::abs;
+
while( y_begin != y_end )
{
*y_begin++ = eps_abs +
- eps_rel * (a_x * std::abs(*x1_begin++) +
- a_dxdt * std::abs(*x2_begin++));
+ eps_rel * (a_x * abs(*x1_begin++) +
+ a_dxdt * abs(*x2_begin++));
}
}
@@ -559,10 +561,13 @@
InputIterator2 x2_begin,
T initial_max )
{
+ using std::abs;
+
while( x1_begin != x1_end )
{
- initial_max = std::max( static_cast<T>(std::abs(*x1_begin++)/(*x2_begin++)),
- initial_max);
+ initial_max = std::max(
+ static_cast<T>( abs(*x1_begin++)/abs(*x2_begin++)),
+ initial_max );
}
return initial_max;
}
Added: sandbox/odeint/boost/numeric/odeint/gram_schmitt.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/boost/numeric/odeint/gram_schmitt.hpp 2010-02-28 05:22:44 EST (Sun, 28 Feb 2010)
@@ -0,0 +1,84 @@
+/*
+ boost header: numeric/odeint/gram_schmitt.hpp
+
+ Copyright 2009 Karsten Ahnert
+ Copyright 2009 Mario Mulansky
+ Copyright 2009 Andre Bergner
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+*/
+
+#ifndef BOOST_NUMERIC_ODEINT_GRAM_SCHMITT_HPP_INCLUDED
+#define BOOST_NUMERIC_ODEINT_GRAM_SCHMITT_HPP_INCLUDED
+
+#include <iterator>
+#include <algorithm>
+#include <numeric>
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+ template< class Iterator , class T >
+ void normalize( Iterator first , Iterator last , T norm )
+ {
+ while( first != last ) *first++ /= norm;
+ }
+
+ template< class Iterator , class T >
+ void substract_vector( Iterator first1 , Iterator last1 ,
+ Iterator first2 , T val )
+ {
+ while( first1 != last1 )
+ *first1++ -= val * *first2++;
+ }
+
+ template< class StateType , class LyapType >
+ void gram_schmitt( StateType &x , LyapType &lyap ,
+ size_t n , size_t num_of_lyap )
+ {
+ if( !num_of_lyap ) return;
+ if( ptrdiff_t( ( num_of_lyap + 1 ) * n ) != std::distance( x.begin() , x.end() ) )
+ throw std::domain_error( "renormalization() : size of state does not match the number of lyapunov exponents." );
+
+ typedef typename StateType::value_type value_type;
+ typedef typename StateType::iterator iterator;
+
+ value_type norm[num_of_lyap] , tmp[num_of_lyap];
+ iterator first = x.begin() + n;
+ iterator beg1 = first , end1 = first + n ;
+
+ std::fill( norm , norm+num_of_lyap , 0.0 );
+
+ // normalize first vector
+ norm[0] = sqrt( std::inner_product( beg1 , end1 , beg1 , 0.0 ) );
+ normalize( beg1 , end1 , norm[0] );
+
+ beg1 += n ;
+ end1 += n;
+
+ for( size_t j=1 ; j<num_of_lyap ; ++j , beg1+=n , end1+=n )
+ {
+ for( size_t k=0 ; k<j ; ++k )
+ tmp[k] = std::inner_product( beg1 , end1 , first + k*n , 0.0 );
+
+ for( size_t k=0 ; k<j ; ++k )
+ substract_vector( beg1 , end1 , first + k*n , tmp[k] );
+
+ // nromalize j-th vector
+ norm[j] = sqrt( std::inner_product( beg1 , end1 , beg1 , 0.0 ) );
+ normalize( beg1 , end1 , norm[j] );
+ }
+
+ for( size_t j=0 ; j<num_of_lyap ; j++ )
+ lyap[j] += log( norm[j] );
+ }
+
+
+} // namespace odeint
+} // namespace numeric
+} // namespace boost
+
+#endif //BOOST_NUMERIC_ODEINT_GRAM_SCHMITT_HPP_INCLUDED
Added: sandbox/odeint/boost/numeric/odeint/stepper_base.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/boost/numeric/odeint/stepper_base.hpp 2010-02-28 05:22:44 EST (Sun, 28 Feb 2010)
@@ -0,0 +1,156 @@
+/*
+ boost header: numeric/odeint/stepper_base.hpp
+
+ Copyright 2009 Karsten Ahnert
+ Copyright 2009 Mario Mulansky
+ Copyright 2009 Andre Bergner
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+*/
+
+#ifndef BOOST_NUMERIC_ODEINT_STEPPER_BASE_HPP_INCLUDED
+#define BOOST_NUMERIC_ODEINT_STEPPER_BASE_HPP_INCLUDED
+
+#include <boost/noncopyable.hpp>
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+ typedef unsigned char order_type;
+
+ template<
+ class Stepper ,
+ class Container ,
+ order_type Order ,
+ class Time ,
+ class Traits
+ >
+ class explicit_stepper_base : private boost::noncopyable
+ {
+ public:
+
+ // some typedef
+
+ typedef Stepper stepper_type;
+ typedef Time time_type;
+ typedef Traits traits_type;
+ typedef typename traits_type::container_type container_type;
+ typedef typename traits_type::value_type value_type;
+ typedef typename traits_type::iterator iterator;
+ typedef typename traits_type::const_iterator const_iterator;
+
+
+
+
+ public:
+
+ // provide some functions
+
+ explicit_stepper_base( stepper_type &stepper )
+ : m_stepper( stepper ) { }
+
+ explicit_stepper_base( stepper_type &stepper , container_type &x )
+ : m_stepper( stepper )
+ {
+ traits_type::adjust_size( x , m_dxdt );
+ }
+
+// order_type order() const { return m_order; }
+ const static order_type order = Order;
+
+ template< class DynamicalSystem >
+ void do_step( DynamicalSystem &system ,
+ container_type &x ,
+ time_type t ,
+ time_type dt )
+ {
+ system( x , m_dxdt , t );
+ stepper.do_step_with_deriv( system , x , m_dxdt , t , dt );
+ }
+
+ void adjust_size( container_type &x )
+ {
+ traits_type::adjust_size( x , m_dxdt );
+ }
+
+ private:
+
+ container_type m_dxdt;
+ stepper_type &m_stepper;
+ };
+
+
+
+
+ template<
+ class ErrorStepper ,
+ class Container ,
+ order_type StepperOrder ,
+ order_type ErrorOrder
+ class Time ,
+ class Traits
+ >
+ class explicit_error_stepper_base : private boost::noncopyable
+ {
+ public:
+
+ // some typedef
+
+ typedef ErrorStepper stepper_type;
+ typedef Time time_type;
+ typedef Traits traits_type;
+ typedef typename traits_type::container_type container_type;
+ typedef typename traits_type::value_type value_type;
+ typedef typename traits_type::iterator iterator;
+ typedef typename traits_type::const_iterator const_iterator;
+
+
+
+
+ public:
+
+ // provide some functions
+
+ explicit_stepper_base( stepper_type &stepper )
+ : m_stepper( stepper ) { }
+
+ explicit_stepper_base( stepper_type &stepper , container_type &x )
+ : m_stepper( stepper )
+ {
+ traits_type::adjust_size( x , m_dxdt );
+ }
+
+// order_type order() const { return m_order; }
+ const static order_type stepper_order = StepperOrder;
+ const static order_type error_order = ErrorOrder;
+
+ template< class DynamicalSystem >
+ void do_step( DynamicalSystem &system ,
+ container_type &x ,
+ time_type t ,
+ time_type dt ,
+ container_type &err )
+ {
+ system( x , m_dxdt , t );
+ stepper.do_step_with_deriv( system , x , m_dxdt , t , dt , err );
+ }
+
+ void adjust_size( container_type &x )
+ {
+ traits_type::adjust_size( x , m_dxdt );
+ }
+
+ private:
+
+ container_type m_dxdt;
+ stepper_type &m_stepper;
+ };
+
+}
+}
+}
+
+#endif //BOOST_NUMERIC_ODEINT_STEPPER_BASE_HPP_INCLUDED
Modified: sandbox/odeint/boost/numeric/odeint/stepper_euler.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepper_euler.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/stepper_euler.hpp 2010-02-28 05:22:44 EST (Sun, 28 Feb 2010)
@@ -21,6 +21,8 @@
#include <boost/numeric/odeint/detail/iterator_algebra.hpp>
#include <boost/numeric/odeint/container_traits.hpp>
+#include <boost/numeric/odeint/stepper_base.hpp>
+
namespace boost {
namespace numeric {
@@ -31,13 +33,26 @@
class Time = double ,
class Traits = container_traits< Container >
>
- class stepper_euler
+ class stepper_euler : public stepper_base<
+ stepper_euler< Container , Time , Traits > ,
+ Container ,
+ 1 ,
+ Time ,
+ Traits >
{
// provide basic typedefs
public:
+ typedef stepper_base< stepper_euler< Container , Time , Traits > ,
+ Container , 1 , Time , Traits > base_type;
+
+ typedef typename base_type::time_type time_type;
+// typedef typename base_type::order_type order_type;
+/*
typedef unsigned short order_type;
typedef Time time_type;
+*/
+
typedef Traits traits_type;
typedef typename traits_type::container_type container_type;
typedef typename traits_type::value_type value_type;
@@ -56,14 +71,13 @@
// public interface
public:
- order_type order() const { return 1; }
-
+// order_type order() const { return 1; }
template< class DynamicalSystem >
void do_step( DynamicalSystem &system ,
- container_type &x ,
- const container_type &dxdt ,
+ container_type &x ,
+ const container_type &dxdt ,
time_type t ,
time_type dt )
{
@@ -74,16 +88,16 @@
dt );
}
- template< class DynamicalSystem >
- void do_step( DynamicalSystem &system ,
- container_type &x ,
- time_type t ,
- time_type dt )
- {
+/* template< class DynamicalSystem >
+ void do_step( DynamicalSystem &system ,
+ container_type &x ,
+ time_type t ,
+ time_type dt )
+ {
traits_type::adjust_size( x , m_dxdt );
system( x , m_dxdt , t );
do_step( system , x , m_dxdt , t , dt );
- }
+ }*/
};
Modified: sandbox/odeint/boost/numeric/odeint/stepper_rk78_fehlberg.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepper_rk78_fehlberg.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/stepper_rk78_fehlberg.hpp 2010-02-28 05:22:44 EST (Sun, 28 Feb 2010)
@@ -61,6 +61,8 @@
order_type order() const { return 7; }
+ order_type order_error() const { return 8; }
+
template< class DynamicalSystem >
void do_step( DynamicalSystem &system ,
container_type &x ,
Modified: sandbox/odeint/libs/numeric/odeint/doc/Jamfile
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/doc/Jamfile (original)
+++ sandbox/odeint/libs/numeric/odeint/doc/Jamfile 2010-02-28 05:22:44 EST (Sun, 28 Feb 2010)
@@ -4,13 +4,11 @@
odeint.qbk
;
-boostbook standalone
- :
+boostbook standalone :
odeint
;
-install html
- :
+install html :
/boost//doc/html/boostbook.css
;
Modified: sandbox/odeint/libs/numeric/odeint/doc/tutorial.qbk
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/doc/tutorial.qbk (original)
+++ sandbox/odeint/libs/numeric/odeint/doc/tutorial.qbk 2010-02-28 05:22:44 EST (Sun, 28 Feb 2010)
@@ -126,11 +126,94 @@
[endsect]
+
+
+
+
+
+
+
+
+
+
+
+
+
[section Solar system]
-define point types via operators
+[section Gravitation and energy conservation]
+
+The next example in this tutorial is a simulation of the solar system. In the
+solar system each planet, and of course also the sun will be represented by
+mass points. The interaction force between each object is the gravitational
+force which can be written as
+Fij = - gamma m M / (qi-qj)^3 * (qi-qj)
+
+where gamma is the gravitational constant, mi and mj are the masses and qi and
+qj are the locations of the two objects.
+
+dqi = pi
+dpi = 1/m sumj\newi Fij
+
+where pi is the momenta of object i. The equations of motion can also be
+derived from the Hamiltonian
+
+H = sum over i pi^2 / 2 + sum j V( qi , qj )
+
+via dqi = dH / dpi, dpi = - dH / dq1. V(qi,qj) is the interaction
+potential.
+
+In time independent Hamiltonian system the energy is conserved and special
+integration methods have to be applied in order to ensure energy
+conservation. The odeint library provides two stepper classes for Hamiltonian
+systems, which are separable and can be written in the form H = sum pi^2/2 +
+Hq.
+
+hamiltonian_stepper_euler
+hamiltonian_stepper_rk
+
+Alltough this functional form might look a bit arbitrary it covers nearly all
+classical mechanical systems with inertia and without dissipation, or where
+the equations of motion can be written in the form dqi=pi dpi=f(qi).
+
+
+
+
+
+[section Define the system function]
+
+To implement this system we define a point type which will represent the space
+as well as the velocity. Therefore, we use the operators from
+<boost/operator.hpp>:
+
+show the code
+
+
+The next step is to define the state type and the system (derivative)
+function. As state type we use std::tr1::array and a state type represents all
+space coordinates q or all momenta coordinates p. As system function we have
+to provide f(q)
+
+show the code
+
+Note, that we have allready define the masses of all planets in the solar
+system.
+
+In general a three body-system is chaotic, hence we can not expect that
+arbitray initial conditions of the system will lead to a dynamic which is
+comparable with the solar system. That is we have to define proper initial
+conditions.
+
+show the code
+
+Now, we use the rk stepper to integrate the solar system. To visualize the
+motion we save the trajectory of each planet in a circular buffer. The output
+can be piped directly into gnuplot and a very nice visualization of the motion
+appears.
+
+[endsect]
+
-define the system
usage of the steppers
Modified: sandbox/odeint/libs/numeric/odeint/examples/Jamfile
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/Jamfile (original)
+++ sandbox/odeint/libs/numeric/odeint/examples/Jamfile 2010-02-28 05:22:44 EST (Sun, 28 Feb 2010)
@@ -15,7 +15,7 @@
;
exe harmonic_oscillator : harmonic_oscillator.cpp ;
-exe test_container_and_stepper : test_container_and_stepper.cpp ;
+# exe test_container_and_stepper : test_container_and_stepper.cpp ;
exe hamiltonian_test : hamiltonian_test.cpp ;
exe lorenz_cmp_rk4_rk_generic : lorenz_cmp_rk4_rk_generic.cpp ;
Added: sandbox/odeint/libs/numeric/odeint/examples/solar_system.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/libs/numeric/odeint/examples/solar_system.cpp 2010-02-28 05:22:44 EST (Sun, 28 Feb 2010)
@@ -0,0 +1,215 @@
+#include <bits/stdc++.h>
+#include <bits/stdtr1c++.h>
+
+#include <boost/operators.hpp>
+#include <boost/circular_buffer.hpp>
+#include <boost/numeric/odeint.hpp>
+
+#define tab "\t"
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+template< class T >
+class point :
+ boost::additive1< point< T > ,
+ boost::additive2< point< T > , T ,
+ boost::multiplicative2< point< T > , T
+ > > >
+{
+public:
+
+ typedef T value_type;
+ typedef point< value_type > point_type;
+
+ value_type x , y , z;
+
+ point( void ) : x(0.0) , y(0.0) , z(0.0) { }
+
+ point( value_type val ) : x(val) , y(val) , z(val) { }
+
+ point( value_type _x , value_type _y , value_type _z )
+ : x(_x) , y(_y) , z(_z) { }
+
+
+ point_type& operator+=( const point_type& p )
+ {
+ x += p.x ; y += p.y ; z += p.z;
+ return *this;
+ }
+
+ point_type& operator-=( const point_type& p )
+ {
+ x -= p.x ; y -= p.y ; z -= p.z;
+ return *this;
+ }
+
+ point_type& operator+=( const value_type& val )
+ {
+ x += val ; y += val ; z += val;
+ return *this;
+ }
+
+ point_type& operator-=( const value_type& val )
+ {
+ x -= val ; y -= val ; z -= val;
+ return *this;
+ }
+
+ point_type& operator*=( const value_type &val )
+ {
+ x *= val ; y *= val ; z *= val;
+ return *this;
+ }
+
+ point_type& operator/=( const value_type &val )
+ {
+ x /= val ; y /= val ; z /= val;
+ return *this;
+ }
+};
+
+template< class T >
+T inner_product( const point< T > &p1 , const point< T > &p2 )
+{
+ return p1.x*p2.x + p1.y*p2.y + p1.z*p2.z;
+}
+
+template< class T >
+T norm( const point< T > &p )
+{
+ return inner_product( p , p );
+}
+
+template< class T >
+T abs( const point< T > &p )
+{
+ return sqrt( norm( p ) );
+}
+
+template< class T >
+ostream& operator<<( ostream &out , const point< T > &p )
+{
+ out << p.x << tab << p.y << tab << p.z;
+ return out;
+}
+
+
+
+const size_t n = 3;
+typedef point< double > point_type;
+typedef std::tr1::array< point_type , n > state_type;
+typedef std::tr1::array< double , n > mass_type;
+
+typedef hamiltonian_stepper_rk< state_type > stepper_type;
+
+typedef boost::circular_buffer< point_type > buffer_type;
+
+
+const mass_type masses = {{ 1.0e9 , 1.0e9 , 1.0e12}};
+const double gravitational_constant = 6.657e-20;
+
+ostream& operator<<( ostream &out , const state_type &x )
+{
+ typedef state_type::value_type value_type;
+ copy( x.begin() , x.end() ,
+ ostream_iterator< value_type >( out , "\n" ) );
+ return out;
+}
+
+point_type get_mean( const state_type &x )
+{
+ point_type mean( 0.0 );
+ if( x.empty() ) return mean;
+ for( size_t i=0 ; i<x.size() ; ++i ) mean += x[i];
+ mean /= double( x.size() );
+ return mean;
+}
+
+point_type get_center_of_mass( const state_type &x ,
+ const mass_type &m )
+{
+ point_type mean( 0.0 );
+ if( x.empty() ) return mean;
+ double overall_mass = 0.0;
+ for( size_t i=0 ; i<x.size() ; ++i )
+ {
+ overall_mass += m[i];
+ mean += m[i] * x[i];
+ }
+ mean /= overall_mass;
+ return mean;
+
+}
+
+void center_system( state_type &x , point_type mean )
+{
+ for( size_t i=0 ; i<x.size() ; ++i ) x[i] -= mean;
+}
+
+
+void solar_system( state_type &q , state_type &dpdt )
+{
+ point_type diff , tmp;
+ fill( dpdt.begin() , dpdt.end() , 0.0 );
+ for( size_t i=0 ; i<n ; ++i )
+ {
+ for( size_t j=i+1 ; j<n ; ++j )
+ {
+ diff = q[j] - q[i];
+ tmp = gravitational_constant * diff / pow( abs( diff ) , 3.0 );
+ dpdt[i] += tmp * masses[j];
+ dpdt[j] -= tmp * masses[i];
+ }
+ }
+}
+
+
+int main( int argc , char **argv )
+{
+ state_type q , p;
+ stepper_type stepper;
+
+ fill( q.begin() , q.end() , 0.0 );
+ fill( p.begin() , p.end() , 0.0 );
+ q[0] = point_type( 0.0 , 1.0 , 0.0 );
+ p[0] = point_type( 0.00001 , 0.0 , 0.0 );
+ q[2] = point_type( 1.0 , 0.0 , 0.0 );
+
+ center_system( q , get_center_of_mass( q , masses ) );
+ center_system( p , get_center_of_mass( p , masses ) );
+
+ const double dt = 1.0;
+
+ buffer_type buffers[n];
+ for( size_t i=0 ; i<n ; ++i ) buffers[i].set_capacity( 100 );
+
+ cout << "unset key\n";
+ const size_t ilen = 1000;
+ double t = 0.0;
+ while( true )
+ {
+ clog << get_mean( p ) << tab << get_mean( q ) << endl;
+ for( size_t i=0 ; i<n ; ++i ) buffers[i].push_back( q[i] );
+
+ cout << "p [-20:20][-20:20] '-' u 1:2 w l\n";
+ for( size_t i=0 ; i<n ; ++i )
+ {
+ copy( buffers[i].begin() , buffers[i].end() ,
+ ostream_iterator<point_type>( cout , "\n" ) );
+ cout << "\n";
+ }
+ cout << "e" << endl;
+
+// getchar();
+
+// cout << "p [-2:2][-2:2] '-' u 1:2 pt 7 ps 5 \n" << q << "e" << endl;
+
+ for( size_t ii=0 ; ii<ilen ; ++ii,t+=dt )
+ {
+ stepper.do_step( solar_system , q , p , dt );
+ }
+ }
+
+ return 0;
+}