#ifndef __JASTRONOMY__
#define __JASTRONOMY__

#include <iostream>

#include "JMath/JConstants.hh"
#include "JMath/JMatrix3D.hh"
#include "JMath/JSVD3D.hh"
#include "JGeometry3D/JAngle3D.hh"

/**
 * \file
 * 
 * Interface methods for slalib and auxiliary classes and methods for astronomy.
 * \author mdejong, fhuang, azegarelli
 */
namespace JASTRONOMY {}
namespace JPP { using namespace JASTRONOMY; }

namespace JASTRONOMY {

  extern "C"
  {
    // To convert Modified Julian Date to Greenwich Mean Sidereal Time (in rads):
    double sla_gmst_(double& ut1);
    double sla_gmsta_(double& ut1, double& part_day);

    // To calculate the difference between the mean & apparant places of the Equinox (in rads):
    double sla_eqeqx_(double& ut1);

    // To convert local equatorial to horizontal coords and vice versa:
    void sla_de2h_(double& hourangle, double& declination, double& observer_latitude,
                   double& azimuth, double& elevation); // azimuth: [0,2pi], elevation: [-pi/2,pi/2]
    void sla_dh2e_(double& azimuth, double& elevation, double& observer_latitude,
                   double& hourangle, double& declination); // HA: [-pi,pi], Dec: [-pi/2,pi/2]

    // To convert equatorial to ecliptic coordinates and vice versa:
    void sla_eqecl_(double& rightascension, double& declination, double& time,
                    double& ecliptic_longitude, double& ecliptic_latitude);
    void sla_ecleq_(double& ecliptic_longitude, double& ecliptic_latitude, double& time,
                    double& rightascension, double& declination);

    // To convert equatorial to galactic coords and vice versa:
    void sla_eqgal_(double& rightascension, double& declination,
                    double& galactic_longitude, double& galactic_latitude);
    void sla_galeq_(double& galactic_longitude, double& galactic_latitude,
                    double& rightascension, double& declination);

    // To convert galactic to supergalactic coords and vice versa:
    void sla_galsup_(double& galactic_longitude, double& galactic_latitude,
                     double& sgalactic_longitude, double& sgalactic_latitude);
    void sla_supgal_(double& sgalactic_longitude, double& sgalactic_latitude,
                     double& galactic_longitude, double& galactic_latitude);

    // To convert Gregorian Calender Date to Modified Julian Date:
    // (Note: MJD = JD - 2400000.5,
    //  i.e.\ JD = # of days since     noon  1st of  januari 4713 BC
    //  and MJD = # of days since midnight 17th of november 1858 AD)
    void sla_caldj_(int& year, int& month, int& day,
                    double& mjd, int& status);

    // calendar to year, days
    void sla_clyd_(int& year, int& month, int& day, int& nyears, int& ndays, int&status);

    // MJD to Gregorian Calendar Date, expressed in single array
    void sla_djcal_(int& precision, double& mjd,  // input : number of decimal places of days in fraction, MJD
                    int result[4],  int& status); // output: GCD, status

    // days to hour, min, sec
    void sla_dd2tf_(int& ndec, double& day,       // input : number of decimal places of seconds, interval in days
                    char sign[1], int result[4]); // output: + or -, hours, minutes, seconds

    // To obtain the difference in dynamical Terrestial Time and the current time (UT1) in secs, at the current time (UT1):
    // (Note: instead of the current time (UT1) one can also use current time (UTC) since UTC-UT1 < 0.9 secs
    double sla_dtt_(double& ut1);

    // To obtain approx. topocentric apparent (RA & Dec) and angular size of planet/moon/sun from an observer:
    void sla_rdplan_(double& dtt, int& object, const double& longitude, const double& latitude, // longitude is +/- if east/west of prime meridian
                     double& rightascension, double& declination, double& diam);

    // To calculate the precession between epochs:
    void sla_preces_(char *system, double& ep0, double& ep1, double& ra, double& dec, int length);

    // To calculate the matrix of precession and nutation (IAU 1976, FK)
    // input: epoch, Julian Epoch for mean coordinates; mjd, modified Julian date (JD -2400000.5) for true coordinates
    // output: combined precission-nutation matrix
    void sla_prenut_(double& epoch, double& mjd, double rmatpn[3][3] );

  }


