NOVASLocate.cxx

Go to the documentation of this file.

#include "CelestialLocator/NOVASLocate.h"
#include "NovaDAQConventions/DAQConventions.h"

#include "cetlib/search_path.h"

#include <cassert>

// TODO - all this NOVAS code should be a ups product

// Old C code triggers some warnings we can't do anything about
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wwrite-strings"
#pragma GCC diagnostic ignored "-Wpedantic"

#include "novas/novas.c"
#include "novas/eph_manager.c"

#include "novas/novascon.c"
#include "novas/solsys1.c"
#include "novas/nutation.c"
#include "novas/readeph0.c"

#pragma  GCC diagnostic pop


namespace shadow
{
  int NOVASLocate::fgRefCount = 0;

  //----------------------------------------------------------------------
  NOVASLocate::NOVASLocate(const std::string& ephemerisFname)
  {
    ++fgRefCount;

    // Dummy temperature and pressure
    make_on_surface(latitude, longitude, height, 0, 0, &fGeoLoc);

    cat_entry dummy_star;
    #pragma GCC diagnostic push
    #pragma GCC diagnostic ignored "-Wpedantic"
    
    std::string str1("DUMMY");
    std::string str2("xxx");
    
    int error = make_cat_entry(str1.c_str(), str2.c_str(),
                               0, 0, 0, 0, 0, 0, 0, &dummy_star);
    assert(error == 0);

    str1 = "Moon";
    str2 = "Sun";
    error = make_object(0, 11, str1.c_str(), &dummy_star, &fMoon);
    assert(error == 0);

    error = make_object(0, 10, str2.c_str(), &dummy_star, &fSun);
    assert(error == 0);
    #pragma GCC diagnostic pop

    double jd_beg, jd_end;
    short int de_num = 0;
    error = ephem_open((char*)ephemerisFname.c_str(), &jd_beg, &jd_end, &de_num);
    assert(error == 0);

     #if 0
     /* Test code: see if you put angles in and convert them back     */
     /* and forth, you get back where you started. You can also       */
     /* externally verify the intermediat angles.                     */

     // zenith is down from the zenith (obviously)
     // azimuth is *east* of *north* (astronomers have it backwards)
     const double zenin = 32.109 * M_PI/180, aziin = 123.45 * M_PI/180;

     const double jd = julian_date(2019, 1, 21, 1);
     printf("Julian date = %f\n", jd);
    
     double ra, dec;
     GetRaDec(zenin, aziin, jd, ra, dec);
     printf("ra = %f, dec = %f\n", ra, dec);
     printf("Or: ra = %fh, dec = %f degrees\n", ra*12/M_PI, dec*180/M_PI);

     // verify ra and dec from another conversion program
     printf("other diff: ra = %fh, dec = %f degrees\n",
             ra*12/M_PI - 4 - 10/60. - 58/60./60.,
             dec*180/M_PI - 23 -  9/60. - 49/60./60.);

     double zenout, aziout;
     GetZenAzi(ra, dec, jd, zenout, aziout);
     printf("zenout = %f (delta %+f), aziout = %f (%+f)\n",
            zenout, zenout - zenin, aziout, aziout - aziin);
     printf("In degrees: zenout = %f (delta %+f), aziout = %f (%+f)\n\n",
            zenout*180/M_PI, (zenout - zenin)*180/M_PI,
            aziout*180/M_PI, (aziout - aziin)*180/M_PI);
     #endif
  }

  //----------------------------------------------------------------------
  NOVASLocate::~NOVASLocate()
  {
    --fgRefCount;

    if(fgRefCount == 0){
      ephem_close();
    }
  }

  //----------------------------------------------------------------------
  void NOVASLocate::GetUT1andDeltat(const double jd_utc,
                                    double & jd_ut1, double & delta_t)
  {
    // XXX This is super fragile.  However, it is also hard to see when
    // on NOvA we'd need enough precision that it matters.
    const double ut1_utc = -0.387845;
    jd_ut1 = jd_utc + ut1_utc / 86400.0;
    delta_t = 32.184 + leap_secs - ut1_utc;
  }

  //----------------------------------------------------------------------
  void NOVASLocate::GetObjectPosition(object& obj,
                                      double jd_utc,
                                      double& zen, double& azi)
  {
    double jd_ut1, delta_t;
    GetUT1andDeltat(jd_utc, jd_ut1, delta_t);
    const double jd_tt = jd_utc + (leap_secs + 32.184) / 86400.0;

    double rat, dect, dist;

    // Coordinate system used in returned ra/dec is "true equator and
    // equinox of date" (a term which I don't really understand)
    int error = topo_planet(jd_tt, &obj, delta_t, &fGeoLoc, 0,
                            &rat, &dect, &dist);
    assert(error == 0);

    // Position of the Moon in local horizon coordinates.  (Polar motion
    // ignored here.)

    double rar, decr;
    equ2hor(jd_ut1,
            delta_t,
            0, // full accuracy
            0, // xp: don't need this much precision
            0, // yp: ditto
            &fGeoLoc,
            rat, dect,
            0, // no refraction
            &zen, &azi,
            &rar, &decr);
  }

