/* 
   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:
      ORSA_ERROR("conversion problem: i = %i",i);
      break;       
    }
  }
  
  std::string label(const InteractionType it);
  
  // bool depends_on_velocities(const InteractionType it);
  
  // 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; }
    
  public:
    void SkipJPLPlanets(const bool b) {
      skip_JPL_planets = b;
    }
    bool IsSkippingJPLPlanets() const {
      return skip_JPL_planets;
    }
  protected:
    bool skip_JPL_planets;
    
  public:	
    virtual InteractionType GetType() const = 0;
  };
  
  void make_new_interaction(Interaction**, const InteractionType);
  
  // 
  
  class MappedTable {
  public:
    void load(const std::vector<Body> & f, const bool skip_JPL_planets);
    
  public:
    inline Vector DistanceVector(const unsigned int i, const unsigned int j) const {
      return (sign_ij(i,j)*distance_vector[ij_to_index(i,j)]);
    }
    
    inline double Distance(const unsigned int i, const unsigned int j) const {
      return d1[ij_to_index(i,j)];
    }
    
    inline double Distance2(const unsigned int i, const unsigned int j) const {
      return d2[ij_to_index(i,j)];
    }
    
    inline double Distance3(const unsigned int i, const unsigned int j) const {
      return d3[ij_to_index(i,j)];
    }
    
    inline double Distance4(const unsigned int i, const unsigned int j) const {
      return d4[ij_to_index(i,j)];
    }
    
    inline double OneOverDistance(const unsigned int i, const unsigned int j) const {
      return one_over_distance[ij_to_index(i,j)];
    }
    
    inline double OneOverDistanceSquare(const unsigned int i, const unsigned int j) const {
      return one_over_distance_square[ij_to_index(i,j)];
    }
    
    inline double OneOverDistanceCube(const unsigned int i, const unsigned int j) const {
      return one_over_distance_cube[ij_to_index(i,j)];
    }
    
    inline Vector DistanceVectorOverDistanceCube(const unsigned int i, const unsigned int j) const {
      return (sign_ij(i,j)*distance_vector_over_distance_cube[ij_to_index(i,j)]);
    }
    
  private:
    // i,j are in the frame domain,
    // while index is in the mapping domain
    
    inline unsigned int ij_to_index(const unsigned int i, const unsigned int j) const {
      if (swap_ij(i,j)) {
	return (mapping[j]*M+mapping[i]);
      } else {
	return (mapping[i]*M+mapping[j]);
      }
    }
    
    inline bool swap_ij(const unsigned int i, const unsigned int j) const {
      return (i < j);
    }
    
    inline double sign_ij(const unsigned int i, const unsigned int j) const {
      if (swap_ij(i,j)) {
	return -1.0;
      } else {
	return +1.0;
      }
    }
    
    inline void index_to_ij(const unsigned int index, unsigned int & i, unsigned int & j) const {
      const unsigned int mi = index / M;
      const unsigned int mj = index % M;
      //
      bool found_i = false;
      bool found_j = false;
      for (unsigned int k=0; k<N; ++k) {
	if (mapping[k] == mi) {
	  i = k;
	  found_i = true;
	}
	if (mapping[k] == mj) {
	  j = k;
	  found_j = true;
	}
	if (found_i && found_j) break;
      }
    }
    
  private:
    // M < N, usually is the number of massive particles
    unsigned int M, N, MN;
    //
    std::vector<unsigned int> mapping;
    std::vector<Vector> distance_vector;
    std::vector<Vector> distance_vector_over_distance_cube;
    std::vector<double> d1;
    std::vector<double> d2;
    std::vector<double> d3;
    std::vector<double> d4;
    std::vector<double> one_over_distance;
    std::vector<double> one_over_distance_square;
    std::vector<double> one_over_distance_cube;
  };
  
  class Legendre {
  public:
    Legendre(const double arg) :
      x(arg), x2(x*x), x3(x2*x), x4(x3*x), tmp1(1.0-x2), tmp2(std::sqrt(tmp1)), tmp3(-x/tmp2),
      P2(0.5*(3.0*x2-1.0)),dP2(3.0*x),
      P3(0.5*(5.0*x3-3.0*x)),dP3(0.5*(15.0*x2-3.0)),
      P4(0.125*(35.0*x4-30.0*x2+3.0)),dP4(0.125*(140.0*x3-60.0*x)),
      P22(3.0*tmp1),dP22(-6.0*x),
      P31(1.5*(1.0-5.0*x2)*tmp2),dP31(1.5*(-10.0*x*tmp2+(1.0-5.0*x2)*tmp3)),
      P32(15.0*x*tmp1),dP32(15.0-45.0*x2),
      P33(-15.0*tmp1*tmp2),dP33(-15.0*(-2.0*x*tmp2+tmp1*tmp3)),
      P41(2.5*x*(3.0-7.0*x2)*tmp2),dP41(2.5*((3.0-21.0*x2)*tmp2+x*(3.0-7.0*x2)*tmp3)),
      P42(7.5*(7.0*x2-1.0)*tmp1),dP42(7.5*(14.0*x*tmp1-14.0*x3+2.0*x)),
      P43(-105.0*x*tmp1*tmp2),dP43(-105.0*(tmp1*tmp2-2*x2*tmp2+x*tmp1*tmp3)),
      P44(105.0*tmp1*tmp1),dP44(-420.0*x*tmp1)
      { }
      
  private:
    const double x,x2,x3,x4;
  private:
    const double tmp1, tmp2, tmp3;
  public:
    const double P2,dP2,P3,dP3,P4,dP4;
    const double P22,dP22,P31,dP31,P32,dP32,P33,dP33,P41,dP41,P42,dP42,P43,dP43,P44,dP44;
  };
  
  // derived classes
  
  class Newton : public Interaction {
  public:
    Newton();
    Newton(const Newton &);
    
  public:
    void Acceleration(const Frame&, std::vector<Vector>&);
    double PotentialEnergy(const Frame&);
    
  public:
    inline bool depends_on_velocity() const { 
      return (include_relativistic_effects || include_fast_relativistic_effects); 
    }
    
  public:
    Interaction * clone() const;
    
    InteractionType GetType() const {
      return NEWTON;
    }
    
  private: 
    void fast_newton_acc(const Frame &, std::vector<Vector> &);    
    
  private:
    std::vector<bool> zero_mass;
    std::vector<Vector> a_newton;
    //
    MappedTable mapped_table;
    
  public:
    void IncludeMultipoleMoments(const bool b) {
      include_multipole_moments = b;
    }
    bool IsIncludingMultipoleMoments() const {
      return include_multipole_moments;
    }
  private:
    bool include_multipole_moments;
    //
    std::vector<Vector> a_multipoles;
    //
    // tmp stuff
    std::vector<Vector> axis, x_axis; // axis is the z axis
    // std::vector<double> phi; // meridian angle
    std::vector<double> R1;
    std::vector<double> R2;
    std::vector<double> R3;
    std::vector<double> R4;
    
  public:
    void IncludeRelativisticEffects(const bool b) {
      include_relativistic_effects = b;
    }
    bool IsIncludingRelativisticEffects() const {
      return include_relativistic_effects;
    }
  private:
    bool include_relativistic_effects;
  public:
    void IncludeFastRelativisticEffects(const bool b) {
      include_fast_relativistic_effects = b;
    }
    bool IsIncludingFastRelativisticEffects() const {
      return include_fast_relativistic_effects;
    }
  private:
    bool include_fast_relativistic_effects;
    //
    const double one_over_c2;
    //
    std::vector<Vector> a_relativity;
    
    /* 
       public:
       void SkipJPLPlanets(const bool b) {
       skip_JPL_planets = b;
       }
       bool IsSkippingJPLPlanets() const {
       return skip_JPL_planets;
       }
       private:
       bool skip_JPL_planets;
    */
    
  private:
    std::vector<bool> skip;
  };
  
#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_