  static const double MJD_EPOCH = 40587.0;               // MJD of the Unix Epoch, i.e.\ 00:00 Jan 1 1970 (UTC)

  static const double NUMBER_OF_SECONDS_PER_HOUR         = 60.0 * 60.0;
  static const double NUMBER_OF_SECONDS_PER_DAY          = NUMBER_OF_SECONDS_PER_HOUR *  24.0; 
  static const double NUMBER_OF_SECONDS_PER_YEAR         = NUMBER_OF_SECONDS_PER_DAY  * 365.0;
  static const double NUMBER_OF_SECONDS_PER_SEDERIAL_DAY = NUMBER_OF_SECONDS_PER_HOUR *  23.9344696;

  /*
  * Bring angle between 0 and 2 pi
  * \param angle          Input angle [rad]
  * \return               angle [rad]  
  */ 
  inline double wrap( double angle )
  {
    int n = int(angle / (2 * JMATH::PI)) - (angle < 0);
    angle -= n * 2 * JMATH::PI;
    return angle;
  }

  /*
  * Get UTM zone from longitude
  * \param angle          Input longitude [rad] 
  * \return               UTM zone 
  */
  inline int get_utm_zone(double lat)
  {
    return 1 + int( ( JMATH::PI + lat ) / ( 6 * JMATH::PI / 180 ) );
  }

  /* Compute longitude of the central meridian given the UTM zone
   * \param utmzone        Input UTM zone [int]
   * \return               Longitude of central meridian [rad] 
   */
  inline double longitude_of_central_meridian(int utmzone)
  {
    const double zone_width = 6 * JMATH::PI / 180;
    return  -JMATH::PI + (utmzone -1)*zone_width + zone_width/2;
  }

  /*
   * Compute meridian convergence angle
   * Angular difference of the UTM Northing direction to the geodetic North. It is positive to the East.
   * For example, for ANTARES, convergence angle = -1.93 deg E, i.e.UTM North pointing a bit towards the geographical west.
   * Code from aanet/evt/Det.cc
   * \param longitude       longitude [rad]
   * \param latitude        latitude  [rad]
   * \return                meridian convergence angle [rad]      
   */
  inline double compute_meridian_convergence_angle( double longitude, double latitude )
  {

    double latitude_deg = latitude * 180/JMATH::PI;

    if ( latitude_deg > 84. || latitude_deg < -80.)
      {
	std::cout << "UTM coordinate system not defined for given Latitude: %f (max 84 deg N or min -80 deg S), setting meridian convergence angle to 0 deg.";
	return 0;
      }

    // detector position, longitude and latitude in rad
    double lambda  = longitude;
    double phi     = latitude;

    // find UTM zone and central meridian
    // longitude of the central meridian of UTM zone in rad
    double lambda0 = longitude_of_central_meridian( get_utm_zone( longitude ) );
    double omega   = lambda-lambda0;

    // parameters of the Earth ellipsoid
    double a  = 6378137.;      // semi-major axis in meters (WGS84)
    double e2 = 0.0066943800;  // eccentricity (WGS84)

    double rho = a*(1.-e2)/pow(1.-e2*pow(sin(phi),2),3./2.);
    double nu  = a/sqrt(1-e2*pow(sin(phi),2.));
    double psi = nu/rho;
    double t = tan(phi);

    // power series for convergence angle 
    double gamma = sin(phi) * omega
      - sin(phi) * pow(omega,3.)/3.  * pow(cos(phi),2.) * (2.*pow(psi,2.)-psi)
      - sin(phi) * pow(omega,5.)/15. * pow(cos(phi),4.) * (pow(psi,4.)*(11.-24.*pow(t,2.))-
							   pow(psi,3.)*(11.-36.*pow(t,2.))+
							   2.*pow(psi,2.)*(1.-7.*pow(t,2.))+
							   psi*pow(t,2.)
							   )
      - sin(phi) * pow(omega,7.)/315. * pow(cos(phi),6.) * (17.-26.*pow(t,2.)+2.*pow(t,4.));

    return gamma;
  }

