/*
 * stepper_details.cpp
 *
 * This example demonstrates some details about the steppers in odeint.
 *
 *
 * Copyright 2011-2012 Karsten Ahnert
 * Copyright 2012 Mario Mulansky
 * Copyright 2013 Pascal Germroth
 *
 * 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 <iostream>
#include <boost/array.hpp>
#include <boost/bind.hpp>
#include <boost/numeric/odeint.hpp>

using namespace std;
using namespace boost::numeric::odeint;

const size_t N = 3;

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

//[ system_function_structure
void sys( const state_type & /*x*/ , state_type & /*dxdt*/ , const double /*t*/ )
{
    // ...
}
//]

void sys1( const state_type &/*x*/ , state_type &/*dxdt*/ , const double /*t*/ )
{
}

void sys2( const state_type &/*x*/ , state_type &/*dxdt*/ , const double /*t*/ )
{
}


//[ symplectic_stepper_detail_system_function
typedef boost::array< double , 1 > vector_type;


struct harm_osc_f1
{
    void operator()( const vector_type &p , vector_type &dqdt )
    {
        dqdt[0] = p[0];
    }
};

struct harm_osc_f2
{
    void operator()( const vector_type &q , vector_type &dpdt )
    {
        dpdt[0] = -q[0];
    }
};
//]

//[ symplectic_stepper_detail_system_class
struct harm_osc
{
    void f1( const vector_type &p , vector_type &dqdt ) const
    {
        dqdt[0] = p[0];
    }

    void f2( const vector_type &q , vector_type &dpdt ) const
    {
        dpdt[0] = -q[0];
    }
};
//]

int main( int argc , char **argv )
{
    using namespace std;

    // Explicit stepper example
    {
        double t( 0.0 ) , dt( 0.1 );
        state_type in , out , dxdtin , inout;
        //[ explicit_stepper_detail_example
        runge_kutta4< state_type > rk;
        rk.do_step( sys1 , inout , t , dt );               // In-place transformation of inout
        rk.do_step( sys2 , inout , t , dt );               // call with different system: Ok
        rk.do_step( sys1 , in , t , out , dt );            // Out-of-place transformation
        rk.do_step( sys1 , inout , dxdtin , t , dt );      // In-place tranformation of inout
        rk.do_step( sys1 , in , dxdtin , t , out , dt );   // Out-of-place transformation
        //]
    }


    // FSAL stepper example
    {
        double t( 0.0 ) , dt( 0.1 );
        state_type in , in2 , in3 , out , dxdtin , dxdtout , inout , dxdtinout;
        //[ fsal_stepper_detail_example
        runge_kutta_dopri5< state_type > rk;
        rk.do_step( sys1 , in , t , out , dt );
        rk.do_step( sys2 , in , t , out , dt );         // DONT do this, sys1 is assumed

        rk.do_step( sys2 , in2 , t , out , dt );
        rk.do_step( sys2 , in3 , t , out , dt );        // DONT do this, in2 is assumed

        rk.do_step( sys1 , inout , dxdtinout , t , dt );
        rk.do_step( sys2 , inout , dxdtinout , t , dt );           // Ok, internal derivative is not used, dxdtinout is updated

        rk.do_step( sys1 , in , dxdtin , t , out , dxdtout , dt );
        rk.do_step( sys2 , in , dxdtin , t , out , dxdtout , dt ); // Ok, internal derivative is not used
        //]
    }


    // Symplectic harmonic oscillator example
    {
        double t( 0.0 ) , dt( 0.1 );
        //[ symplectic_stepper_detail_example
        pair< vector_type , vector_type > x;
        x.first[0] = 1.0; x.second[0] = 0.0;
        symplectic_rkn_sb3a_mclachlan< vector_type > rkn;
        rkn.do_step( make_pair( harm_osc_f1() , harm_osc_f2() ) , x , t , dt );
        //]

        //[ symplectic_stepper_detail_system_class_example
        harm_osc h;
        rkn.do_step( make_pair( boost::bind( &harm_osc::f1 , h , _1 , _2 ) , boost::bind( &harm_osc::f2 , h , _1 , _2 ) ) ,
                x , t , dt );
        //]
    }

    // Simplified harmonic oscillator example
    {
        double t = 0.0, dt = 0.1;
        //[ simplified_symplectic_stepper_example
        pair< vector_type , vector_type > x;
        x.first[0] = 1.0; x.second[0] = 0.0;
        symplectic_rkn_sb3a_mclachlan< vector_type > rkn;
        rkn.do_step( harm_osc_f1() , x , t , dt );
        //]

        vector_type q = {{ 1.0 }} , p = {{ 0.0 }};
        //[ symplectic_stepper_detail_ref_usage
        rkn.do_step( harm_osc_f1() , make_pair( boost::ref( q ) , boost::ref( p ) ) , t , dt );
        rkn.do_step( harm_osc_f1() , q , p , t , dt );
        rkn.do_step( make_pair( harm_osc_f1() , harm_osc_f2() ) , q , p , t , dt );
        //]
    }
    
    // adams_bashforth_moulton stepper example
    {
        double t = 0.0 , dt = 0.1;
        state_type inout;
        //[ multistep_detail_example
        adams_bashforth_moulton< 5 , state_type > abm;
        abm.initialize( sys , inout , t , dt );
        abm.do_step( sys , inout , t , dt );
        //]
        
        //[ multistep_detail_own_stepper_initialization
        abm.initialize( runge_kutta_fehlberg78< state_type >() , sys , inout , t , dt );
        //]
    }


    // dense output stepper examples
    {
        double t = 0.0 , dt = 0.1;
        state_type in;
        //[ dense_output_detail_example
        dense_output_runge_kutta< controlled_runge_kutta< runge_kutta_dopri5< state_type > > > dense;
        dense.initialize( in , t , dt );
        pair< double , double > times = dense.do_step( sys );
        (void)times;
        //]

        state_type inout;
        double t_start = 0.0 , t_end = 1.0;
        //[ dense_output_detail_generation1
        typedef boost::numeric::odeint::result_of::make_dense_output<
            runge_kutta_dopri5< state_type > >::type dense_stepper_type;
        dense_stepper_type dense2 = make_dense_output( 1.0e-6 , 1.0e-6 , runge_kutta_dopri5< state_type >() );
        (void)dense2;
        //]

        //[ dense_output_detail_generation2
        integrate_const( make_dense_output( 1.0e-6 , 1.0e-6 , runge_kutta_dopri5< state_type >() ) , sys , inout , t_start , t_end , dt );
        //]
    }


    return 0;
}