ceph/src/boost/libs/numeric/odeint/examples/solar_system.cpp

   1 /* Boost libs/numeric/odeint/examples/solar_system.cpp
   2
   3  Copyright 2010-2012 Karsten Ahnert
   4  Copyright 2011 Mario Mulansky
   5
   6  Solar system example for Hamiltonian stepper
   7
   8  Distributed under the Boost Software License, Version 1.0.
   9 (See accompanying file LICENSE_1_0.txt or
  10  copy at http://www.boost.org/LICENSE_1_0.txt)
  11  */
  12
  13
  14 #include <iostream>
  15 #include <boost/array.hpp>
  16
  17 #include <boost/numeric/odeint.hpp>
  18
  19 #include "point_type.hpp"
  20
  21 //[ container_type_definition
  22 // we simulate 5 planets and the sun
  23 const size_t n = 6;
  24
  25 typedef point< double , 3 > point_type;
  26 typedef boost::array< point_type , n > container_type;
  27 typedef boost::array< double , n > mass_type;
  28 //]
  29
  30
  31
  32
  33
  34
  35
  36 //[ coordinate_function
  37 const double gravitational_constant = 2.95912208286e-4;
  38
  39 struct solar_system_coor
  40 {
  41     const mass_type &m_masses;
  42
  43     solar_system_coor( const mass_type &masses ) : m_masses( masses ) { }
  44
  45     void operator()( const container_type &p , container_type &dqdt ) const
  46     {
  47         for( size_t i=0 ; i<n ; ++i )
  48             dqdt[i] = p[i] / m_masses[i];
  49     }
  50 };
  51 //]
  52
  53
  54 //[ momentum_function
  55 struct solar_system_momentum
  56 {
  57     const mass_type &m_masses;
  58
  59     solar_system_momentum( const mass_type &masses ) : m_masses( masses ) { }
  60
  61     void operator()( const container_type &q , container_type &dpdt ) const
  62     {
  63         const size_t n = q.size();
  64         for( size_t i=0 ; i<n ; ++i )
  65         {
  66             dpdt[i] = 0.0;
  67             for( size_t j=0 ; j<i ; ++j )
  68             {
  69                 point_type diff = q[j] - q[i];
  70                 double d = abs( diff );
  71                 diff *= ( gravitational_constant * m_masses[i] * m_masses[j] / d / d / d );
  72                 dpdt[i] += diff;
  73                 dpdt[j] -= diff;
  74
  75             }
  76         }
  77     }
  78 };
  79 //]
  80
  81
  82
  83
  84
  85
  86
  87 //[ some_helpers
  88 point_type center_of_mass( const container_type &x , const mass_type &m )
  89 {
  90     double overall_mass = 0.0;
  91     point_type mean( 0.0 );
  92     for( size_t i=0 ; i<x.size() ; ++i )
  93     {
  94         overall_mass += m[i];
  95         mean += m[i] * x[i];
  96     }
  97     if( !x.empty() ) mean /= overall_mass;
  98     return mean;
  99 }
 100
 101
 102 double energy( const container_type &q , const container_type &p , const mass_type &masses )
 103 {
 104     const size_t n = q.size();
 105     double en = 0.0;
 106     for( size_t i=0 ; i<n ; ++i )
 107     {
 108         en += 0.5 * norm( p[i] ) / masses[i];
 109         for( size_t j=0 ; j<i ; ++j )
 110         {
 111             double diff = abs( q[i] - q[j] );
 112             en -= gravitational_constant * masses[j] * masses[i] / diff;
 113         }
 114     }
 115     return en;
 116 }
 117 //]
 118
 119
 120 //[ streaming_observer
 121 struct streaming_observer
 122 {
 123     std::ostream& m_out;
 124
 125     streaming_observer( std::ostream &out ) : m_out( out ) { }
 126
 127     template< class State >
 128     void operator()( const State &x , double t ) const
 129     {
 130         container_type &q = x.first;
 131         m_out << t;
 132         for( size_t i=0 ; i<q.size() ; ++i ) m_out << "\t" << q[i];
 133         m_out << "\n";
 134     }
 135 };
 136 //]
 137
 138
 139 int main( int argc , char **argv )
 140 {
 141
 142     using namespace std;
 143     using namespace boost::numeric::odeint;
 144
 145     mass_type masses = {{
 146             1.00000597682 ,      // sun
 147             0.000954786104043 ,  // jupiter
 148             0.000285583733151 ,  // saturn
 149             0.0000437273164546 , // uranus
 150             0.0000517759138449 , // neptune
 151             1.0 / ( 1.3e8 )      // pluto
 152     }};
 153
 154     container_type q = {{
 155             point_type( 0.0 , 0.0 , 0.0 ) ,                        // sun
 156             point_type( -3.5023653 , -3.8169847 , -1.5507963 ) ,   // jupiter
 157             point_type( 9.0755314 , -3.0458353 , -1.6483708 ) ,    // saturn
 158             point_type( 8.3101420 , -16.2901086 , -7.2521278 ) ,   // uranus
 159             point_type( 11.4707666 , -25.7294829 , -10.8169456 ) , // neptune
 160             point_type( -15.5387357 , -25.2225594 , -3.1902382 )   // pluto
 161     }};
 162
 163     container_type p = {{
 164             point_type( 0.0 , 0.0 , 0.0 ) ,                        // sun
 165             point_type( 0.00565429 , -0.00412490 , -0.00190589 ) , // jupiter
 166             point_type( 0.00168318 , 0.00483525 , 0.00192462 ) ,   // saturn
 167             point_type( 0.00354178 , 0.00137102 , 0.00055029 ) ,   // uranus
 168             point_type( 0.00288930 , 0.00114527 , 0.00039677 ) ,   // neptune
 169             point_type( 0.00276725 , -0.00170702 , -0.00136504 )   // pluto
 170     }};
 171
 172     point_type qmean = center_of_mass( q , masses );
 173     point_type pmean = center_of_mass( p , masses );
 174     for( size_t i=0 ; i<n ; ++i )
 175     {
 176         q[i] -= qmean ;
 177         p[i] -= pmean;
 178     }
 179
 180     for( size_t i=0 ; i<n ; ++i ) p[i] *= masses[i];
 181
 182     //[ integration_solar_system
 183     typedef symplectic_rkn_sb3a_mclachlan< container_type > stepper_type;
 184     const double dt = 100.0;
 185
 186     integrate_const(
 187             stepper_type() ,
 188             make_pair( solar_system_coor( masses ) , solar_system_momentum( masses ) ) ,
 189             make_pair( boost::ref( q ) , boost::ref( p ) ) ,
 190             0.0 , 200000.0 , dt , streaming_observer( cout ) );
 191     //]
 192
 193
 194     return 0;
 195 }
 196
 197
 198 /*
 199 Plot with gnuplot:
 200 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
 201  */