  /** 
   * Convert (Dec, RA) to J2000. Adapted from aanet/astro/Astro.cc
   *
   * \param mjd          Modified Julian Date (= # of days since midnight 17th of november 1858 AD)
   * \param in_dec       Input dec
   * \param in_ra        Input RA
   * \param out_dec      Output dec
   * \param out_ra       Output RA
   * \param  reverse     if reverse = false, this function converts from ra,dec at a certain time to J2000, if reverse = true, it goes the other way (i.e. if you have a catalog position and want to compute a track, use reverse = true) 
  */

  inline void correct_to_j2000(double& in_dec, double& in_ra, double& mjd, double& out_dec, double& out_ra, bool reverse=false)
  {
    // --- get the rotation matrix ----
    // NOTE: (from seatray Astro): Julian epoch of J2000 time definition
    // is 2000, corrsponding by definition to 2451545.0 TT (unmodified)
    // Julian date!

    double epoch = 2000;
    double rmatpn[3][3] {};
    sla_prenut_(epoch, mjd, rmatpn);         // SLALIB function to get the matrix of precession and nutation
 
    JMATH::JMatrix3D M = JMATH::JMatrix3D(rmatpn[0][0], rmatpn[0][1], rmatpn[0][2], rmatpn[1][0], rmatpn[1][1], rmatpn[1][2], rmatpn[2][0], rmatpn[2][1], rmatpn[2][2]);
    JMATH::JMatrix3D& M_T = M.transpose();    // transpose 'cause of fortan matrix convention 
    JMATH::JSVD3D V(M_T);

    double theta = JMATH::PI/2 - in_dec;
    double phi = in_ra;             // get the (theta, phi) from (in_dec, in_ra)
    double v[3];                    // get the vector in the form (x, y, z) from (theta, phi)
    v[0] = sin (theta) * cos(phi);
    v[1] = sin (theta) * sin(phi);
    v[2] = cos (theta);
    // multiplying the matrix with v and get the transformed vector 
    if (reverse){
      M_T.transform(v[0],v[1],v[2]);
    }
    else{
      JMATH::JMatrix3D M_T_Inv;
      M_T_Inv = V.invert();
      M_T_Inv.transform(v[0],v[1],v[2]);
    }

    // get the (phi, theta) from v2
    double phi_2 = atan2( v[1], v[0] ); double theta_2 =  acos(v[2]);
    // convert (phi,theta) to declination and right ascention
    out_ra = wrap(phi_2);
    out_dec = JMATH::PI/2 - theta_2;
  }


  /**
   * Convert angle to radians.
   *
   * \param  angle         angle [deg]
   * \return               angle [rad]
   */
  inline double getRadians(const double angle)
  {
    return JMATH::PI * angle / 180.0;
  }

  
  /**
   * Convert angle to radians.
   *
   * \param  angle         angle [deg]
   * \param  amin          arcminutes
   * \param  asec          arcseconds
   * \return               angle [rad]
   */
  inline double getRadians(const int    angle,
			   const int    amin,
			   const double asec)
  {
    return JMATH::PI * ((double) angle          +
			(double) amin  /   60.0 + 
			(double) asec  / 3600.0)  /  180.0;
  }

  
  /**
   * Convert hour angle to radians.
   *
   * \param  hour          hour
   * \param  min           minutes
   * \param  sec           seconds
   * \return               angle [rad]
   */
  inline double getHourAngle(const int    hour,
			     const int    min,
			     const double sec)
  {
    double ha = (JMATH::PI * (double) hour /    12.0 + 
		 JMATH::PI * (double) min  /   720.0 + 
		 JMATH::PI * (double) sec  / 43200.0);

    if (ha > JMATH::PI) {
      ha -= 2*JMATH::PI;
    }

    return ha;
  }


  /**
   * Location of astrophysical source.
   */
  class JSourceLocation {
  public:

    /**
     * Default constructor.
     */
    JSourceLocation() :
      __declination    (0.0),
      __right_ascension(0.0)
    {}


    /**
     * Constructor.
     *
     * \param    declination      declination 
     * \param    right_ascension  right ascension
     */
    JSourceLocation(const double declination,
		    const double right_ascension) :
      __declination    (declination),
      __right_ascension(right_ascension)
    {}

