JAstronomy.hh

Go to the documentation of this file.

#ifndef __JASTRONOMY__JASTRONOMY__
#define __JASTRONOMY__JASTRONOMY__
 
#include <istream>
#include <ostream>
#include <limits>
#include <time.h>
 
#include "JLang/JException.hh"
#include "JLang/JManip.hh"
#include "JLang/JEquals.hh"
 
#include "JMath/JConstants.hh"
#include "JMath/JMatrix3D.hh"
#include "JMath/JSVD3D.hh"
 
#include "JUTM/JUTMGrid.hh"
 
#include "JGeometry3D/JAngle3D.hh"
#include "JGeometry3D/JDirection3D.hh"
 
#include "slalib/slalib.h"
 
/**
 * \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 {
 
  using JLANG::JValueOutOfRange;
  using JLANG::JEquals;
  using JMATH::PI;
  using JGEOMETRY3D::JAngle3D;                               //!< polar angles [rad]
  using JGEOMETRY3D::JDirection3D;                           //!< Cartesian coordinates
  using JUTM::getUTMZone;
  using JUTM::getUTMLongitude;
 
  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;
 
  static const double LATITUDE_MIN_DEG  = -80.0;             //!< Minimal latitude angle [deg] of UTM for meridian convergence angle calculation
  static const double LATITUDE_MAX_DEG  =  84.0;             //!< Maximal latitude angle [deg] of UTM for meridian convergence angle calculation
 
  static const double MJD_EPOCH   = 40587.0;                                  //!< MJD of the Unix Epoch, i.e.\ 00:00 Jan 1 1970 (UTC) [d]
  static const double MJD_EPOCH_S = MJD_EPOCH * NUMBER_OF_SECONDS_PER_DAY;    //!< MJD of the Unix Epoch, i.e.\ 00:00 Jan 1 1970 (UTC) [s]
 
 
  /**
   * Ellipsoid.
   */
  struct JEllipsoid {
    /**
     * Get ellipsoid.
     *
     * \return               Earth ellipsoid
     */
    const JEllipsoid& getEllipsoid() const
    {
      return static_cast<const JEllipsoid&>(*this);
    }
    
    double radius_m;         //!< Ellipsoid radius [m];
    double eccentricity;     //!< Ellipsoid eccentricity squared
  };
 
 
  /**
   * Earth ellipsoid according WGS84.
   */
  static const JEllipsoid EARTH_WGS84  = { 6378137.0, 0.00669438 };
  
 
  /**
   * Get MJD time given UTC time.
   *
   * \param  t1_s          UTC time [s]
   * \return               MJD time [s]
   */
  double getTimeMJD(time_t t1_s)
  {
    return t1_s + MJD_EPOCH_S;
  }
 
 
  /**
   * Get UTC time given MJD time.
   *
   * \param  t1_s          MJD time [s]
   * \return               UTC time [s]
   */
  double getTimeUTC(time_t t1_s)
  {
    return t1_s - MJD_EPOCH_S;
  }
 
 
  /**
   * Convert angle to radians.
   *
   * \param  angle         angle [deg]
   * \return               angle [rad]
   */
  inline double getRadians(const double angle)
  {
    return PI * angle / 180.0;
  }
 
 
  /**
   * Convert angle to degrees.
   *
   * \param  angle         angle [rad]
   * \return               angle [deg]
   */
  inline double getDegrees(const double angle)
  {
    return 180.0 * angle / PI;
  }
 
  
  /**
   * 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 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 = PI * ((double) hour /    12.0 + 
                      (double) min  /   720.0 + 
                      (double) sec  / 43200.0);
 
    if (ha > PI) {
      ha -= 2*PI;
    }
 
    return ha;
  }
 
 
  /**
   * Forward declarations.
   */
  struct angle_type_rad;
  struct angle_type_deg;
  struct J2000;
  struct JSourceLocation;
  struct JSourceLocationJ2000;
  struct JGalacticCoordinates;
 
  
  /**
   * Auxiliary data structure for pair of angles.
   */
  struct angle_type :
    public JEquals<angle_type>
  {
    /**
     * Check equality.
     *
     * \param  angle              pair of angles
     * \param  precision          precision
     * \return                    true if angles are equal; else false
     */
    bool equals(const angle_type& angle,
                const double      precision = std::numeric_limits<double>::min()) const
    {
      return (fabs(this->_theta_ - angle._theta_) <= precision &&
              fabs(this->_phi_   - angle._phi_)   <= precision);
    }
 
    
    /**
     * Read pair of angles from input stream.
     *
     * \param  in                 input stream
     * \param  object             pair of angles
     * \return                    input stream
     */
    friend inline std::istream& operator>>(std::istream& in, angle_type& object)
    {
      return in >> object._theta_
                >> object._phi_;
    }
    
    /**
     * Write pair of angles to output stream.
     *
     * \param  out                output stream
     * \param  object             pair of angles
     * \return                    output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const angle_type& object)
    {
      return out << FIXED(12,7) << object._theta_ << ' '
                 << FIXED(12,7) << object._phi_;
    }
    
  protected:
    /**
     * Default constructor.
     */
    angle_type() :
      _theta_(0.0),
      _phi_  (0.0)
    {}
 
 
    /**
     * Constructor.
     *
     * \param  theta              theta
     * \param  phi                phi
     */
    angle_type(const double theta,
               const double phi) :
      _theta_(theta),
      _phi_  (phi)
    {}
 
    double _theta_;
    double _phi_;
  };
 
 
  /**
   * Auxiliary data structure for pair of angles.
   */
  struct angle_type_rad :
    public angle_type
  {
    /**
     * Default constructor.
     */
    angle_type_rad()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  theta              theta [rad]
     * \param  phi                phi   [rad]
     */
    angle_type_rad(const double theta,
                   const double phi) :
      angle_type(theta, phi)
    {}
 
 
    /**
     * Conversion constructor.
     *
     * \param  angle              angle [deg]
     */
    angle_type_rad(const angle_type_deg& angle);
 
 
    /**
     * Convert angle.
     *
     * \param  angle              angle [rad]
     */
    void set(const angle_type_deg& angle)
    {
      static_cast<angle_type_rad&>(*this) = angle_type_rad(angle);
    }
 
 
    friend struct angle_type_deg;
    friend struct J2000;
  };
 
 
  /**
   * Auxiliary data structure for pair of angles.
   */
  struct angle_type_deg :
    public angle_type
  {
    /**
     * Default constructor.
     */
    angle_type_deg()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  theta              theta [deg]
     * \param  phi                phi   [deg]
     */
    angle_type_deg(const double theta,
                   const double phi) :
      angle_type(theta, phi)
    {}
 
 
    /**
     * Conversion constructor.
     *
     * \param  angle              angle [rad]
     */
    angle_type_deg(const angle_type_rad& angle);
 
 
    /**
     * Convert angle.
     *
     * \param  angle              angle [rad]
     */
    void set(const angle_type_rad& angle)
    {
      static_cast<angle_type_deg&>(*this) = angle_type_deg(angle);
    }
 
 
    friend struct angle_type_rad;
  };
 
 
  /**
   * Conversion constructor.
   *
   * \param  angle              angle [deg]
   */
  angle_type_rad::angle_type_rad(const angle_type_deg& angle) :
    angle_type(getRadians(angle._theta_), getRadians(angle._phi_))
  {}
 
 
  /**
   * Conversion constructor.
   *
   * \param  angle              angle [rad]
   */
  angle_type_deg::angle_type_deg(const angle_type_rad& angle) :
    angle_type(getDegrees(angle._theta_), getDegrees(angle._phi_))
  {}
 
 
  /**
   * Conversion of source location.
   */
  struct J2000 {
    /**
     * Conversion options.
     */
    enum CONVERSION {
      FROM = -1,
      TO   = +1
    };
 
 
    /**
     * Convert source location according J2000 at given time.
     * Adapted from <tt>aanet/astro/Astro.cc</tt>
     *
     * \param   t1_s          time since MJD [s]
     * \param   location      source location
     * \param   option        option
     * \return                source location
     */
    static inline angle_type_rad convert(const double t1_s, const angle_type_rad& location, const J2000::CONVERSION option)
    {
      using namespace JPP;
 
      // --- get the rotation matrix ----
      // NOTE: (from seatray Astro): Julian epoch of J2000 time definition is 2000,
      // corresponding by definition to 2451545.0 TT (unmodified) Julian date!
 
      double mjd   = t1_s / NUMBER_OF_SECONDS_PER_DAY;   // [days]
      double epoch = 2000;
      double R[3][3] {};
 
      slaPrenut(epoch, mjd, R);                        // SLALIB function to get the matrix of precession and nutation
 
      JMatrix3D  M  = JMatrix3D(R[0][0], R[0][1], R[0][2],
                                R[1][0], R[1][1], R[1][2],
                                R[2][0], R[2][1], R[2][2]);
 
      M.transpose();                                     // transpose 'cause of fortan matrix convention 
 
      if (option == J2000::TO) {
 
        JSVD3D V(M);
 
        M = V.invert();
      }
 
      // convert (declination, right ascension) to polar coordinates (theta, phi)
 
      JDirection3D v(JAngle3D(PI/2 - location._theta_, location._phi_));
 
      // apply transformation
      
      v.transform(M);
 
      // convert polar coordinates (theta, phi) to (declination, right ascension)
 
      return angle_type_rad(PI/2 - v.getTheta(), v.getPhi());
    }
  };
 
 
  /**
   * Location of astrophysical source.
   */
  struct JSourceLocation :
    public angle_type_rad
  {
    /**
     * Default constructor.
     */
    JSourceLocation()
    {}
 
 
    /**
     * Constructor.
     *
     * \param    dec              declination     [rad]
     * \param    ra               right ascension [rad]
     */
    JSourceLocation(const double dec,
                    const double ra) :
      angle_type_rad(dec, ra)
    {}
 
 
    /**
     * Conversion constructor.
     *
     * \param   t1_s              time since MJD [s]
     * \param   location          source location
     */
    JSourceLocation(const double t1_s, const JSourceLocationJ2000& location);
 
 
    /**
     * Constructor.
     *
     * \param    angle            polar angles [rad]
     */
    JSourceLocation(const JAngle3D& angle) :
      angle_type_rad(PI/2 - angle.getTheta(), angle.getPhi())
    {}
 
 
    /**
     * Get source location.
     *
     * \return                    source location
     */
    const JSourceLocation& getSourceLocation() const
    {
      return static_cast<const JSourceLocation&>(*this);
    }
 
 
    /**
     * Type conversion operator.
     *
     * \return                    polar angles [rad]
     */
    operator JAngle3D() const
    {
      return JAngle3D(PI/2 - _theta_, _phi_);
    }
 
 
    double getDeclination()    const { return _theta_; }        //!< Get declination.
    double getRightAscension() const { return _phi_; }          //!< Get right ascension.
 
 
    /**
     * Dot product.
     *
     * \param  location           source location
     * \return                    dot product
     */
    inline double getDot(const JSourceLocation& location) const
    {
      return 
        cos(this->_theta_) * cos(location._theta_) * cos(this->_phi_ - location._phi_) +
        sin(this->_theta_) * sin(location._theta_);
    }
  };
 
 
  /**
   * Location of astrophysical source.
   */
  struct JSourceLocationJ2000 :
    public angle_type_rad
  {
    /**
     * Default constructor.
     */
    JSourceLocationJ2000()
    {}
 
 
    /**
     * Constructor.
     *
     * \param    dec              declination     [rad]
     * \param    ra               right ascension [rad]
     */
    JSourceLocationJ2000(const double dec,
                         const double ra) :
      angle_type_rad(dec, ra)
    {}
 
 
    /**
     * Conversion constructor.
     *
     * \param   t1_s              time since MJD [s]
     * \param   location          source location
     */
    JSourceLocationJ2000(const double t1_s, const JSourceLocation& location);
 
 
    /**
     * Conversion constructor.
     *
     * \param   location          source location
     */
    JSourceLocationJ2000(const JGalacticCoordinates& location);
 
 
    /**
     * Constructor.
     *
     * \param    angle            polar angles [rad]
     */
    JSourceLocationJ2000(const JAngle3D& angle) :
      angle_type_rad(PI/2 - angle.getTheta(), angle.getPhi())
    {}
 
 
    /**
     * Get source location.
     *
     * \return                    source location
     */
    const JSourceLocationJ2000& getSourceLocation() const
    {
      return static_cast<const JSourceLocationJ2000&>(*this);
    }
 
 
    /**
     * Type conversion operator.
     *
     * \return                    polar angles [rad]
     */
    operator JAngle3D() const
    {
      return JAngle3D(PI/2 - _theta_, _phi_);
    }
 
 
    double getDeclination()    const { return _theta_; }        //!< Get declination.
    double getRightAscension() const { return _phi_; }          //!< Get right ascension.
 
 
    /**
     * Dot product.
     *
     * \param  location           source location
     * \return                    dot product
     */
    inline double getDot(const JSourceLocationJ2000& location) const
    {
      return 
        cos(this->_theta_) * cos(location._theta_) * cos(this->_phi_ - location._phi_) +
        sin(this->_theta_) * sin(location._theta_);
    }
  };
 
 
  /**
   * Conversion constructor.
   *
   * \param   t1_s              time since MJD [s]
   * \param   location          source location
   */
  JSourceLocation::JSourceLocation(const double t1_s, const JSourceLocationJ2000& location) :
    angle_type_rad(J2000::convert(t1_s, location, J2000::FROM))
  {}
 
 
  /**
   * Conversion constructor.
   *
   * \param   t1_s              time since MJD [s]
   * \param   location          source location
   */
  JSourceLocationJ2000::JSourceLocationJ2000(const double t1_s, const JSourceLocation& location) :
    angle_type_rad(J2000::convert(t1_s, location, J2000::TO))
  {}
 
 
  /**
   * Location of astrophysical source in Galactic coordinates.
   */
  struct JGalacticCoordinates :
    public angle_type_rad
  {
    /**
     * Default constructor.
     */
    JGalacticCoordinates()
    {}
 
 
    /**
     * Constructor.
     *
     * \param    latitude         Galactic latitude  [rad]
     * \param    longitude        Galactic longitude [rad]
     */
    JGalacticCoordinates(const double latitude,
                         const double longitude) :
      angle_type_rad(latitude, longitude)
    {}
 
 
    /** 
     * Constructor.
     *
     * \param   location          source location
     */
    JGalacticCoordinates(const JSourceLocationJ2000& location);
 
 
    /**
     * Constructor.
     *
     * \param    angle            polar angles [rad]
     */
    JGalacticCoordinates(const JAngle3D& angle) :
      angle_type_rad(PI/2 - angle.getTheta(), angle.getPhi())
    {}
 
 
    /**
     * Get Galactic coordinates.
     *
     * \return                    Galactic coordinates
     */
    const JGalacticCoordinates& getGalacticCoordinates() const
    {
      return static_cast<const JGalacticCoordinates&>(*this);
    }
 
 
    /**
     * Type conversion operator.
     *
     * \return                    polar angles [rad]
     */
    operator JAngle3D() const
    {
      return JAngle3D(PI/2 - _theta_, _phi_);
    }
 
 
    double getLatitude()  const { return _theta_; }        //!< Get latitude.
    double getLongitude() const { return _phi_; }          //!< Get longitude.
 
 
    /**
     * Dot product.
     *
     * \param  location           source location
     * \return                    dot product
     */
    inline double getDot(const JGalacticCoordinates& location) const
    {
      return 
        cos(this->_theta_) * cos(location._theta_) * cos(this->_phi_ - location._phi_) +
        sin(this->_theta_) * sin(location._theta_);
    }
  };
 
 
  /**
   * Conversion constructor.
   *
   * \param   location          source location
   */
  JSourceLocationJ2000::JSourceLocationJ2000(const JGalacticCoordinates& location)
  {
    double lat = location.getLatitude();
    double lon = location.getLongitude();
 
    slaGaleq(lon, lat, &this->_phi_, &this->_theta_);
  }
 
 
  /** 
   * Constructor.
   *
   * \param   location          source location
   */
  JGalacticCoordinates::JGalacticCoordinates(const JSourceLocationJ2000& location)
  {
    double dec = location.getDeclination();
    double ra  = location.getRightAscension();
 
    slaEqgal(ra, dec, &this->_phi_, &this->_theta_);
  }
 
 
  /**
   * Direction of incident neutrino.
   */
  struct JNeutrinoDirection :
    public angle_type_rad
  {
    /**
     * Default constructor.
     */
    JNeutrinoDirection()
    {}
 
 
    /**
     * Constructor.
     *
     * \param    zenith           zenith  [rad]
     * \param    azimuth          azimuth [rad]
     */
    JNeutrinoDirection(const double zenith,
                       const double azimuth) :
      angle_type_rad(zenith, azimuth)
    {}
 
 
    /**
     * Constructor.
     *
     * \param    angle            polar angles [rad]
     */
    JNeutrinoDirection(const JAngle3D& angle) :
      angle_type_rad(angle.getTheta(), angle.getPhi())
    {}
 
 
    /**
     * Get neutrino direction.
     *
     * \return                    neutrino direction
     */
    const JNeutrinoDirection& getNeutrinoDirection() const
    {
      return static_cast<const JNeutrinoDirection&>(*this);
    }
 
 
    /**
     * Type conversion operator.
     *
     * \return                    polar angles [rad]
     */
    operator JAngle3D() const
    {
      return JAngle3D(_theta_, _phi_);
    }
 
 
    double getZenith()  const { return _theta_; }        //!< Get zenith.
    double getAzimuth() const { return _phi_; }          //!< Get azimuth.
 
 
    /**
     * Dot product.
     *
     * \param  dir                neutrino direction
     * \return                    dot product
     */
    inline double getDot(const JNeutrinoDirection& dir) const
    {
      return 
        sin(this->_theta_) * sin(dir._theta_) * cos(this->_phi_ - dir._phi_) +
        cos(this->_theta_) * cos(dir._theta_);
    }
  };
 
 
  /**
   * Location of detector.
   */
  struct JGeographicalLocation :
    public angle_type_rad
  {
    /**
     * Default constructor.
     */
    JGeographicalLocation()
    {}
 
 
    /**
     * Constructor.
     *
     * \param    latitude         latitude 
     * \param    longitude        longitude
     */
    JGeographicalLocation(const double latitude,
                          const double longitude) :
      angle_type_rad(latitude, 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) :
      angle_type_rad(getRadians((double) degreesNorth + (double) minutesNorth / 60.0),
                     getRadians((double) degreesEast  + (double) minutesEast  / 60.0))
    {}
 
 
    /**
     * Constructor.
     *
     * \param    angle            polar angles [rad]
     */
    JGeographicalLocation(const JAngle3D& angle) :
      angle_type_rad(PI/2 - angle.getTheta(), angle.getPhi())
    {}
 
 
    /**
     * Get geographical location.
     *
     * \return                    geographical location
     */
    const JGeographicalLocation& getGeographicalLocation() const
    {
      return static_cast<const JGeographicalLocation&>(*this);
    }
 
 
    /**
     * Type conversion operator.
     *
     * \return                    polar angles [rad]
     */
    operator JAngle3D() const
    {
      return JAngle3D(PI/2 - _theta_, _phi_);
    }
 
 
    double getLatitude()  const { return _theta_; }        //!< Get latitude.
    double getLongitude() const { return _phi_; }          //!< Get longitude.
 
 
    /**
     * Dot product.
     *
     * \param  location           location
     * \return                    dot product
     */
    inline double getDot(const JGeographicalLocation& location) const
    {
      return
        sin(this->_theta_) * sin(location._theta_) +
        cos(this->_theta_) * cos(location._theta_) * cos(this->_phi_ - location._phi_);
    }
  };
 
 
  /**
   * Get Meridian convergence angle
   *
   * The meridian convergence angle is the angle between the true meridian of a point on the surface of the earth's ellipsoid and
   * the Prime Meridian of the projection zone where the point is located, that is, the angle between the coordinate north and
   * the true north of the Gaussian plane rectangular coordinate system.
   *  see https://gis.stackexchange.com/questions/115531/calculating-grid-convergence-true-north-to-grid-north
   *
   * \param  location           geographical location
   * \return                    meridian convergence angle [rad]      
   */
  inline double getMeridianConvergenceAngle(const JGeographicalLocation& location)
  {
    const double omega = location.getLongitude() - getUTMLongitude(getUTMZone(location.getLongitude()));
    const double gamma = atan(tan(omega) * sin(location.getLatitude()));
 
    return gamma;
  }
 
 
  /**
   * Get 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 <tt>aanet/evt/Det.cc</tt>.
   *
   * \param  location           geographical location
   * \param  ellipsoid          Earth ellipsoid
   * \return                    meridian convergence angle [rad]      
   */
  inline double getMeridianConvergenceAngle(const JGeographicalLocation& location, const JEllipsoid& ellipsoid)
  {
    const double latitude_deg = getDegrees(location.getLatitude());
 
    if (latitude_deg > LATITUDE_MAX_DEG ||
        latitude_deg < LATITUDE_MIN_DEG) {
      THROW(JValueOutOfRange,
            "UTM coordinate system not defined for given Latitude: "
            << latitude_deg << " [deg] <> [" << FIXED(5,2) << LATITUDE_MAX_DEG << "," << FIXED(5,2) << LATITUDE_MIN_DEG << "]");
    }
 
    // detector position, longitude and latitude in rad
    const double lambda  = location.getLongitude();
    const double phi     = location.getLatitude();
 
    // find UTM zone and central meridian
    // longitude of the central meridian of UTM zone in rad
    const double omega   = lambda - getUTMLongitude(getUTMZone(location.getLongitude()));
 
    const double rho = ellipsoid.radius_m*(1.-ellipsoid.eccentricity)/pow(1.-ellipsoid.eccentricity*pow(sin(phi),2),3./2.);
    const double nu  = ellipsoid.radius_m/sqrt(1-ellipsoid.eccentricity*pow(sin(phi),2.));
    const double psi = nu/rho;
    const double t   = tan(phi);
 
    // power series for convergence angle 
    const 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;
  }
  
 
  /**
   * Auxiliary class to make coordinate transformations for a specific geographical location of the detector.
   *
   * Note that SLALIB reference system corresponds to (x,y,z) = (N,E,up).
   */
  class JAstronomy : 
    public JGeographicalLocation,
    public JEllipsoid
  {
  public:
 
    /**
     * Constructor.
     *
     * \param  location           location of detector
     * \param  ellipsoid          Earth ellipsoid
     */
    JAstronomy(const JGeographicalLocation& location, const JEllipsoid& ellipsoid = EARTH_WGS84) :
      JGeographicalLocation(location),
      JEllipsoid(ellipsoid)
    {}
 
 
    /**
     * Get direction pointing to source given time and source location.
     *
     * \param   t1_s              time since MJD [s]
     * \param   pos               source location
     * \return                    direction of neutrino
     */
    JNeutrinoDirection getNeutrinoDirection(const double t1_s, const JSourceLocation& pos) const
    {
      double longitude = this->getLongitude();
      double latitude  = this->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 mjd   = t1_s / NUMBER_OF_SECONDS_PER_DAY;   // [days]
 
      double gmst  = slaGmst(mjd);                       // gmst = Greenwich mean sidereal time
      double eqeqx = slaEqeqx(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;
 
      slaDe2h(ha, dec, latitude, &azimuth, &elevation);
      
      double theta = -elevation + PI/2.0;
      double phi   = -azimuth   + PI/2.0;
 
      phi += getMeridianConvergenceAngle(this->getGeographicalLocation(), this->getEllipsoid());
 
      // invert direction
 
      theta = PI  - theta;
      phi   = phi + PI;
 
      if (phi > PI) {
        phi -= 2.0*PI;
      }
 
      return JNeutrinoDirection(theta, phi);
    }
 
 
    /**
     *  Get location of source given time and neutrino direction.
     *
     * \param   t1_s              time since MJD [s]
     * \param   dir               direction of neutrino
     * \return                    source location  
     */
    JSourceLocation getSourceLocation(const double t1_s, const JNeutrinoDirection& dir) const
    {
      double longitude = this->getLongitude();
      double latitude  = this->getLatitude();
 
      double theta = dir.getZenith(); 
      double phi   = dir.getAzimuth(); 
      
      // invert direction
 
      theta = PI  - theta;
      phi   = phi - PI;
 
      if (phi < PI) {
        phi += 2.0*PI;
      }
 
      phi -= getMeridianConvergenceAngle(this->getGeographicalLocation(), this->getEllipsoid());
 
      double elevation = -theta + PI/2.0;
      double azimuth   = -phi   + PI/2.0;
 
      double ha;
      double dec;
 
      slaDh2e(azimuth, elevation, latitude, &ha, &dec);
 
      double mjd   = t1_s / NUMBER_OF_SECONDS_PER_DAY;   // [days]
 
      double gmst  = slaGmst(mjd);                       // gmst = Greenwich mean sidereal time
      double eqeqx = slaEqeqx(mjd);                      // Note: Equation of the Equinoxes < 1.15 secs
      double lst   = gmst + longitude + eqeqx;
      double ra    = lst - ha;
      
      return JSourceLocation(dec, ra);
    }
  };
 
 
  // 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   (0.63342, 0.27883);
  static const JGeographicalLocation ORCA   (0.74702, 0.10511);
 
  // source locations
 
  static const JSourceLocation GALACTIC_CENTER(-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