#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