/* 
   ORSA - Orbit Reconstruction, Simulation and Analysis
   Copyright (C) 2002-2004 Pasquale Tricarico
   
   This program is free software; you can redistribute it and/or
   modify it under the terms of the GNU General Public License
   as published by the Free Software Foundation; either version 2
   of the License, or (at your option) any later version.
   
   As a special exception, Pasquale Tricarico gives permission to
   link this program with Qt commercial edition, and distribute the
   resulting executable, without including the source code for the Qt
   commercial edition in the source distribution.
   
   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
   
   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*/

#ifndef _ORSA_INTEGRATOR_H_
#define _ORSA_INTEGRATOR_H_

#include "orsa_frame.h"
#include "orsa_error.h"
#include "orsa_interaction.h"

#include <string>

namespace orsa {
  
  enum IntegratorType {
    STOER=1,
    BULIRSCHSTOER=2,
    RUNGEKUTTA=3,
    DISSIPATIVERUNGEKUTTA=4,
    RA15=5,
    LEAPFROG=6
  };
  
  inline void convert(IntegratorType &it, const unsigned int i)  {
    switch(i) {
    case 1: it = STOER;                 break;
    case 2: it = BULIRSCHSTOER;         break;
    case 3: it = RUNGEKUTTA;            break;
    case 4: it = DISSIPATIVERUNGEKUTTA; break;
    case 5: it = RA15;                  break;
    case 6: it = LEAPFROG;              break;
      //
    default:
      ORSA_ERROR("conversion problem: i = %i",i);
      break;       
    }
  }
  
  inline std::string label(const IntegratorType it) {
    std::string s;
    switch (it) {
    case STOER:	                 s="Stoer";                              break;
    case BULIRSCHSTOER:          s="Bulirsch-Stoer";                     break;
    case RUNGEKUTTA:             s="Runge-Kutta 4th order";              break;
    case DISSIPATIVERUNGEKUTTA:  s="Dissipative Runge-Kutta 4th order";  break;
    case RA15:                   s="Everhart's RADAU 15th order";        break;
    case LEAPFROG:               s="Leapfrog 2nd order";                 break;
    }	
    return s;
  }
  
  //! This is the interface for all the Integrator classes.
  class Integrator {
  public:
    virtual void Step(const Frame &, Frame &, Interaction *) = 0;
    virtual ~Integrator() { };
    
  public:
    virtual Integrator * clone() const = 0;
    
  public:
    // double timestep,timestep_done;
    UniverseTypeAwareTimeStep timestep;
    
  protected:
    UniverseTypeAwareTimeStep timestep_done;
    
  public:
    // these vars should be moved inside each single class
    double accuracy; //! used only with variable step size integrators
    unsigned int m;  //! substeps for multisteps integrators
    
  public:
    virtual bool can_handle_velocity_dependant_interactions() const { return false; }
    
  public:
    inline IntegratorType GetType() const { return type; }
    
  protected: 
    IntegratorType type;
  };
  
  void make_new_integrator(Integrator**, const IntegratorType);
  
  // 
  // First classification: fixed timestep and variable timestep
  //
  
  class FixedTimestepIntegrator : public Integrator {
    
  };
  
  class MultistepIntegrator : public FixedTimestepIntegrator {
    
  };
  
  class VariableTimestepIntegrator : public Integrator {
    
  };
  
  //
  // derived classes
  //
  
  //! Advances using a number of substeps (midpoints). For conservative and non-conservative Interaction.
  class ModifiedMidpoint : public MultistepIntegrator {
    
  };
  
  //! Like the ModifiedMidpoint, but for conservative Interaction.
  class Stoer : public MultistepIntegrator {
  public:
    Stoer();
    Stoer(int);
    Stoer(const Stoer &);
    ~Stoer();
    
  public:
    void Step(const Frame&, Frame&, Interaction*);
    
  public:
    Integrator * clone() const;
  };
  
  /* 
     class BulirschStoer : public VariableTimestepIntegrator {
     public:
     BulirschStoer();
     // BulirschStoer(Interaction*);
     
     ~BulirschStoer();
     
     private:
     void init();
     
     public:
     void Step(const Frame&, Frame&, Interaction*);
     
     
     // for internal computations only
     private:
     double h, h_done, h_next, h_try;
     Frame main_frame, frame_save, frame_seq, frame_out_stoerm;
     std::vector<unsigned int> midpoint_sequence; // {2,4,6,8,10,12,14,16,18,20};
     const unsigned int max_counter_midpoint; // = 8; // ==> midpoint_sequence[]={2,4,...,14,16}
     const unsigned int max_n_points; //= max_counter_midpoint;  // to delete in future
     Stoer *stoer;
     //
     bool flag_first;
     unsigned int k_optimum, k_max;
     double old_accuracy;
     std::vector<double> a;
     std::vector< std::vector <double> > alpha;
     double x_new;
     };
  */
  
  class RungeKutta : public FixedTimestepIntegrator {
   public:
    RungeKutta();
    RungeKutta(const RungeKutta &);
    ~RungeKutta();
    
  public:
    void Step(const Frame&, Frame&, Interaction*);
    
  public:
    Integrator * clone() const;
  };
  
  class DissipativeRungeKutta : public FixedTimestepIntegrator {
  public:
    DissipativeRungeKutta();
    DissipativeRungeKutta(const DissipativeRungeKutta &);
    ~DissipativeRungeKutta();
    
  public:
    bool can_handle_velocity_dependant_interactions() const { return true; }
    
  public:
    void Step(const Frame&, Frame&, Interaction*);
    
  public:
    Integrator * clone() const;
  };
  
  /* 
     class Symplectic : public FixedTimestepIntegrator {
     
     };
  */
  
  class Radau15 : public VariableTimestepIntegrator {
  public:
    Radau15();
    Radau15(const Radau15 &);
    ~Radau15();
    
  private:
    void init();
    void Bodies_Mass_or_N_Bodies_Changed(const Frame&);
    
  public:
    bool can_handle_velocity_dependant_interactions() const { return true; }
    
  public:
    void Step(const Frame&, Frame&, Interaction*);
    
  public:
    Integrator * clone() const;
    
  private:
    // Gauss-Radau spacings for substeps within a sequence, for the 15th order 
    // integrator. The sum of the H values should be 3.733333333333333 
    // const vector<double> h;
    // vector<double> h;
    double h[8];
    
    // Constant coefficients used in series expansions for X and V
    //  XC: 1/2,  1/6,  1/12, 1/20, 1/30, 1/42, 1/56, 1/72
    //  VC: 1/2,  1/3,  1/4,  1/5,  1/6,  1/7,  1/8
    // vector<double> xc,vc;
    double xc[8],vc[7];
    
    // vector<double> r,c,d,s;
    double r[28],c[21],d[21],s[9];
    
    std::vector< std::vector<double> > g,b,e;
    // double g[7][3000],b[7][3000],e[7][3000];
    
    unsigned int nv,niter;
    
    std::vector<double> x,v,a,x1,v1,a1;
    // double x[3000],v[3000],a[3000],x1[3000],v1[3000],a1[3000];
    
    std::vector<double> mass;
    // double mass[1000];
    
    std::vector<Vector> acc;
    // Vector acc[1000];
    
    unsigned int size;
  };  
  
  class Leapfrog : public FixedTimestepIntegrator {
  public:
    Leapfrog(); 
    Leapfrog(const Leapfrog &);
    ~Leapfrog();
    
  public:
    void Step(const Frame&, Frame&, Interaction*);
    
  public:
    Integrator * clone() const;
  };
  
} // namespace orsa

#endif // _ORSA_INTEGRATOR_H_

