odeint v2 - Solving ordinary differential equations in C++

/* Boost libs/numeric/odeint/examples/point_type.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 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