    const double& getDeclination()    const { return __declination; }
    const double& getRightAscension() const { return __right_ascension; }

    /**
     * Get declination in J2000.
     *
     * \param   t1         number of seconds since MJD
     */
    double getDeclinationJ2000(const double t1) const
    {
      double mjd = t1 / NUMBER_OF_SECONDS_PER_DAY;
      double dec_j2000; double ra_j2000;
      double &dec = const_cast<double&> (__declination);
      double &ra = const_cast<double&> (__right_ascension);
      correct_to_j2000(dec, ra, mjd, dec_j2000, ra_j2000, false);
      const double const_dec_j2000 = dec_j2000;
      return const_dec_j2000;
    }

    /**
     * Get Right Ascension in J2000.
     *
     * \param   t1         number of seconds since MJD
     */
    double getRightAscensionJ2000(const double t1) const
    {
      double mjd = t1 / NUMBER_OF_SECONDS_PER_DAY;
      double dec_j2000; double ra_j2000;
      double &dec = const_cast<double&> (__declination);
      double &ra = const_cast<double&> (__right_ascension);
      correct_to_j2000(dec, ra, mjd, dec_j2000, ra_j2000, false);
      const double const_ra_j2000 = ra_j2000;
      return const_ra_j2000;
    }

  protected:
    double __declination;
    double __right_ascension;
  };

  /**
   * Location of astrophysical source in Galactic coordinates.
   */
  class JGalacticCoordinates {
  public:

    /**
     * Default constructor.
     */
    JGalacticCoordinates() :
      __gal_latitude (0.0),
      __gal_longitude(0.0)
    {}

    /**
     * Constructor.
     *
     * \param    gal_latitude       Galactic latitude  [rad]
     * \param    gal_longitude      Galactic longitude [rad]
     */
    JGalacticCoordinates(const double gal_latitude,
                         const double gal_longitude) :
      __gal_latitude (gal_latitude),
      __gal_longitude(gal_longitude)
    {}

    const double& getGalacticLatitude()  const { return __gal_latitude; }
    const double& getGalacticLongitude() const { return __gal_longitude; }

  protected:
    double __gal_latitude;
    double __gal_longitude;
  };

  /**
   * Direction of incident neutrino.
   */
  class JNeutrinoDirection {
  public:

    /**
     * Default constructor.
     */
    JNeutrinoDirection() :
      __zenith (0.0),
      __azimuth(0.0)
    {}


    /**
     * Constructor.
     *
     * \param    zenith      zenith 
     * \param    azimuth     azimuth
     */
    JNeutrinoDirection(const double& zenith,
		       const double& azimuth) :
      __zenith (zenith),
      __azimuth(azimuth)
    {}

    const double& getZenith()  const { return __zenith; }
    const double& getAzimuth() const { return __azimuth; }

  protected:
    double __zenith;
    double __azimuth;
  };


  /**
   * Location of detector.
   */
  class JGeographicalLocation {
  public:

    /**
     * Default constructor.
     */
    JGeographicalLocation() :
      __latitude (0.0),
      __longitude(0.0)
    {}


    /**
     * Constructor.
     *
     * \param    latitude       latitude 
     * \param    longitude      longitude
     */
    JGeographicalLocation(const double latitude,
			  const double longitude) :
      __latitude (latitude),
      __longitude(longitude)
    {}


    /**
     * Constructor.
     *
     * \param    degreesNorth   degrees North
     * \param    minutesNorth   minutes North
     * \param    degreesEast    degrees East
     * \param    minutesEast    minutes East
     */
    JGeographicalLocation(const int degreesNorth,
			  const int minutesNorth,
			  const int degreesEast,
			  const int minutesEast) 
    {
      __latitude  = (degreesNorth + minutesNorth / 60.0) * JMATH::PI / 180.0;
      __longitude = (degreesEast  + minutesEast  / 60.0) * JMATH::PI / 180.0;  
    }

    const double& getLatitude()  const { return __latitude; }
    const double& getLongitude() const { return __longitude; }

  protected:
    double __latitude;
    double __longitude;
  };



