#ifndef SHUFFLED_H
#define SHUFFLED_H

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

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

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

namespace MM
{

class Shuffled : public MTS
{
    int     inner_steps_2_;
    int     inner_steps_3_;
    int     update_list_;
    double  R0_, R1_, R2_;     // A
    double  skin1_, skin2_, skin3_;

public:
    explicit Shuffled (Prototype const &) {adopt_prototype (this);}
    explicit Shuffled () 
    :
        update_list_(5),
        inner_steps_2_(4), inner_steps_3_(2),
        R0_(1000.), R1_(4.5), R2_(8), 
        skin1_(0.5), skin2_(0.5), skin3_(0.5) 
    {inner_steps_ = 2;}


    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_;}
    void        set_inner_steps_3 (int v)   {inner_steps_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;}

protected:
    char const*     class_name () const {return "Shuffled";}
    char const*     title      () const {return "Shuffled";}
    char const*     alias      () const {return "Shuffled";}

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

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

        step_N_ = 0;
        list (MD_subject::Long3,    skin3_); 
        list (MD_subject::Long2,    R2_, skin2_); 
        list (MD_subject::Long1,    R1_, skin1_); 

        F (MD_subject::Long3); 
        F (MD_subject::Long2, R2_, S_);
        F (MD_subject::Long1, R1_, S_);
        F (MD_subject::Short);
    }

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

		int i, j, k;
        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_;

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

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

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

                    if (i == i_last_step)   
                        F (MD_subject::Long1, R1_, S_);
                    F (MD_subject::Short);

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

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

            F (MD_subject::Long2, R2_, S_);
            //F (MD_subject::Long1, R1_, S_);

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

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

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

                    if (i == i_last_step)   
                    {
                        //if (j == j_last_step) 
                        {
                            if (k == k_last_step)   
                            {
                                if (step_N_ % update_list_ == 0)
                                {
                                    list (MD_subject::Long3,    skin3_); 
                                    list (MD_subject::Long2,    R2_, skin2_); 
                                    list (MD_subject::Long1,    R1_, skin1_); 
                                }
                                zero_virial();
                                F (MD_subject::Long3);
                            }
                            F (MD_subject::Long2, R2_, S_);
                        }
                        F (MD_subject::Long1, R1_, S_);
                    }
                    F (MD_subject::Short);

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

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

        ++step_N_;
    }
};

}//MM

#endif //SHUFFLED_H