  //----------------------------------------------------------------------
  void NOVASLocate::GetRaDec(const double zen, const double azi,
                             double jd_utc,
                             double& ra, double& dec)
  {
    double jd_ut1, delta_t;
    GetUT1andDeltat(jd_utc, jd_ut1, delta_t);

    // Turn zenith and azimuth into (x,y,z), where z means *up* (not north)
    const double dir[3] = { cos(M_PI - azi)*cos(M_PI_2 - zen),
                            sin(M_PI - azi)*cos(M_PI_2 - zen),
                                            sin(M_PI_2 - zen) };

    // Change in angles between local and terrestrial.  If the detector
    // were at the north pole with the x-axis pointed towards Greenwich,
    // these would be zero.
    const double dtheta = (90-latitude)*M_PI/180, dphi = -longitude*M_PI/180;

    // Turn (x,y,z), where z is relative to the surface, into (x,y,z) in 
    // ITRS coordinates, where z is parallel to the north pole and x is in
    // the direction of the prime meridian (?).
    double itrs[3] = {
     cos(dtheta)*cos(dphi)*dir[0] + sin(dphi)*dir[1] + sin(dtheta)*cos(dphi)*dir[2],
    -cos(dtheta)*sin(dphi)*dir[0] + cos(dphi)*dir[1] - sin(dtheta)*sin(dphi)*dir[2],
    -sin(dtheta)*          dir[0]                    + cos(dtheta)          *dir[2]
    };

    // Finally have NOVAS turn this into ra and dec in GCRS.
    double gcrs[3] = {0};
    ter2cel(jd_ut1,
            0, // no need for very precise time
            delta_t,
            1, // Method: 0 or 1 makes little difference
            0, // full accuracy
            1, // Want equator and equinox of date
            0, // xp: don't need high precision
            0, // yp: ditto
            itrs, gcrs);
    ra  = atan2(gcrs[1], gcrs[0]);
    dec = asin(gcrs[2]);
  }

  //----------------------------------------------------------------------
  void NOVASLocate::GetTrackRaDec(const TVector3 & trkdir,
                                  double jd_utc,
                                  double & ra, double & dec)
  {
    const double zen_trk = acos(trkdir.Y());      

    // Chris Backhouse says:
    // docdb 5485 / numix 17 says FD is oriented -27o51'26'' from true North
    //      gTrkAzi = atan2(-dir.X(), dir.Z()) * 180 / M_PI - 27.857222;
    // Or: 332o03'58.071769" (clockwise from North). Email from Virgil Bocean
    // to Alec Habig.
    const double azi_trk = atan2(-trkdir.X(), trkdir.Z()) + point * M_PI/180;

    GetRaDec(zen_trk, azi_trk, jd_utc, ra, dec);
  }

  //----------------------------------------------------------------------
  void NOVASLocate::GetZenAzi(const double rat, const double dect,
                              double jd_utc,
                              double& zen, double& azi)
  {
    double jd_ut1, delta_t;
    GetUT1andDeltat(jd_utc, jd_ut1, delta_t);

    double rar, decr;
    equ2hor(jd_ut1, delta_t,
            0, // full accuracy
            0, // xp: set to zero because we do not need arcsecond precision
            0, // yp: ditto
            &fGeoLoc,
            rat  *  12/M_PI, // in *hours*
            dect * 180/M_PI, // in *degrees*
            0, // no refraction
            &zen, &azi,
            &rar, &decr);
    zen *= M_PI/180;
    azi *= M_PI/180;
  }

  //----------------------------------------------------------------------
  double NOVASLocate::JulianDate(time_t t)
  {
    const tm* t2 = gmtime(&t);
    const short int year  = t2->tm_year + 1900;
    const short int month = t2->tm_mon + 1;
    const short int day   = t2->tm_mday;
    const double    hour  = t2->tm_hour + t2->tm_min / 60. + t2->tm_sec / (60. * 60.);

    return julian_date(year, month, day, hour);
  }

  //----------------------------------------------------------------------
  void NOVASLocate::GetMoonPosition(time_t t,
                                    double& zen, double& azi)
  {
    GetMoonPosition(JulianDate(t), zen, azi);
  }

  //----------------------------------------------------------------------
  void NOVASLocate::GetMoonPosition(double jd_utc,
                                    double& zen, double& azi)
  {
    GetObjectPosition(fMoon, jd_utc, zen, azi);
  }

  //----------------------------------------------------------------------
  void NOVASLocate::GetSunPosition(time_t t,
                                   double& zen, double& azi)
  {
    GetSunPosition(JulianDate(t), zen, azi);
  }

  //----------------------------------------------------------------------
  void NOVASLocate::GetSunPosition(double jd_utc,
                                    double& zen, double& azi)
  {
    GetObjectPosition(fSun, jd_utc, zen, azi);
  }

  void NOVASLocate::SetDetector(const novadaq::cnv::DetId det)
  {
    switch(det){
      case novadaq::cnv::DetId::kFARDET:
        latitude = FDlatitude;
        longitude = FDlongitude;
        height = FDheight;
        point = FDpoint;
        break;
      case novadaq::cnv::DetId::kNEARDET:
        latitude = NDlatitude;
        longitude = NDlongitude;
        height = NDheight;
        point = NDpoint;
        break;
      case novadaq::cnv::DetId::kTESTBEAM:
        latitude = TBlatitude;
        longitude = TBlongitude;
        height = TBheight;
        point = TBpoint;
        break;
      default:
        fprintf(stderr, "CelestialLocator: unknown detector %d\n", det);
        abort();
    }
    make_on_surface(latitude, longitude, height, 0, 0, &fGeoLoc);
  }
} // namespace