/* Boost numeric test of the symplectic steppers test file

 Copyright 2012 Mario Mulansky
 Copyright 2012 Karsten Ahnert

 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)
*/

// disable checked iterator warning for msvc
#include <boost/config.hpp>
#ifdef BOOST_MSVC
    #pragma warning(disable:4996)
#endif

#define BOOST_TEST_MODULE numeric_symplectic

#include <iostream>
#include <cmath>

#include <boost/array.hpp>

#include <boost/test/unit_test.hpp>

#include <boost/mpl/vector.hpp>

#include <boost/numeric/odeint.hpp>

using namespace boost::unit_test;
using namespace boost::numeric::odeint;
namespace mpl = boost::mpl;

typedef double value_type;

typedef boost::array< double ,1 > state_type;

// harmonic oscillator, analytic solution x[0] = sin( t )
struct osc
{
    void operator()( const state_type &q , state_type &dpdt ) const
    {
        dpdt[0] = -q[0];
    }
};

BOOST_AUTO_TEST_SUITE( numeric_symplectic_test )


/* generic test for all symplectic steppers */
template< class Stepper >
struct perform_symplectic_test
{
    void operator()( void )
    {
   
        Stepper stepper;
        const int o = stepper.order()+1; //order of the error is order of approximation + 1

        const state_type q0 = {{ 0.0 }};
        const state_type p0 = {{ 1.0 }};
        state_type q1,p1;
        std::pair< state_type , state_type >x1( q1 , p1 );
        const double t = 0.0;
        /* do a first step with dt=0.1 to get an estimate on the prefactor of the error dx = f * dt^(order+1) */
        double dt = 0.5;
        stepper.do_step( osc() , std::make_pair( q0 , p0 ) , t , x1 , dt );
        const double f = 2.0 * std::abs( sin(dt) - x1.first[0] ) / std::pow( dt , o );

        std::cout << o << " , " << f << std::endl;

        /* as long as we have errors above machine precision */
        while( f*std::pow( dt , o ) > 1E-16 )
        {
            stepper.do_step( osc() , std::make_pair( q0 , p0 ) , t , x1 , dt );
            std::cout << "Testing dt=" << dt << std::endl;
            BOOST_CHECK_SMALL( std::abs( sin(dt) - x1.first[0] ) , f*std::pow( dt , o ) );
            dt *= 0.5;
        }
    }
};


typedef mpl::vector<
    symplectic_euler< state_type > ,
    symplectic_rkn_sb3a_mclachlan< state_type >
    > symplectic_steppers;

BOOST_AUTO_TEST_CASE_TEMPLATE( symplectic_test , Stepper, symplectic_steppers )
{
    perform_symplectic_test< Stepper > tester;
    tester();
}

BOOST_AUTO_TEST_SUITE_END()