  /**
   * Auxiliary class to make coordinate transformations for a specific geographical location of the detector.
   *
   * Note: SLALIB  ref. system                 : (x,y,z) = (N,E,up)
   *       ANTARES ref. system for d10_c00_s00 : (x,y,z) = (N,W,up)
   *       ANTARES ref. system for real det.   : (x,y,z) = (E,N,up)
   */
  class JAstronomy : 
    public JGeographicalLocation 
  {
  public:

    /**
     * Constructor.
     *
     * \param    location       location of detector
     */
    JAstronomy(const JGeographicalLocation& location) :
      JGeographicalLocation(location)
    {}


    /**
     * Get direction pointing to source.
     *
     * \param   t1         number of seconds since MJD
     * \param   pos        source location
     * \return             direction of neutrino
     */
    JNeutrinoDirection getDirectionOfNeutrino(const double& t1,
					      const JSourceLocation& pos) const
    {
      double mjd       = t1 / NUMBER_OF_SECONDS_PER_DAY;     // [days]

      double longitude = getLongitude();
      double latitude  = getLatitude();

      // Convert current mjd (UTC) to local sidereal time, 
      // taking into account the Equation of the Equinoxes:
      // Note: LST = GMST + local longitude, where l.l. is +/- if east/west of prime meridian

      double gmst      = sla_gmst_ (mjd);                    // gmst = Greenwich mean sidereal time
      double eqeqx     = sla_eqeqx_(mjd);                    // Note: Equation of the Equinoxes < 1.15 secs
      double lst       = gmst + longitude + eqeqx;
      
      // Transform time-independent equatorial coordinates to time-dependent equatorial coordinates (i.e.\ ra->ha):

      double dec = pos.getDeclination();
      double ra  = pos.getRightAscension();
      double ha  = lst - ra;
      
      // Convert time-dependent equatorial coordinates to local horizontal coordinates:
      // Note: azimuth: [0,2pi] and elevation: [-pi,pi]

      double azimuth;
      double elevation;
      
      sla_de2h_(ha, dec, latitude, azimuth, elevation);

      double theta = -elevation + JMATH::PI/2.0;
      double phi   = -azimuth   + JMATH::PI/2.0;

      // invert direction

      theta = JMATH::PI  - theta;
      phi   = phi + JMATH::PI;

      if (phi > JMATH::PI)
	phi -= 2.0*JMATH::PI;

      return JNeutrinoDirection(theta,phi);
    }

    /**
     *  Get location of source given a neutrino direction (zenith,azimuth) and time. 
     *
     * \param   t1         number of seconds since MJD 
     * \param   dir        direction of neutrino                                                                                                                      
     * \return             source location  
    */

    JSourceLocation getLocationOfSourceFromZenithAzimuth(const double t1,
							 const JGEOMETRY3D::JAngle3D& dir) const
    {
      double longitude = getLongitude();
      double latitude  = getLatitude();

      double zenith  = JMATH::PI  - dir.getTheta();
      double azimuth_utm = wrap(dir.getPhi() - JMATH::PI);

      double meridian_convergence_angle = compute_meridian_convergence_angle(longitude, latitude);
      double azimuth_geo = azimuth_utm - meridian_convergence_angle;

      double elevation   = -zenith  + JMATH::PI/2.0;
      double azimuth_sla = wrap(JMATH::PI/2.0 - azimuth_geo);

      double ha;
      double dec;

      sla_dh2e_(azimuth_sla, elevation, latitude, ha, dec);

      double mjd = t1 / NUMBER_OF_SECONDS_PER_DAY;       // [days]
      double gmst  = sla_gmst_ (mjd);                    // gmst = Greenwich mean sidereal time
      double eqeqx = sla_eqeqx_(mjd);                    // Note: Equation of the Equinoxes < 1.15 secs
      double lst   = gmst + longitude + eqeqx;
      double ra = lst - ha;

      ra = wrap(ra);  // Bring RA in [0, 2pi]
      return JSourceLocation(dec, ra);

    }

    /** 
     * Get location of source in galactic coordinates given a neutrino direction and time.
     *
     * \param   t1         number of seconds since MJD
     * \param   dir        direction of neutrino                                                                            
     * \return             source galactic coordinates
     */

    JGalacticCoordinates getGalacticCoordinatesOfSource(const double t1,
							const JGEOMETRY3D::JAngle3D& dir) const
    {
      const JSourceLocation& source_location = getLocationOfSourceFromZenithAzimuth(t1, dir);
      const double dec_j2000 = source_location.getDeclinationJ2000(t1);
      const double ra_j2000 =  source_location.getRightAscensionJ2000(t1);
      double &dec_j2000_n = const_cast<double&> (dec_j2000);
      double &ra_j2000_n  = const_cast<double&> (ra_j2000);
      double galactic_longitude; double galactic_latitude;
      sla_eqgal_(ra_j2000_n, dec_j2000_n, galactic_longitude, galactic_latitude);
      return JGalacticCoordinates(galactic_latitude, galactic_longitude);
    }

    /**
     * Get location of source.
     *
     * \param   t1         number of seconds since MJD
     * \param   dir        direction of neutrino
     * \return             source location
     */
    JSourceLocation getLocationOfSource(const double t1,
					const JNeutrinoDirection& dir) const
    {
      double longitude = getLongitude();
      double latitude  = getLatitude();
      
      // invert direction

      double theta = JMATH::PI  - dir.getZenith();
      double phi   = dir.getAzimuth() - JMATH::PI;

      if (phi < 0.0) {
	phi += 2.0*JMATH::PI;
      }

      double elevation = -theta  + JMATH::PI/2.0;
      double azimuth   = -phi    + JMATH::PI/2.0;
      
      // Convert time-dependent equatorial coordinates to local horizontal coordinates:
      // Note: azimuth: [0,2pi] and elevation: [-pi,pi]

      double ha;
      double dec;
      
      sla_dh2e_(azimuth, elevation, latitude, ha, dec);
      
      // Get current mjd (UTC):

      double mjd = t1 / NUMBER_OF_SECONDS_PER_DAY;       // [days]
      
      // Convert current mjd (UTC) to local sidereal time, 
      // taking into account the Equation of the Equinoxes:
      // Note: LST = GMST + local longitude, where l.l. is +/- if east/west of prime meridian

      double gmst  = sla_gmst_ (mjd);                    // gmst = Greenwich mean sidereal time
      double eqeqx = sla_eqeqx_(mjd);                    // Note: Equation of the Equinoxes < 1.15 secs
      double lst   = gmst + longitude + eqeqx;
      
      // Transform time-independent equatorial coordinates to time-dependent equatorial coordinates (i.e.\ ra->ha):

      double ra = lst - ha;

      if (ra > JMATH::PI) {
	ra -= 2.0*JMATH::PI;
      }

      return JSourceLocation(dec, ra);
    }
  };


