NOVASLocate.cxx
Go to the documentation of this file.
2 #include "NovaDAQConventions/DAQConventions.h"
3 
4 #include "cetlib/search_path.h"
5 
6 #include <cassert>
7 
8 // TODO - all this NOVAS code should be a ups product
9 
10 // Old C code triggers some warnings we can't do anything about
11 #pragma GCC diagnostic push
12 #pragma GCC diagnostic ignored "-Wwrite-strings"
13 #pragma GCC diagnostic ignored "-Wpedantic"
14 
15 #include "novas/novas.c"
16 #include "novas/eph_manager.c"
17 
18 #include "novas/novascon.c"
19 #include "novas/solsys1.c"
20 #include "novas/nutation.c"
21 #include "novas/readeph0.c"
22 
23 #pragma GCC diagnostic pop
24 
25 
26 namespace shadow
27 {
29 
30  //----------------------------------------------------------------------
31  NOVASLocate::NOVASLocate(const std::string& ephemerisFname)
32  {
33  ++fgRefCount;
34 
35  // Dummy temperature and pressure
37 
38  cat_entry dummy_star;
39  #pragma GCC diagnostic push
40  #pragma GCC diagnostic ignored "-Wpedantic"
41 
42  std::string str1("DUMMY");
43  std::string str2("xxx");
44 
45  int error = make_cat_entry(str1.c_str(), str2.c_str(),
46  0, 0, 0, 0, 0, 0, 0, &dummy_star);
47  assert(error == 0);
48 
49  str1 = "Moon";
50  str2 = "Sun";
51  error = make_object(0, 11, str1.c_str(), &dummy_star, &fMoon);
52  assert(error == 0);
53 
54  error = make_object(0, 10, str2.c_str(), &dummy_star, &fSun);
55  assert(error == 0);
56  #pragma GCC diagnostic pop
57 
58  double jd_beg, jd_end;
59  short int de_num = 0;
60  error = ephem_open((char*)ephemerisFname.c_str(), &jd_beg, &jd_end, &de_num);
61  assert(error == 0);
62 
63  #if 0
64  /* Test code: see if you put angles in and convert them back */
65  /* and forth, you get back where you started. You can also */
66  /* externally verify the intermediat angles. */
67 
68  // zenith is down from the zenith (obviously)
69  // azimuth is *east* of *north* (astronomers have it backwards)
70  const double zenin = 32.109 * M_PI/180, aziin = 123.45 * M_PI/180;
71 
72  const double jd = julian_date(2019, 1, 21, 1);
73  printf("Julian date = %f\n", jd);
74 
75  double ra, dec;
76  GetRaDec(zenin, aziin, jd, ra, dec);
77  printf("ra = %f, dec = %f\n", ra, dec);
78  printf("Or: ra = %fh, dec = %f degrees\n", ra*12/M_PI, dec*180/M_PI);
79 
80  // verify ra and dec from another conversion program
81  printf("other diff: ra = %fh, dec = %f degrees\n",
82  ra*12/M_PI - 4 - 10/60. - 58/60./60.,
83  dec*180/M_PI - 23 - 9/60. - 49/60./60.);
84 
85  double zenout, aziout;
86  GetZenAzi(ra, dec, jd, zenout, aziout);
87  printf("zenout = %f (delta %+f), aziout = %f (%+f)\n",
88  zenout, zenout - zenin, aziout, aziout - aziin);
89  printf("In degrees: zenout = %f (delta %+f), aziout = %f (%+f)\n\n",
90  zenout*180/M_PI, (zenout - zenin)*180/M_PI,
91  aziout*180/M_PI, (aziout - aziin)*180/M_PI);
92  #endif
93  }
94 
95  //----------------------------------------------------------------------
97  {
98  --fgRefCount;
99 
100  if(fgRefCount == 0){
101  ephem_close();
102  }
103  }
104 
105  //----------------------------------------------------------------------
106  void NOVASLocate::GetUT1andDeltat(const double jd_utc,
107  double & jd_ut1, double & delta_t)
108  {
109  // XXX This is super fragile. However, it is also hard to see when
110  // on NOvA we'd need enough precision that it matters.
111  const double ut1_utc = -0.387845;
112  jd_ut1 = jd_utc + ut1_utc / 86400.0;
113  delta_t = 32.184 + leap_secs - ut1_utc;
114  }
115 
116  //----------------------------------------------------------------------
118  double jd_utc,
119  double& zen, double& azi)
120  {
121  double jd_ut1, delta_t;
122  GetUT1andDeltat(jd_utc, jd_ut1, delta_t);
123  const double jd_tt = jd_utc + (leap_secs + 32.184) / 86400.0;
124 
125  double rat, dect, dist;
126 
127  // Coordinate system used in returned ra/dec is "true equator and
128  // equinox of date" (a term which I don't really understand)
129  int error = topo_planet(jd_tt, &obj, delta_t, &fGeoLoc, 0,
130  &rat, &dect, &dist);
131  assert(error == 0);
132 
133  // Position of the Moon in local horizon coordinates. (Polar motion
134  // ignored here.)
135 
136  double rar, decr;
137  equ2hor(jd_ut1,
138  delta_t,
139  0, // full accuracy
140  0, // xp: don't need this much precision
141  0, // yp: ditto
142  &fGeoLoc,
143  rat, dect,
144  0, // no refraction
145  &zen, &azi,
146  &rar, &decr);
147  }
148 
149  //----------------------------------------------------------------------
150  void NOVASLocate::GetRaDec(const double zen, const double azi,
151  double jd_utc,
152  double& ra, double& dec)
153  {
154  double jd_ut1, delta_t;
155  GetUT1andDeltat(jd_utc, jd_ut1, delta_t);
156 
157  // Turn zenith and azimuth into (x,y,z), where z means *up* (not north)
158  const double dir[3] = { cos(M_PI - azi)*cos(M_PI_2 - zen),
159  sin(M_PI - azi)*cos(M_PI_2 - zen),
160  sin(M_PI_2 - zen) };
161 
162  // Change in angles between local and terrestrial. If the detector
163  // were at the north pole with the x-axis pointed towards Greenwich,
164  // these would be zero.
165  const double dtheta = (90-latitude)*M_PI/180, dphi = -longitude*M_PI/180;
166 
167  // Turn (x,y,z), where z is relative to the surface, into (x,y,z) in
168  // ITRS coordinates, where z is parallel to the north pole and x is in
169  // the direction of the prime meridian (?).
170  double itrs[3] = {
171  cos(dtheta)*cos(dphi)*dir[0] + sin(dphi)*dir[1] + sin(dtheta)*cos(dphi)*dir[2],
172  -cos(dtheta)*sin(dphi)*dir[0] + cos(dphi)*dir[1] - sin(dtheta)*sin(dphi)*dir[2],
173  -sin(dtheta)* dir[0] + cos(dtheta) *dir[2]
174  };
175 
176  // Finally have NOVAS turn this into ra and dec in GCRS.
177  double gcrs[3] = {0};
178  ter2cel(jd_ut1,
179  0, // no need for very precise time
180  delta_t,
181  1, // Method: 0 or 1 makes little difference
182  0, // full accuracy
183  1, // Want equator and equinox of date
184  0, // xp: don't need high precision
185  0, // yp: ditto
186  itrs, gcrs);
187  ra = atan2(gcrs[1], gcrs[0]);
188  dec = asin(gcrs[2]);
189  }
190 
191  //----------------------------------------------------------------------
192  void NOVASLocate::GetTrackRaDec(const TVector3 & trkdir,
193  double jd_utc,
194  double & ra, double & dec)
195  {
196  const double zen_trk = acos(trkdir.Y());
197 
198  // Chris Backhouse says:
199  // docdb 5485 / numix 17 says FD is oriented -27o51'26'' from true North
200  // gTrkAzi = atan2(-dir.X(), dir.Z()) * 180 / M_PI - 27.857222;
201  // Or: 332o03'58.071769" (clockwise from North). Email from Virgil Bocean
202  // to Alec Habig.
203  const double azi_trk = atan2(-trkdir.X(), trkdir.Z()) + point * M_PI/180;
204 
205  GetRaDec(zen_trk, azi_trk, jd_utc, ra, dec);
206  }
207 
208  //----------------------------------------------------------------------
209  void NOVASLocate::GetZenAzi(const double rat, const double dect,
210  double jd_utc,
211  double& zen, double& azi)
212  {
213  double jd_ut1, delta_t;
214  GetUT1andDeltat(jd_utc, jd_ut1, delta_t);
215 
216  double rar, decr;
217  equ2hor(jd_ut1, delta_t,
218  0, // full accuracy
219  0, // xp: set to zero because we do not need arcsecond precision
220  0, // yp: ditto
221  &fGeoLoc,
222  rat * 12/M_PI, // in *hours*
223  dect * 180/M_PI, // in *degrees*
224  0, // no refraction
225  &zen, &azi,
226  &rar, &decr);
227  zen *= M_PI/180;
228  azi *= M_PI/180;
229  }
230 
231  //----------------------------------------------------------------------
232  double NOVASLocate::JulianDate(time_t t)
233  {
234  const tm* t2 = gmtime(&t);
235  const short int year = t2->tm_year + 1900;
236  const short int month = t2->tm_mon + 1;
237  const short int day = t2->tm_mday;
238  const double hour = t2->tm_hour + t2->tm_min / 60. + t2->tm_sec / (60. * 60.);
239 
240  return julian_date(year, month, day, hour);
241  }
242 
243  //----------------------------------------------------------------------
245  double& zen, double& azi)
246  {
247  GetMoonPosition(JulianDate(t), zen, azi);
248  }
249 
250  //----------------------------------------------------------------------
251  void NOVASLocate::GetMoonPosition(double jd_utc,
252  double& zen, double& azi)
253  {
254  GetObjectPosition(fMoon, jd_utc, zen, azi);
255  }
256 
257  //----------------------------------------------------------------------
259  double& zen, double& azi)
260  {
261  GetSunPosition(JulianDate(t), zen, azi);
262  }
263 
264  //----------------------------------------------------------------------
265  void NOVASLocate::GetSunPosition(double jd_utc,
266  double& zen, double& azi)
267  {
268  GetObjectPosition(fSun, jd_utc, zen, azi);
269  }
270 
272  {
273  switch(det){
277  height = FDheight;
278  point = FDpoint;
279  break;
283  height = NDheight;
284  point = NDpoint;
285  break;
289  height = TBheight;
290  point = TBpoint;
291  break;
292  default:
293  fprintf(stderr, "CelestialLocator: unknown detector %d\n", det);
294  abort();
295  }
297  }
298 } // namespace
const double TBpoint
Definition: NOVASLocate.h:93
short int ephem_open(char *ephem_name, double *jd_begin, double *jd_end, short int *de_number)
const double TBlatitude
Definition: NOVASLocate.h:90
short int ephem_close(void)
short int make_object(short int type, short int number, const char name[SIZE_OF_OBJ_NAME], cat_entry *star_data, object *cel_obj)
void GetUT1andDeltat(const double jd_utc, double &jd_ut1, double &deltat)
const double TBlongitude
Definition: NOVASLocate.h:91
const double FDlatitude
Definition: NOVASLocate.h:69
short int topo_planet(double jd_tt, object *ss_body, double delta_t, on_surface *position, short int accuracy, double *ra, double *dec, double *dis)
T acos(T number)
Definition: d0nt_math.hpp:54
#define M_PI
Definition: SbMath.h:34
double dist
Definition: runWimpSim.h:113
short int ter2cel(double jd_ut_high, double jd_ut_low, double delta_t, short int method, short int accuracy, short int option, double xp, double yp, double *vec1, double *vec2)
void SetDetector(const novadaq::cnv::DetId det)
const double TBheight
Definition: NOVASLocate.h:92
#define M_PI_2
Definition: SbMath.h:35
const double NDlongitude
Definition: NOVASLocate.h:84
short int make_cat_entry(const char star_name[SIZE_OF_OBJ_NAME], const char catalog[SIZE_OF_CAT_NAME], long int star_num, double ra, double dec, double pm_ra, double pm_dec, double parallax, double rad_vel, cat_entry *star)
void GetMoonPosition(time_t t, double &zen, double &azi)
double julian_date(short int year, short int month, short int day, double hour)
const double NDpoint
Definition: NOVASLocate.h:86
void GetTrackRaDec(const TVector3 &trkdir, double jd_utc, double &ra, double &dec)
printf("%d Experimental points found\n", nlines)
void equ2hor(double jd_ut1, double delta_t, short int accuracy, double xp, double yp, on_surface *location, double ra, double dec, short int ref_option, double *zd, double *az, double *rar, double *decr)
void make_on_surface(double latitude, double longitude, double height, double temperature, double pressure, on_surface *obs_surface)
const double NDheight
Definition: NOVASLocate.h:85
double t2
static double JulianDate(time_t t)
const double FDheight
Definition: NOVASLocate.h:72
const double FDlongitude
Definition: NOVASLocate.h:70
static const short int leap_secs
Definition: NOVASLocate.h:111
const double FDpoint
Definition: NOVASLocate.h:79
void GetSunPosition(time_t t, double &zen, double &azi)
on_surface fGeoLoc
Definition: NOVASLocate.h:61
NOVASLocate(const std::string &ephemerisFname)
Definition: NOVASLocate.cxx:31
T sin(T number)
Definition: d0nt_math.hpp:132
void GetRaDec(const double zen, const double azi, double jd_utc, double &ra, double &dec)
static int fgRefCount
Definition: NOVASLocate.h:65
TDirectory * dir
Definition: macro.C:5
void GetZenAzi(const double ra, const double dec, double jd_utc, double &zen, double &azi)
T cos(T number)
Definition: d0nt_math.hpp:78
assert(nhit_max >=nhit_nbins)
void GetObjectPosition(object &obj, double jd_utc, double &zen, double &azi)
Definition: novas.h:78
const double NDlatitude
Definition: NOVASLocate.h:83
T asin(T number)
Definition: d0nt_math.hpp:60
T atan2(T number)
Definition: d0nt_math.hpp:72
static constexpr Double_t year
Definition: Munits.h:185
enum BeamMode string