/* 
   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_ORBIT_H_
#define _ORSA_ORBIT_H_

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include "orsa_coord.h"
#include "orsa_secure_math.h"
#include "orsa_frame.h"

#include <string>
#include <vector>
#include <list>

namespace orsa {
  
  class Observation {
  public:
    // used to sort observations
    inline bool operator < (const Observation &obs) const {
      return (date.GetJulian() < obs.date.GetJulian());
    }
    
  public:
    std::string designation;
    std::string discovery_asterisk;
    Date   date;
    Angle  ra, dec;         
    double mag; // magnitude
    std::string obscode;
  };
  
  class Sky {
  public:
    void Set(Angle & ra, Angle & dec) {
      _ra    = ra;
      _dec   = dec;
    }
    
  public:
    inline Angle ra()  const { return _ra; }
    inline Angle dec() const { return _dec; }
    
  public:
    void Compute_J2000(const Vector & relative_position);
    
  public:
    double delta_arcsec(const Observation & obs) const;
    
  private:
    Angle _ra, _dec;
  };
  
  class Orbit {
  public:
    Orbit() { };
    
  public:
    inline double Period() const { return (secure_sqrt(4*pisq*a*a*a/mu)); }
    
    // eccentric anomaly
    double GetE() const;
    
  public:
    void RelativePosVel(Vector & relative_position, Vector & relative_velocity) const;
    
  public:
    void Compute(const Vector & dr, const Vector & dv, const double mu);
    void Compute(const Body & b, const Body & ref_b);
    
  public:
    double a,e,i,omega_node,omega_pericenter,M;
    double mu; // G*(m+M)
  };
  
  class OrbitWithEpoch : public Orbit {
  public:
    inline OrbitWithEpoch() : Orbit() { }
    inline OrbitWithEpoch(const Orbit & orbit) : Orbit(orbit) { }
    inline OrbitWithEpoch(const OrbitWithEpoch & orbit) : Orbit(orbit), epoch(orbit.epoch) { }
    
  public:
    /* 
       inline void Compute(const Vector & dr, const Vector & dv, const double mu) {
       Orbit::Compute(dr,dv,mu);
       }
    */
    
    inline void Compute(const Vector & dr, const Vector & dv, const double mu, const UniverseTypeAwareTime &epoch_in) {
      epoch = epoch_in;
      Orbit::Compute(dr,dv,mu);
    }
    
    /* 
       inline void Compute(const Body & b, const Body & ref_b) {
       Orbit::Compute(b, ref_b);
       }
    */
    
    inline void Compute(const Body & b, const Body & ref_b, const UniverseTypeAwareTime & epoch_in) {
      epoch = epoch_in;
      Orbit::Compute(b, ref_b);
    }
    
  public:    
    inline void RelativePosVel(Vector & relative_position, Vector & relative_velocity) const {
      Orbit::RelativePosVel(relative_position,relative_velocity);
    }
    
  public:
    inline void RelativePosVel(Vector & relative_position, Vector & relative_velocity, const UniverseTypeAwareTime &epoch_in) const {
      OrbitWithEpoch o(*this);
      o.M += twopi*(epoch_in.Time() - epoch.Time())/Period();
      o.M  = std::fmod(10*twopi+std::fmod(o.M,twopi),twopi);
      o.RelativePosVel(relative_position,relative_velocity);
    }
    
  public:
    UniverseTypeAwareTime epoch;
    double libration_angle; // to be removed soon!!
  };
  
  //! The RMS of the residuals in arcsec units
  double residual(const Observation&, const OrbitWithEpoch&);
  double RMS_residuals(const std::vector<Observation>&, const OrbitWithEpoch&);
  
  /* 
     Vector Position(const Orbit&);
     Vector Position(const Orbit&, const double);
     
     Vector Velocity(const Orbit&);
     Vector Velocity(const Orbit&, const double);
     
     void RelativePosVel(const Orbit&, Vector&, Vector&);
     void RelativePosVel(const Orbit&, Vector&, Vector&, const double);
     void RelativePosVel(const Orbit&, Vector&, Vector&, const Date);
  */
  
  // sky position of an object
  // Sky PropagatedSky(const OrbitWithEpoch&, const UniverseTypeAwareTime&, const std::string &obscode, bool integrate=true);
  // Sky PropagatedSky_J2000(const OrbitWithEpoch&, const UniverseTypeAwareTime&, const std::string &obscode, bool integrate=true);
  Sky PropagatedSky_J2000(const OrbitWithEpoch&, const UniverseTypeAwareTime&, const std::string &obscode, const bool integrate=true, const bool light_time_corrections=true);
  
  // an optimized class for frequent position computations of an object with a given orbit
  class OptimizedOrbitPositions {
  public:
    OptimizedOrbitPositions(const OrbitWithEpoch &orbit);
    
  public:
    // Sky PropagatedSky(const UniverseTypeAwareTime&, const std::string &obscode, const bool integrate=true, const bool light_time_corrections=true);
    Sky PropagatedSky_J2000(const UniverseTypeAwareTime&, const std::string &obscode, const bool integrate=true, const bool light_time_corrections=true);
    void PropagatedPosVel(const UniverseTypeAwareTime&, Vector&, Vector&, bool integrate=true);
    OrbitWithEpoch PropagatedOrbit(const UniverseTypeAwareTime&, bool integrate=true);
    
  private:
    const OrbitWithEpoch _orbit;
    std::list<JPL_planets> l;
    std::vector<Frame> frames;
  };
  
#ifndef HAVE_GSL // else: definitions in orsa_orbit_gsl.h
  class Asteroid {
  public:
    int n;    
    OrbitWithEpoch orb;
    std::string name;
    double mag;
  };
  
  class AsteroidDatabase : public std::vector<Asteroid> {
    
  };
  
  class ObservationCandidate { 
  public:
    Asteroid        a;
    Observation     o;
    Angle           delta;
  };
#endif // HAVE_GSL
  
  class Location {
  public:
    double lon,pxy,pz;
    std::string name;
    // bool present;  
  };
  
} // namespace orsa

#endif // _ORSA_ORBIT_H_
