/* Boost libs/numeric/odeint/examples/point_type.hpp 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) */ #ifndef POINT_TYPE_HPP_INCLUDED #define POINT_TYPE_HPP_INCLUDED #include <boost/operators.hpp> #include <ostream> //[ point_type /*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; // ... // constructors //<- 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; } //-> // ... // operators //<- 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]; }; //... // more operators //] // // 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 //POINT_TYPE_HPP_INCLUDED