[Boost-commit] svn:boost r61327 - sandbox/odeint/libs/numeric/odeint/examples

Subject: [Boost-commit] svn:boost r61327 - sandbox/odeint/libs/numeric/odeint/examples
From: karsten.ahnert_at_[hidden]
Date: 2010-04-16 16:00:21

Author: karsten
Date: 2010-04-16 16:00:21 EDT (Fri, 16 Apr 2010)
New Revision: 61327
URL: http://svn.boost.org/trac/boost/changeset/61327

Log:
solar system example edited but still not finished
Added:
   sandbox/odeint/libs/numeric/odeint/examples/solar_system.hpp   (contents, props changed)
Text files modified: 
   sandbox/odeint/libs/numeric/odeint/examples/Jamfile          |     2                                         
   sandbox/odeint/libs/numeric/odeint/examples/solar_system.cpp |   104 ++++++--------------------------------- 
   2 files changed, 18 insertions(+), 88 deletions(-)

Modified: sandbox/odeint/libs/numeric/odeint/examples/Jamfile
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/Jamfile	(original)
+++ sandbox/odeint/libs/numeric/odeint/examples/Jamfile	2010-04-16 16:00:21 EDT (Fri, 16 Apr 2010)
@@ -26,5 +26,7 @@
 exe pendulum_vibrating_pivot : pendulum_vibrating_pivot.cpp ;
 exe dnls : dnls.cpp ;
 
+exe solar_system : solar_system.cpp ;
+
 exe doc_harm_osc : doc_harm_osc.cpp ;
 exe doc_integrate : doc_integrate.cpp ;

Modified: sandbox/odeint/libs/numeric/odeint/examples/solar_system.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/solar_system.cpp	(original)
+++ sandbox/odeint/libs/numeric/odeint/examples/solar_system.cpp	2010-04-16 16:00:21 EDT (Fri, 16 Apr 2010)
@@ -1,108 +1,36 @@
+/* Boost libs/numeric/odeint/examples/solar_system.cpp
+
+ Copyright 2009 Karsten Ahnert
+ Copyright 2009 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 <bits/stdc++.h>
 #include <bits/stdtr1c++.h>
 
-#include <boost/operators.hpp>
 #include <boost/circular_buffer.hpp>
 #include <boost/numeric/odeint.hpp>
 
+#include "solar_system.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 point< double ,3 > 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;
 
 

Added: sandbox/odeint/libs/numeric/odeint/examples/solar_system.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/libs/numeric/odeint/examples/solar_system.hpp	2010-04-16 16:00:21 EDT (Fri, 16 Apr 2010)
@@ -0,0 +1,166 @@
+/* Boost libs/numeric/odeint/examples/solar_system.hpp
+
+ Copyright 2009 Karsten Ahnert
+ Copyright 2009 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)
+*/
+
+#ifndef SOLAR_SYSTEM_HPP_INCLUDED
+#define SOLAR_SYSTEM_HPP_INCLUDED
+
+
+#include <boost/operators.hpp>
+#include <ostream>
+
+//
+// the point type
+//               
+template< class T , size_t Dim >
+class point :                   
+    boost::additive1< point< T , Dim > ,
+    boost::additive2< point< T , Dim  > , T ,
+    boost::multiplicative2< point< T , Dim > , T
+    > > >                                       
+{                                               
+public:                                         
+
+    const static size_t dim = Dim;
+    typedef T value_type;         
+    typedef point< value_type , dim > point_type;
+
+    point( void )
+        {        
+            for( size_t i=0 ; i<dim ; ++i ) m_val[i] = 0.0;
+        }                                                  
+    point( value_type val )                                
+        {                                                  
+            for( size_t i=0 ; i<dim ; ++i ) m_val[i] = val;
+        }                                                  
+    point( value_type x , value_type y , value_type z = 0.0 )
+        {                                                    
+            if( dim > 0 ) m_val[0] = x;                      
+            if( dim > 1 ) m_val[1] = y;                      
+            if( dim > 2 ) m_val[2] = z;                      
+        }                                                    
+
+    template< size_t i > T get( void ) const
+        {                                   
+            BOOST_STATIC_ASSERT( i < dim ); 
+            return m_val[i];                
+        }                                   
+    template< size_t i > T& get( void )     
+        {                                   
+            BOOST_STATIC_ASSERT( i < dim ); 
+            return m_val[i];                
+        }                                   
+
+    T operator[]( size_t i ) const { return m_val[i]; }
+    T& operator[]( size_t i ) { return m_val[i]; }     
+                                                       
+
+    point_type& operator+=( const point_type& p ) 
+        {                                         
+            for( size_t i=0 ; i<dim ; ++i )       
+                m_val[i] += p[i];                 
+            return *this;                         
+        }                                         
+
+    point_type& operator-=( const point_type& p )
+        {                                        
+            for( size_t i=0 ; i<dim ; ++i )      
+                m_val[i] -= p[i];                
+            return *this;                        
+        }                                        
+
+    point_type& operator+=( const value_type& val )
+        {                                          
+            for( size_t i=0 ; i<dim ; ++i )        
+                m_val[i] += val;                   
+            return *this;                          
+        }                                          
+
+    point_type& operator-=( const value_type& val )
+        {                                          
+            for( size_t i=0 ; i<dim ; ++i )        
+                m_val[i] -= val;                   
+            return *this;                          
+        }                                          
+
+    point_type& operator*=( const value_type &val )
+        {                                          
+            for( size_t i=0 ; i<dim ; ++i )        
+                m_val[i] *= val;                   
+            return *this;                          
+        }                                          
+
+    point_type& operator/=( const value_type &val )
+        {                                          
+            for( size_t i=0 ; i<dim ; ++i )        
+                m_val[i] /= val;                   
+            return *this;                          
+        }                                          
+
+private:
+
+    T m_val[dim];
+};               
+
+//
+// the - operator
+//               
+template< class T , size_t Dim >
+point< T , Dim > operator-( const point< T , Dim > &p )
+{                                                      
+    point< T , Dim > tmp;                              
+    for( size_t i=0 ; i<Dim ; ++i ) tmp[i] = -p[i];    
+    return tmp;                                        
+}                                                      
+
+//
+// scalar product
+//               
+template< class T , size_t Dim >
+T scalar_prod( const point< T , Dim > &p1 , const point< T , Dim > &p2 )
+{                                                                       
+    T tmp = 0.0;                                                        
+    for( size_t i=0 ; i<Dim ; ++i ) tmp += p1[i] * p2[i];               
+    return tmp;                                                         
+}                                                                       
+
+//
+// norm
+//
+template< class T , size_t Dim >
+T norm( const point< T , Dim > &p1 )
+{
+    return scalar_prod( p1 , p1 );
+}
+
+//
+// absolute value
+//
+template< class T , size_t Dim >
+T abs( const point< T , Dim > &p1 )
+{
+    return sqrt( norm( p1 ) );
+}
+
+//
+// output operator
+//
+template< class T , size_t Dim >
+std::ostream& operator<<( std::ostream &out , const point< T , Dim > &p )
+{
+    if( Dim > 0 ) out << p[0];
+    for( size_t i=1 ; i<Dim ; ++i ) out << " " << p[i];
+    return out;
+}
+
+
+
+#endif //SOLAR_SYSTEM_HPP_INCLUDED