/* 
   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_UNITS_H
#define _ORSA_UNITS_H

#include <cmath>
#include <string>
#include <cstdio>

#include "orsa_secure_math.h"
#include "orsa_coord.h"
#include "orsa_error.h"

namespace orsa {
  
  enum time_unit {
    YEAR=1,
    DAY=2,
    HOUR=3,
    MINUTE=4,
    SECOND=5
  };
  
  inline void convert(time_unit &tu, const unsigned int i)  {
    switch(i) {
    case 1: tu = YEAR;   break;
    case 2: tu = DAY;    break;
    case 3: tu = HOUR;   break;
    case 4: tu = MINUTE; break;
    case 5: tu = SECOND; break;
      //
    default:
      ORSA_ERROR("conversion problem: i = %i",i);
      break;       
    }
  }
  
  enum length_unit {
    MPARSEC=1,
    KPARSEC=2,
    PARSEC=3,
    LY=4,
    AU=5,
    EARTHMOON=6,
    REARTH=7,
    RMOON=8,
    KM=9,
    M=10,
    CM=11,
    // aliases
    LD=EARTHMOON,
    ER=REARTH,
    MR=RMOON
  };
  
  inline void convert(length_unit &lu, const unsigned int i)  {
    switch(i) {
    case 1:  lu = MPARSEC;   break;
    case 2:  lu = KPARSEC;   break;
    case 3:  lu = PARSEC;    break;
    case 4:  lu = LY;        break;
    case 5:  lu = AU;        break;
    case 6:  lu = EARTHMOON; break;
    case 7:  lu = REARTH;    break;
    case 8:  lu = RMOON;     break;
    case 9:  lu = KM;        break;
    case 10: lu = M;         break;
    case 11: lu = CM;        break;
      //
    default:
      ORSA_ERROR("conversion problem: i = %i",i);
      break;       
    }
  }
  
  enum mass_unit {
    MSUN=1,
    MJUPITER=2,
    MEARTH=3,
    MMOON=4,
    KG=5,
    GRAM=6
  };
  
  inline void convert(mass_unit &mu, const unsigned int i)  {
    switch(i) {
    case 1: mu = MSUN;     break;
    case 2: mu = MJUPITER; break;
    case 3: mu = MEARTH;   break;
    case 4: mu = MMOON;    break;
    case 5: mu = KG;       break;
    case 6: mu = GRAM;     break;    
      //
    default:
      ORSA_ERROR("conversion problem: i = %i",i);
      break;       
    }
  }
  
  template <class UNIT> class UnitBaseScale {
  public:
    UnitBaseScale() {}
    UnitBaseScale(UNIT unit) { base_unit = unit; }
 
  public:
    void Set(UNIT unit) { base_unit = unit; }
    const UNIT & GetBaseUnit() const { return base_unit; }
    
  private:
    UNIT base_unit;
  };
  
  class Units {
    
  private:
    UnitBaseScale<time_unit>   Time;
    UnitBaseScale<length_unit> Length;
    UnitBaseScale<mass_unit>   Mass;
    
  public:
    Units();
    Units(time_unit, length_unit, mass_unit);
    
  private:
    void init_base();
    
  private:
    void TryToSetUnitsFromJPLFile();
    
  public:
    void SetSystem(time_unit tu, length_unit lu, mass_unit mu);
    
    // constants
    inline double GetG()    const { return G; };
    inline double GetMSun() const { return MSun; };
    inline double GetC()    const { return c; }; // c = speed of light
    
    double GetG_MKS() const;
    
  public:
    std::string label(time_unit)   const;
    std::string label(length_unit) const;
    std::string label(mass_unit)   const;
    
  public:
    inline std::string TimeLabel()   const { return label(  Time.GetBaseUnit()); };
    inline std::string LengthLabel() const { return label(Length.GetBaseUnit()); };
    inline std::string MassLabel()   const { return label(  Mass.GetBaseUnit()); };
    
    // conversions
  public:
    /* 
       inline double FromUnits(double x, time_unit   t_in, double power=1.0) const { return (x*secure_pow(GetTimeScale(t_in)/GetTimeScale(),power));     }
       inline double FromUnits(double x, length_unit l_in, double power=1.0) const { return (x*secure_pow(GetLengthScale(l_in)/GetLengthScale(),power)); }
       inline double FromUnits(double x, mass_unit   m_in, double power=1.0) const { return (x*secure_pow(GetMassScale(m_in)/GetMassScale(),power));     }
    */
    
  private:
    inline static double __int_pow__(const double x, const int p) {
      if (p == 0) return 1.0;
      double _pow = x;
      const unsigned int max_k = static_cast<unsigned int>(std::abs(p));
      for (unsigned int k=1; k < max_k; ++k) {
	_pow *= x;
      }
      if (p < 0) _pow = 1.0/_pow;
      return _pow;
    }
    
    // conversions
  public:
    inline double FromUnits(const double x, const time_unit   t_in, const int power = 1) const { return (x*__int_pow__(GetTimeScale(t_in)/GetTimeScale(),power));     }
    inline double FromUnits(const double x, const length_unit l_in, const int power = 1) const { return (x*__int_pow__(GetLengthScale(l_in)/GetLengthScale(),power)); }
    inline double FromUnits(const double x, const mass_unit   m_in, const int power = 1) const { return (x*__int_pow__(GetMassScale(m_in)/GetMassScale(),power));     }
    
  public:
    inline const time_unit   & GetTimeBaseUnit()   const { return Time.GetBaseUnit(); };
    inline const length_unit & GetLengthBaseUnit() const { return Length.GetBaseUnit(); };
    inline const mass_unit   & GetMassBaseUnit()   const { return Mass.GetBaseUnit(); };
    
  protected:
    double GetTimeScale(  const   time_unit tu) const; 
    double GetLengthScale(const length_unit lu) const; 
    double GetMassScale(  const   mass_unit mu) const; 
    
    double GetTimeScale()   const;
    double GetLengthScale() const; 
    double GetMassScale()   const; 
    
    void Recompute();
    
  private:
    double G,G_base;
    double MSun,MSun_base;
    double MJupiter_base, MEarth_base, MMoon_base;
    double AU_base;
    double c,c_base;
    double r_earth_base;
    double r_moon_base;
    double parsec_base;
  };
  
  // world visible units // defined in orsa_universe.cc
  extern Units * units;
  
  // defined in orsa_universe.cc
  double GetG();
  double GetG_MKS();
  double GetMSun();
  double GetC();
  
  // everything in orsa_universe.cc
  /* 
     double FromUnits(double, time_unit,   double = 1.0);
     double FromUnits(double, length_unit, double = 1.0);
     double FromUnits(double, mass_unit,   double = 1.0); 
  */
  //
  double FromUnits(const double, const time_unit,   const int = 1);
  double FromUnits(const double, const length_unit, const int = 1);
  double FromUnits(const double, const mass_unit,   const int = 1); 
  //
  std::string TimeLabel();
  std::string LengthLabel();
  std::string MassLabel();
  
  //! TimeScale enum, useful only when using a Real Universe. More information can be obtained here: http://www.hartrao.ac.za/nccsdoc/slalib/sun67.htx/node217.html
  enum TimeScale {
    UTC=1,
    UT=2,
    TAI=3,
    TDT=4,
    GPS=5,
    // aliases
    UT1=UT,
    ET=TDT,
    TT=TDT
  };
  
  inline void convert(TimeScale &ts, const unsigned int i)  {
    switch(i) {
    case 1: ts = UTC; break;
    case 2: ts = UT;  break;
    case 3: ts = TAI; break;
    case 4: ts = TDT; break;
    case 5: ts = GPS; break;
      //
    default:
      ORSA_ERROR("conversion problem: i = %i",i);
      break;       
    }
  }
  
  extern TimeScale default_Date_timescale;
  
  std::string TimeScaleLabel(TimeScale);
  
  class UniverseTypeAwareTime;
  
  class UniverseTypeAwareTimeStep;
  
  class TimeStep {
  public:
    TimeStep();
    TimeStep(const unsigned int days, const unsigned int day_fraction, const int sign);
    TimeStep(const double t);
    TimeStep(const TimeStep &);
    
  public:
    double GetDouble() const;
    
  public:
    TimeStep & operator += (const TimeStep &);
    TimeStep & operator -= (const TimeStep &);
    
    TimeStep operator + (const TimeStep &) const;
    TimeStep operator - (const TimeStep &) const;
    
    TimeStep & operator *= (const int);
    TimeStep & operator *= (const double);
    
  public:
    void AddDays(const unsigned int, const int);
    void AddDayFractions(const unsigned int, const int);
    
    // sign operators 
    TimeStep operator + () const;
    TimeStep operator - () const;
    
    bool operator == (const TimeStep &) const;
    bool operator != (const TimeStep &) const;
    
    bool operator > (const TimeStep &) const;
    bool operator < (const TimeStep &) const;
    
    inline bool operator >= (const TimeStep & ts) const {
      if ((*this) == ts) return true;
      if ((*this) >  ts) return true;
      return false;
    }
    
    inline bool operator <= (const TimeStep & ts) const {
      if ((*this) == ts) return true;
      if ((*this) <  ts) return true;
      return false;      
    }
    
    TimeStep absolute() const;
    
    bool IsZero() const;
    
    static unsigned int max_day_fraction() {
      return 864000000; // 86400 ( = second per day) * 10000 ( = 1.0e-4 seconds resolution)
    }
    
  private:
    void internal_check();
    
  public:
    unsigned int days()         const { return _days;         }	
    unsigned int day_fraction() const { return _day_fraction; }
    int          sign()         const { return _sign;         }
    
  private:
    unsigned int _days;         // always positive
    unsigned int _day_fraction; // always positive
    int          _sign;         // +1 or -1
  };
  
  /*! A note on Gregorian and Julian calendars:
   * since its introduction in 1582, the Gregorian calendar
   * has been adopted in different years from different countries,
   * so the date obtained from the Date class can be different from 
   * the real 'old' one which used the Julian date, and usually
   * the difference is of a few days.
   *
   * In particular, ORSA applies the Gregorian calendar in any epoch,
   * i.e. even before 1582. This appears to be the simplest solution,
   * at least for the moment.
   * 
   * For more info, check out i.e. http://www.dome-igm.com/convers.htm
   */
  
  class Date {
  public:
    Date();
    Date(const Date &);
    Date(const UniverseTypeAwareTime &);
    
  public:
    void SetGregor(int   y, int   m, double  d,                              TimeScale ts = default_Date_timescale);
    void SetGregor(int   y, int   m, int     d, int H, int M, int S, int ms, TimeScale ts = default_Date_timescale);
    //
    void GetGregor(int & y, int & m, int &   d, TimeScale ts = default_Date_timescale) const;
    double GetDayFraction(                      TimeScale ts = default_Date_timescale) const;
    unsigned int GetDayFraction_unsigned_int(   TimeScale ts = default_Date_timescale) const;
    //
    void   SetJulian(double,               TimeScale ts = default_Date_timescale);
    void   GetJulian(double &,             TimeScale ts = default_Date_timescale) const;
    double GetJulian(                      TimeScale ts = default_Date_timescale) const;
    // 
    double Time()    const;
    double GetTime() const;
    void   SetTime(const double);
    
    Date & operator += (const UniverseTypeAwareTimeStep &);
    Date & operator -= (const UniverseTypeAwareTimeStep &);
    
    bool operator == (const Date &) const;
    bool operator != (const Date &) const;
    
    bool operator < (const Date &) const;
    bool operator > (const Date &) const;
    
  public:
    void SetNow();
    void SetToday();
    void SetJ2000();
    
  public:
    //! delta_seconds = to - from   ('to' minus 'from')
    static double delta_seconds(int y, int m, int d, const TimeScale from, const TimeScale to=default_Date_timescale);
    
    //! returns a TimeStep as big as the Date
    TimeStep GetTimeStep() const {
      return TimeStep(sdn,df,+1);
    }
    
  private:
    unsigned int sdn;  // internal day format
    unsigned int  df;  // day fraction; compare it with TimeStep::max_day_fraction()
  };
  
  /* 
     inline bool operator == (const Date &x, const Date &y) {
     return (x.GetTimeStep() == y.GetTimeStep());
     }
     
     inline bool operator!= (const Date &x, const Date &y) {
     return !(x==y);
     }
  */
  
  // UniverseTypeAwareTime
  
  //! A Date like class which represents time
  //! as Date in Real Universes and as a double
  //! in Simulated Universes  
  class UniverseTypeAwareTime {
  public:
    UniverseTypeAwareTime();
    UniverseTypeAwareTime(const double);
    UniverseTypeAwareTime(const Date&);
    UniverseTypeAwareTime(const UniverseTypeAwareTime&);
    
  public:
    double Time()    const;
    double GetTime() const;
    Date   GetDate() const;
    virtual void SetTime(const double);
    virtual void SetTime(const Date &);
    virtual void SetDate(const Date &);
    virtual void SetTime(const UniverseTypeAwareTime &);
    
    UniverseTypeAwareTime & operator += (const UniverseTypeAwareTimeStep &);
    UniverseTypeAwareTime & operator -= (const UniverseTypeAwareTimeStep &);
    
    // ***
    // UniverseTypeAwareTimeStep operator + (const UniverseTypeAwareTime &) const;
    UniverseTypeAwareTimeStep operator - (const UniverseTypeAwareTime &) const;
    
    // UniverseTypeAwareTimeStep operator + (const UniverseTypeAwareTimeStep &) const;
    // UniverseTypeAwareTimeStep operator - (const UniverseTypeAwareTimeStep &) const;
    
    UniverseTypeAwareTime operator + (const UniverseTypeAwareTimeStep &) const;
    UniverseTypeAwareTime operator - (const UniverseTypeAwareTimeStep &) const;
    
    bool operator == (const UniverseTypeAwareTime &) const;
    bool operator != (const UniverseTypeAwareTime &) const;
    
    bool operator < (const UniverseTypeAwareTime &) const;
    bool operator > (const UniverseTypeAwareTime &) const;
    
    bool operator <= (const UniverseTypeAwareTime &) const;
    bool operator >= (const UniverseTypeAwareTime &) const;
    
  protected:
    double time;
    Date   date;
  };
  
  class UniverseTypeAwareTimeStep {
  public:
    UniverseTypeAwareTimeStep();
    UniverseTypeAwareTimeStep(const double);
    UniverseTypeAwareTimeStep(const TimeStep &);
    UniverseTypeAwareTimeStep(const int days, const unsigned int day_fraction, const int sign = +1);
    UniverseTypeAwareTimeStep(const UniverseTypeAwareTimeStep &);
    
  public:
    UniverseTypeAwareTimeStep & operator += (const UniverseTypeAwareTimeStep &);
    UniverseTypeAwareTimeStep & operator -= (const UniverseTypeAwareTimeStep &);
    UniverseTypeAwareTimeStep & operator *= (const int);
    UniverseTypeAwareTimeStep & operator *= (const double);
    
  public:
    double GetDouble() const;
    
  public:
    UniverseTypeAwareTimeStep operator + () const;
    UniverseTypeAwareTimeStep operator - () const;
  
  public:
    UniverseTypeAwareTimeStep operator + (const UniverseTypeAwareTimeStep &) const;
    UniverseTypeAwareTimeStep operator - (const UniverseTypeAwareTimeStep &) const;
    
    UniverseTypeAwareTimeStep operator + (const UniverseTypeAwareTime &) const;
    UniverseTypeAwareTimeStep operator - (const UniverseTypeAwareTime &) const;
    
    UniverseTypeAwareTimeStep absolute() const;
    
    bool IsZero() const;
    
    bool operator < (const UniverseTypeAwareTimeStep &) const;
    bool operator > (const UniverseTypeAwareTimeStep &) const;
    
    bool operator < (const UniverseTypeAwareTime &) const;
    bool operator > (const UniverseTypeAwareTime &) const;
    
    bool operator < (const double) const;
    bool operator > (const double) const;
    
    bool operator == (const UniverseTypeAwareTimeStep &) const;
    bool operator != (const UniverseTypeAwareTimeStep &) const;
    
  public:
    unsigned int days()         const { return ts.days();         }	
    unsigned int day_fraction() const { return ts.day_fraction(); }
    int          sign()         const { return ts.sign();         }
    
    TimeStep GetTimeStep() const { return ts; };
    
    void     SetTimeStep(const TimeStep &ts_in) { ts = ts_in; }
    
    void SetDouble(const double d) { dts = d; }
    
  protected:
    TimeStep  ts;    
    double   dts;
  };
  
  UniverseTypeAwareTimeStep operator * (const int, const UniverseTypeAwareTimeStep &);
  UniverseTypeAwareTimeStep operator * (const UniverseTypeAwareTimeStep &, const int);
  
  UniverseTypeAwareTimeStep operator * (const double, const UniverseTypeAwareTimeStep &);
  UniverseTypeAwareTimeStep operator * (const UniverseTypeAwareTimeStep &, const double);
  
  // Angle
  
  class Angle {
  public:
    Angle() { radians = 0; }
    Angle(double x) { radians = x; }
    
  public:
    void   SetRad(double);
    void   GetRad(double&) const;
    double GetRad() const;
    //
    void   SetDPS(double  d, double  p, double  s);
    void   GetDPS(double &d, double &p, double &s) const;
    //
    void   SetHMS(double  h, double  m, double  s);
    void   GetHMS(double &h, double &m, double &s) const;

  private:
    double radians;
  };
  
  inline double sin(const Angle & alpha) {
    return std::sin(alpha.GetRad());
  }
  
  inline double cos(const Angle & alpha) {
    return std::cos(alpha.GetRad());
  }
  
  inline double tan(const Angle & alpha) {
    return std::tan(alpha.GetRad());
  }
  
  inline void sincos(const Angle & alpha, double & s, double & c) {
#ifdef HAVE_SINCOS
    ::sincos(alpha.GetRad(),&s,&c); 
#else // HAVE_SINCOS
    s = std::sin(alpha.GetRad());
    c = std::cos(alpha.GetRad());
#endif // ! HAVE_SINCOS
  }
  
  // reference systems
  
  enum ReferenceSystem {
    EQUATORIAL=1,
    ECLIPTIC=2
  };
  
  inline void convert(ReferenceSystem &rs, const unsigned int i)  {
    switch(i) {
    case 1: rs = EQUATORIAL; break;
    case 2: rs = ECLIPTIC; break;
      //
    default:
      ORSA_ERROR("conversion problem: i = %i",i);
      break;
    }
  }
  
  Angle obleq(const Date &);
  // Angle delta_obleq(const Date&);
  // Angle delta_longitude(const Date&);
  Angle gmst(const Date &);
  
  // rotations  
  void EclipticToEquatorial(Vector&, const Date&);
  void EquatorialToEcliptic(Vector&, const Date&);
  
  // Angle gmst(const Date&);
  // void EquatorialToEcliptic(Vector&, const Date&);
  //
  // void EquatorialToDate_To_EquatorialJ2000(Vector&, const Date&);
  
  // J2000 versions
  Angle obleq_J2000();
  // Angle delta_obleq_J2000();
  // Angle delta_longitude_J2000();
  //
  void EclipticToEquatorial_J2000(Vector&);
  void EquatorialToEcliptic_J2000(Vector&);
  
} // namespace orsa

#include "orsa_file_jpl.h"

namespace orsa {
  
  // another, improved orientation/rotation
  // formula, from the "Report of the IAU/IAG working group on
  // cartographic coordinates and rotational elements of the planets
  // and satellites: 2000", P.K. Seidelmann et al,
  // Celestial Mechanics and Dynamical Astronomy 82: 83-110 (2002).
  // data from Table I, pag. 86.
  void alpha_delta_meridian(const JPL_planets, const Date &,
			    Angle & alpha_zero, Angle & delta_zero, Angle & W);
  
} // namespace orsa

#endif // _ORSA_UNITS_H
