/* 
   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_INTERACTION_H_
#define _ORSA_INTERACTION_H_

#include <vector>
#include <string>

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif // HAVE_CONFIG_H

#include "orsa_coord.h"
#include "orsa_common.h"
#include "orsa_body.h"
#include "orsa_error.h"
#include "orsa_frame.h"

namespace orsa {
  
  // class Frame;
  
  // == Unique id == if you need to modify this list, don't change the numbers, but add new numbers
  enum InteractionType {
    NEWTON=1,
    ARMONICOSCILLATOR=2,
    GALACTIC_POTENTIAL_ALLEN=3,
    GALACTIC_POTENTIAL_ALLEN_PLUS_NEWTON=4,
    JPL_PLANETS_NEWTON=5,
    GRAVITATIONALTREE=6,
    NEWTON_MPI=7,
    RELATIVISTIC=8
  };
  
  inline void convert(InteractionType &it, const unsigned int i)  {
    switch(i) {
    case 1: it = NEWTON;                               break;
    case 2: it = ARMONICOSCILLATOR;                    break;
    case 3: it = GALACTIC_POTENTIAL_ALLEN;             break;
    case 4: it = GALACTIC_POTENTIAL_ALLEN_PLUS_NEWTON; break;
    case 5: it = JPL_PLANETS_NEWTON;                   break;
    case 6: it = GRAVITATIONALTREE;                    break;
    case 7: it = NEWTON_MPI;                           break;
    case 8: it = RELATIVISTIC;                         break;
      //
    default: {
      char err_msg[1024];
      sprintf(err_msg,"conversion problem: i = %i",i);
      ORSA_ERROR_MACRO(err_msg);    
      break;       
    }
    }
  }
  
  std::string label(const InteractionType it);
  
  // bool depends_on_velocities(const InteractionType it);
  
  // TODO: distinguish between conservative and non-conservative interactions
  
  // base class Interaction
  
  class Interaction {
  public:
    virtual ~Interaction() { };
    
  public:
    virtual void Acceleration(const Frame&, std::vector<Vector>&) = 0;  
    virtual double PotentialEnergy(const Frame&) = 0;
    
  public:
    virtual Interaction * clone() const = 0;
    
  public:
    virtual bool depends_on_velocity() const { return false; }
    
    //
    // protected:
    // InteractionType type;
  public:	
    virtual InteractionType GetType() const = 0;
  };
  
  void make_new_interaction(Interaction**, const InteractionType);
  
  // derived classes
  
  class Newton : public Interaction {
  public:
    Newton();
    Newton(const Newton &);
    
  public:
    void Acceleration(const Frame&, std::vector<Vector>&);
    double PotentialEnergy(const Frame&);
    
  public:
    Interaction * clone() const;
    
    InteractionType GetType() const {
      return NEWTON;
    }
    
  private:
    const double g;
  };
  
#ifdef HAVE_MPI
  class Newton_MPI : public Interaction {
  public:
    Newton_MPI();
    Newton_MPI(const Newton_MPI &);
    ~Newton_MPI();
    
  public:
    void Acceleration(const Frame&, std::vector<Vector>&);
    double PotentialEnergy(const Frame&);
    
  public:
    Interaction * clone() const;
    
    InteractionType GetType() const {
      return NEWTON_MPI;
    }
    
  private:
    double g;
    //
    double *bv;
    double *av;
    double *av_local;
  };
#endif // HAVE_MPI
  
  class GravitationalTree : public Interaction {
  public:
    GravitationalTree();
    GravitationalTree(const GravitationalTree &);
    
  public:
    void Acceleration(const Frame&, std::vector<Vector>&);
    double PotentialEnergy(const Frame&);
  
  public:
    Interaction * clone() const;
    
    InteractionType GetType() const {
      return GRAVITATIONALTREE;
    }
    
  private:
    double g;
    double theta;
  };
  
  class Relativistic : public Interaction {
  public:
    Relativistic();
    Relativistic(const Relativistic &);
    
  public:
    void Acceleration(const Frame&, std::vector<Vector>&);
    double PotentialEnergy(const Frame&);
    
  public:
    Interaction * clone() const;
    
    InteractionType GetType() const {
      return RELATIVISTIC;
    }
    
  public:
    bool depends_on_velocity() const { return true; }
    
  private:
    const double g;
    const double c_2;
  };
  
  class ArmonicOscillator : public Interaction {
  public:
    ArmonicOscillator(const double free_length_in, const double k_in);
    ArmonicOscillator(const ArmonicOscillator &);
    
  public:
    void Acceleration(const Frame&, std::vector<Vector>&);
    double PotentialEnergy(const Frame&);
    
  public:
    Interaction * clone() const;
    
    InteractionType GetType() const {
      return ARMONICOSCILLATOR;
    }
    
  private:
    const double free_length;
    const double k; // 'spring' constant
  };
  
  class GalacticPotentialAllen : public Interaction {
  public:
    GalacticPotentialAllen();
    GalacticPotentialAllen(const GalacticPotentialAllen &);
    
  public:
    void Acceleration(const Frame&, std::vector<Vector>&);
    double PotentialEnergy(const Frame&);
    
  public:
    Interaction * clone() const;
    
    InteractionType GetType() const {
      return GALACTIC_POTENTIAL_ALLEN;
    }
    
  private:
    double g;
    double mb; // bulge mass
    double bb; // bulge length scale
    double md; // disk mass
    double ad; // disk length scale
    double bd; // disk height scale
    double mh; // halo mass (different from the reference article, rescaled to obtain a circular velocity at the sun radius of about 220 km/s)
    double ah; // halo length scale
  };
  
  class GalacticPotentialAllenPlusNewton : public Interaction {
  public:
    GalacticPotentialAllenPlusNewton() : Interaction() {
      
    }
      
    GalacticPotentialAllenPlusNewton(const GalacticPotentialAllenPlusNewton &) : Interaction() {
      
    }
    
    InteractionType GetType() const {
      return GALACTIC_POTENTIAL_ALLEN_PLUS_NEWTON;
    }
    
  public:
    Interaction * clone() const {
      return new GalacticPotentialAllenPlusNewton(*this);
    }
    
    inline void Acceleration(const Frame &f, std::vector<Vector> &a) {
      
      tmp_a.resize(a.size());
      
      unsigned int i;
      
      for (i=0;i<a.size();++i)
	a[i].Set(0,0,0);
      
      gpa_itg.Acceleration(f,tmp_a);
      for (i=0;i<a.size();++i)
	a[i] += tmp_a[i];
      
      newton_itg.Acceleration(f,tmp_a);
      for (i=0;i<a.size();++i)
	a[i] += tmp_a[i];
      
    }	
    
    inline double PotentialEnergy(const Frame &f) {
      return (gpa_itg.PotentialEnergy(f)+newton_itg.PotentialEnergy(f));
    }
    
  private:
    GalacticPotentialAllen    gpa_itg;
    Newton                 newton_itg;
    //
    std::vector<Vector> tmp_a;
  };
  
  class JPLPlanetsNewton : public Interaction {
  public:
    JPLPlanetsNewton(std::list<JPL_planets> &);
    JPLPlanetsNewton(const JPLPlanetsNewton &);
    
  public:
    void Acceleration(const orsa::Frame &, std::vector<orsa::Vector>&);
    double PotentialEnergy(const orsa::Frame &);
    
  public:
    Interaction * clone() const;
    
    InteractionType GetType() const {
      return JPL_PLANETS_NEWTON;
    }
    
  private:
    Newton newton_itg;
    std::list<JPL_planets> l;
    orsa::Frame planets_frame;
    double g;
  };
  
} // namespace orsa

#endif // _ORSA_INTERACTION_H_
