#ifndef RRESPA_4_H
#define RRESPA_4_H

#ifndef MTS_H
#include "MTS.h"
#endif

#ifndef MD_H
#include "MD.h"
#endif

#ifndef LOCK_H
#include "Lock.h"
#endif  

#ifndef TIMER_H
#include "Timer.h"
#define TIME(ref) ;
#endif  

namespace MM
{


class rRESPA_4 : public MTS
{
protected:
    int     inner_steps_2_;
    int     inner_steps_3_;
    int     /*update_list_1_,*/ update_list_2_, update_list_3_;
    double  /*R0_, R1_,*/ R2_;     // A
    double  /*skin1_,*/ skin2_, skin3_;

public:
    explicit rRESPA_4 (Prototype const &) {adopt_prototype (this);}
    explicit rRESPA_4 () 
    :   
    // precise
        //update_list_(10),
        //inner_steps_2_(4), inner_steps_3_(2),
        //R0_(1000.), R1_(4), R2_(8), 
        ////skin1_(1.5), skin2_(1.5), skin3_(10) 

    // fast
        //update_list_(5),
        //inner_steps_2_(4), inner_steps_3_(2),
        //R0_(1000.), R1_(4), R2_(8), 
        //skin1_(0.5), skin2_(0.5), skin3_(0.5) 

        //update_list_1_(4), update_list_2_(4), update_list_3_(4),
        //inner_steps_2_(4), inner_steps_3_(2),
        //R0_(1000.), R1_(4.5), R2_(8), 
        ////R0_(1000.), R1_(6), R2_(15), 
        //skin1_(0.5), skin2_(0.5), skin3_(0.5) 

        //update_list_2_(4), update_list_3_(4),
        //inner_steps_2_(4), inner_steps_3_(2),
        //R2_(8), 
        ////R0_(1000.), R1_(6), R2_(15), 
        //skin2_(0.5), skin3_(0.5) 

        update_list_2_(8), update_list_3_(16),
        inner_steps_2_(2), inner_steps_3_(2),
        R2_(8), 
        //R0_(1000.), R1_(6), R2_(15), 
        skin2_(0.7), skin3_(1.) 

    {
        update_list_ = 4;
        inner_steps_ = 2; 
        R_ = 4.5; 
        skin_ = 0.5;
    }

    int         inner_steps_2 ()  const     {return inner_steps_2_;}
    void        set_inner_steps_2 (int v)   {inner_steps_2_= v;}

    int         inner_steps_3 ()  const     {return inner_steps_3_;}
    void        set_inner_steps_3 (int v)   {inner_steps_3_= v;}

    int         update_list_2 ()  const     {return update_list_2_;}
    void        set_update_list_2 (int v)   {update_list_2_= v;}

    int         update_list_3 ()  const     {return update_list_3_;}
    void        set_update_list_3 (int v)   {update_list_3_= v;}


    //double      R0 ()  const                {return R0_;}
    //void        set_R0 (double v)           {R0_= v;}

    //double      R1 ()  const                {return R1_;}
    //void        set_R1 (double v)           {R1_= v;}

    double      R2 ()  const                {return R2_;}
    void        set_R2 (double v)           {R2_= v;}


    //double      skin_1 ()  const            {return skin1_;}
    //void        set_skin_1 (double v)       {skin1_= v;}

    double      skin_2 ()  const            {return skin2_;}
    void        set_skin_2 (double v)       {skin2_= v;}

    double      skin_3 ()  const            {return skin3_;}
    void        set_skin_3 (double v)       {skin3_= v;}

protected:
    char const*     class_name () const {return "rRESPA_4";}
    char const*     title      () const {return "rRESPA 4-level";}
    char const*     alias      () const {return "rRESPA 4";}

    MD_method *     clone () const      {return new rRESPA_4;}

    virtual void start ()
    {
        Lock<MD> lock(simulation());
        MD_method::start ();
        TIME ("      ");
        step_N_ = 0;

        step_N_ = 0;
        list (MD_subject::Long3,    skin3_); 
        TIME ("list 3");
        list (MD_subject::Long2,    R2_, skin2_); 
        TIME ("list 2");
        list (MD_subject::Long1,    R_, skin_); 
        TIME ("list 1");

        //F (MD_subject::Whole);
        //TIME ("F W   ");
        F (MD_subject::Long3); 
        TIME ("F 3   ");
        F (MD_subject::Long2, R2_, S_);
        TIME ("F 2   ");
        F (MD_subject::Long1, R_, S_);
        TIME ("F 1   ");
        F (MD_subject::Short);
        TIME ("F 0   ");
    }


