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
solar_system_coorsolar_system_coor43 solar_system_coor( const mass_type &masses ) : m_masses( masses ) { }
44
operator ()solar_system_coor45 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
solar_system_momentumsolar_system_momentum59 solar_system_momentum( const mass_type &masses ) : m_masses( masses ) { }
60
operator ()solar_system_momentum61 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
center_of_mass(const container_type & x,const mass_type & m)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
energy(const container_type & q,const container_type & p,const mass_type & masses)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
streaming_observerstreaming_observer125 streaming_observer( std::ostream &out ) : m_out( out ) { }
126
127 template< class State >
operator ()streaming_observer128 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
main(int argc,char ** argv)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 */
202