  /**
   * Dot product.
   *
   * \param  first     neutrino direction
   * \param  second    neutrino direction
   * \return           dot product
   */
  inline double getDot(const JNeutrinoDirection& first,
		       const JNeutrinoDirection& second)
  {
    return 
      cos(first.getZenith()) * cos(second.getZenith()) +
      sin(first.getZenith()) * sin(second.getZenith()) * 
      cos(first.getAzimuth() - second.getAzimuth());
  }


  /**
   * Dot product.
   *
   * \param  first     source location
   * \param  second    source location
   * \return           dot product
   */
  inline double getDot(const JSourceLocation& first,
		       const JSourceLocation& second)
  {
    return 
      sin(first.getDeclination()) * sin(second.getDeclination()) +
      cos(first.getDeclination()) * cos(second.getDeclination()) * 
      cos(first.getRightAscension() - second.getRightAscension());
  }


  // detector locations

  static const JGeographicalLocation Antares(42, 48, 06, 10);
  static const JGeographicalLocation Sicily (36, 16, 16, 06);
  static const JGeographicalLocation Pylos  (36, 33, 16, 06);
  static const JGeographicalLocation ARCA   (36, 17, 15, 58);
  static const JGeographicalLocation ORCA   (42, 48, 06, 02);

  // source locations

  static const JSourceLocation galacticCenter(-0.5062816, -1.633335);
  static const JSourceLocation RXJ1713       (getRadians(-39, -46,   0.0), getHourAngle(17, 13, 7));
  static const JSourceLocation VELAX         (getRadians(-45, -10, -35.2), getHourAngle( 8, 35, 20.66));
}

#endif