Squashed 'third_party/boostorg/odeint/' content from commit 6ff2719
Change-Id: If4892e29c1a5e6cf3a7aa51486a2725c251b0c7d
git-subtree-dir: third_party/boostorg/odeint
git-subtree-split: 6ff2719b6907b86596c3d43e88c1bcfdf29df560
diff --git a/examples/solar_system.cpp b/examples/solar_system.cpp
new file mode 100644
index 0000000..3586d5b
--- /dev/null
+++ b/examples/solar_system.cpp
@@ -0,0 +1,201 @@
+/* Boost libs/numeric/odeint/examples/solar_system.cpp
+
+ Copyright 2010-2012 Karsten Ahnert
+ Copyright 2011 Mario Mulansky
+
+ Solar system example for Hamiltonian stepper
+
+ 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)
+ */
+
+
+#include <iostream>
+#include <boost/array.hpp>
+
+#include <boost/numeric/odeint.hpp>
+
+#include "point_type.hpp"
+
+//[ container_type_definition
+// we simulate 5 planets and the sun
+const size_t n = 6;
+
+typedef point< double , 3 > point_type;
+typedef boost::array< point_type , n > container_type;
+typedef boost::array< double , n > mass_type;
+//]
+
+
+
+
+
+
+
+//[ coordinate_function
+const double gravitational_constant = 2.95912208286e-4;
+
+struct solar_system_coor
+{
+ const mass_type &m_masses;
+
+ solar_system_coor( const mass_type &masses ) : m_masses( masses ) { }
+
+ void operator()( const container_type &p , container_type &dqdt ) const
+ {
+ for( size_t i=0 ; i<n ; ++i )
+ dqdt[i] = p[i] / m_masses[i];
+ }
+};
+//]
+
+
+//[ momentum_function
+struct solar_system_momentum
+{
+ const mass_type &m_masses;
+
+ solar_system_momentum( const mass_type &masses ) : m_masses( masses ) { }
+
+ void operator()( const container_type &q , container_type &dpdt ) const
+ {
+ const size_t n = q.size();
+ for( size_t i=0 ; i<n ; ++i )
+ {
+ dpdt[i] = 0.0;
+ for( size_t j=0 ; j<i ; ++j )
+ {
+ point_type diff = q[j] - q[i];
+ double d = abs( diff );
+ diff *= ( gravitational_constant * m_masses[i] * m_masses[j] / d / d / d );
+ dpdt[i] += diff;
+ dpdt[j] -= diff;
+
+ }
+ }
+ }
+};
+//]
+
+
+
+
+
+
+
+//[ some_helpers
+point_type center_of_mass( const container_type &x , const mass_type &m )
+{
+ double overall_mass = 0.0;
+ point_type mean( 0.0 );
+ for( size_t i=0 ; i<x.size() ; ++i )
+ {
+ overall_mass += m[i];
+ mean += m[i] * x[i];
+ }
+ if( !x.empty() ) mean /= overall_mass;
+ return mean;
+}
+
+
+double energy( const container_type &q , const container_type &p , const mass_type &masses )
+{
+ const size_t n = q.size();
+ double en = 0.0;
+ for( size_t i=0 ; i<n ; ++i )
+ {
+ en += 0.5 * norm( p[i] ) / masses[i];
+ for( size_t j=0 ; j<i ; ++j )
+ {
+ double diff = abs( q[i] - q[j] );
+ en -= gravitational_constant * masses[j] * masses[i] / diff;
+ }
+ }
+ return en;
+}
+//]
+
+
+//[ streaming_observer
+struct streaming_observer
+{
+ std::ostream& m_out;
+
+ streaming_observer( std::ostream &out ) : m_out( out ) { }
+
+ template< class State >
+ void operator()( const State &x , double t ) const
+ {
+ container_type &q = x.first;
+ m_out << t;
+ for( size_t i=0 ; i<q.size() ; ++i ) m_out << "\t" << q[i];
+ m_out << "\n";
+ }
+};
+//]
+
+
+int main( int argc , char **argv )
+{
+
+ using namespace std;
+ using namespace boost::numeric::odeint;
+
+ mass_type masses = {{
+ 1.00000597682 , // sun
+ 0.000954786104043 , // jupiter
+ 0.000285583733151 , // saturn
+ 0.0000437273164546 , // uranus
+ 0.0000517759138449 , // neptune
+ 1.0 / ( 1.3e8 ) // pluto
+ }};
+
+ container_type q = {{
+ point_type( 0.0 , 0.0 , 0.0 ) , // sun
+ point_type( -3.5023653 , -3.8169847 , -1.5507963 ) , // jupiter
+ point_type( 9.0755314 , -3.0458353 , -1.6483708 ) , // saturn
+ point_type( 8.3101420 , -16.2901086 , -7.2521278 ) , // uranus
+ point_type( 11.4707666 , -25.7294829 , -10.8169456 ) , // neptune
+ point_type( -15.5387357 , -25.2225594 , -3.1902382 ) // pluto
+ }};
+
+ container_type p = {{
+ point_type( 0.0 , 0.0 , 0.0 ) , // sun
+ point_type( 0.00565429 , -0.00412490 , -0.00190589 ) , // jupiter
+ point_type( 0.00168318 , 0.00483525 , 0.00192462 ) , // saturn
+ point_type( 0.00354178 , 0.00137102 , 0.00055029 ) , // uranus
+ point_type( 0.00288930 , 0.00114527 , 0.00039677 ) , // neptune
+ point_type( 0.00276725 , -0.00170702 , -0.00136504 ) // pluto
+ }};
+
+ point_type qmean = center_of_mass( q , masses );
+ point_type pmean = center_of_mass( p , masses );
+ for( size_t i=0 ; i<n ; ++i )
+ {
+ q[i] -= qmean ;
+ p[i] -= pmean;
+ }
+
+ for( size_t i=0 ; i<n ; ++i ) p[i] *= masses[i];
+
+ //[ integration_solar_system
+ typedef symplectic_rkn_sb3a_mclachlan< container_type > stepper_type;
+ const double dt = 100.0;
+
+ integrate_const(
+ stepper_type() ,
+ make_pair( solar_system_coor( masses ) , solar_system_momentum( masses ) ) ,
+ make_pair( boost::ref( q ) , boost::ref( p ) ) ,
+ 0.0 , 200000.0 , dt , streaming_observer( cout ) );
+ //]
+
+
+ return 0;
+}
+
+
+/*
+Plot with gnuplot:
+p "solar_system.dat" u 2:4 w l,"solar_system.dat" u 5:7 w l,"solar_system.dat" u 8:10 w l,"solar_system.dat" u 11:13 w l,"solar_system.dat" u 14:16 w l,"solar_system.dat" u 17:19 w l
+ */