JAstronomy.hh

Go to the documentation of this file.

#ifndef __JASTRONOMY__
#define __JASTRONOMY__

#include "JMath/JConstants.hh"

/**
 * \file
 * 
 * Interface methods for slalib and auxiliary classes and methods for astronomy.
 * \author mdejong
 */
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);
  }


  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;

  
  /**
   * 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; }

  protected:
    double __declination;
    double __right_ascension;
  };


  /**
   * 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.
     *
     * \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);


  // 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