    virtual void step ()
    {
        TIME (" { ### ");
        Lock<MD> lock(simulation());
        TIME (" lock  ");

        int i_last_step = inner_steps_  -1;
        int j_last_step = inner_steps_2_-1;
        int k_last_step = inner_steps_3_-1;
        int i_steps = inner_steps_  * inner_steps_2_ * inner_steps_3_;
        int j_steps = inner_steps_2_* inner_steps_3_;
        int k_steps = inner_steps_3_;

        //if (step_N_ % update_list_ == 0)
        //    list (MD_subject::Long3,    skin3_); 

        kick (.5, MD_subject::Long3);
        simulation_->algorithm().description ("{");

        //if (step_N_ % update_list_ == 0)
        //    list (MD_subject::Long2,    R2_, skin2_); 

        for (int k=0;  k<inner_steps_3_;  ++k)
        {
            kick (.5 / k_steps, MD_subject::Long2);
            simulation_->algorithm().description ("  {");

            //if (step_N_ % update_list_ == 0)
            //    list (MD_subject::Long1,    R1_, skin1_); 

            for (int j=0;  j<inner_steps_2_;  ++j)
            {
                kick (.5 / j_steps, MD_subject::Long1);
                simulation_->algorithm().description ("    {");

                for (int i=0;  i<inner_steps_;  ++i)
                {
                    kick  (.5 / i_steps, MD_subject::Short);
                    TIME ("kick  .5 ");
                    drift (1. / i_steps);
                    TIME ("drift 1. ");

                    if (i == i_last_step)
                    {
                        if (j == j_last_step) 
                        {
                            if (k == k_last_step)
                            {
                                if (step_N_ % update_list_3_ == 0)
                                    list (MD_subject::Long3,    skin3_); 
                                if (step_N_ % update_list_2_ == 0)
                                    list (MD_subject::Long2,    R2_, skin2_); 
                                if (step_N_ % update_list_   == 0)
                                    list (MD_subject::Long1,    R_, skin_); 

                                zero_virial();
                                F (MD_subject::Long3);
                                TIME ("F     3");
                            }
                            F (MD_subject::Long2, R2_, S_);
                            TIME ("F     2");
                        }
                        F (MD_subject::Long1, R_, S_);
                        TIME ("F     1");
                    }

                    F (MD_subject::Short);
                    TIME ("F     0");

                    kick  (.5 / i_steps, MD_subject::Short);
                    TIME ("kick  .5 ");

                    if (i<i_last_step)
                        simulation_->algorithm().description ("");
                }
                simulation_->algorithm().description ("    }");
                kick (.5 / j_steps, MD_subject::Long1);
            }
            simulation_->algorithm().description ("  }");
            kick (.5 / k_steps, MD_subject::Long2);
        }
        simulation_->algorithm().description ("}");
        kick (.5, MD_subject::Long3);

        ++step_N_;

        TIME (" } ### ");
    }

    //virtual void step ()
    //{
    //    inner_steps_ = 10;

    //    int last_step = inner_steps_  -1;

    //    kick (.5, MD_subject::Long2);
    //    simulation_->algorithm().description ("{");

    //    for (int i=0;  i<inner_steps_;  ++i)
    //    {
    //        kick  (.5 / inner_steps_, MD_subject::Long);
    //        kick  (.5 / inner_steps_, MD_subject::Short);
    //        drift (1. / inner_steps_);

    //        if (i == last_step) 
    //            F (MD_subject::Long2, 1000);

    //        F (MD_subject::Long, 6);
    //        F (MD_subject::Short, R0_);

    //        kick  (.5 / inner_steps_, MD_subject::Short);
    //        kick  (.5 / inner_steps_, MD_subject::Long);

    //        if (i<last_step)
    //            simulation_->algorithm().description ("");
    //    }
    //    simulation_->algorithm().description ("}");
    //    kick (.5, MD_subject::Long2);
    //}
};

}//MM

#endif //RRESPA_4_H
