GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCOMOads.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCOMOads.cpp - COMPTEL Orbit Aspect Data container class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2017-2019 by Juergen Knodlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GCOMOads.hpp
23  * @brief COMPTEL Orbit Aspect Data container class implementation
24  * @author Juergen Knodlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GException.hpp"
32 #include "GMath.hpp"
33 #include "GFits.hpp"
34 #include "GFitsTable.hpp"
35 #include "GCOMOads.hpp"
36 #include "GCOMSupport.hpp"
37 
38 /* __ Method name definitions ____________________________________________ */
39 #define G_AT "GCOMOads::at(int&)"
40 #define G_INSERT "GCOMOads::insert(int&, GCOMOad&)"
41 #define G_REMOVE "GCOMOads::remove(int&)"
42 #define G_READ "GCOMOads::read(GFitsTable&)"
43 
44 /* __ Macros _____________________________________________________________ */
45 
46 /* __ Coding definitions _________________________________________________ */
47 //#define G_GEORAD_WARNING
48 
49 /* __ Debug definitions __________________________________________________ */
50 
51 
52 /*==========================================================================
53  = =
54  = Constructors/destructors =
55  = =
56  ==========================================================================*/
57 
58 /***********************************************************************//**
59  * @brief Void constructor
60  *
61  * Constructs an empty Orbit Aspect Data container
62  ***************************************************************************/
64 {
65  // Initialise class members
66  init_members();
67 
68  // Return
69  return;
70 }
71 
72 
73 /***********************************************************************//**
74  * @brief Filename constructor
75  *
76  * @param[in] filename Orbit Aspect Data FITS file
77  *
78  * Constructs an Orbit Aspect Data container from an OAD FITS file.
79  ***************************************************************************/
80 GCOMOads::GCOMOads(const GFilename& filename)
81 {
82  // Initialise class members
83  init_members();
84 
85  // Load data
86  load(filename);
87 
88  // Return
89  return;
90 }
91 
92 
93 /***********************************************************************//**
94  * @brief Copy constructor
95  *
96  * @param[in] oads COMPTEL Orbit Aspect Data container.
97  ***************************************************************************/
99 {
100  // Initialise class members
101  init_members();
102 
103  // Copy members
104  copy_members(oads);
105 
106  // Return
107  return;
108 }
109 
110 
111 /***********************************************************************//**
112  * @brief Destructor
113  ***************************************************************************/
115 {
116  // Free members
117  free_members();
118 
119  // Return
120  return;
121 }
122 
123 
124 /*==========================================================================
125  = =
126  = Operators =
127  = =
128  ==========================================================================*/
129 
130 /***********************************************************************//**
131  * @brief Assignment operator
132  *
133  * @param[in] oads COMPTEL Orbit Aspect Data container.
134  * @return COMPTEL Orbit Aspect Data container.
135  ***************************************************************************/
137 {
138  // Execute only if object is not identical
139  if (this != &oads) {
140 
141  // Free members
142  free_members();
143 
144  // Initialise private members
145  init_members();
146 
147  // Copy members
148  copy_members(oads);
149 
150  } // endif: object was not identical
151 
152  // Return this object
153  return *this;
154 }
155 
156 
157 /*==========================================================================
158  = =
159  = Public methods =
160  = =
161  ==========================================================================*/
162 
163 /***********************************************************************//**
164  * @brief Clear COMPTEL Orbit Aspect Data container
165  ***************************************************************************/
166 void GCOMOads::clear(void)
167 {
168  // Free members
169  free_members();
170 
171  // Initialise private members
172  init_members();
173 
174  // Return
175  return;
176 }
177 
178 
179 /***********************************************************************//**
180  * @brief Clone COMPTEL Orbit Aspect Data container
181  *
182  * @return Pointer to deep copy of COMPTEL Orbit Aspect Data container.
183  ***************************************************************************/
185 {
186  return new GCOMOads(*this);
187 }
188 
189 
190 /***********************************************************************//**
191  * @brief Return reference to Orbit Aspect Data
192  *
193  * @param[in] index Orbit Aspect Data index [0,...,size()-1].
194  *
195  * @exception GException::out_of_range
196  * Orbit Aspect Data index is out of range.
197  *
198  * Returns a reference to the Orbit Aspect Data with the specified @p index.
199  ***************************************************************************/
200 GCOMOad& GCOMOads::at(const int& index)
201 {
202  // Raise exception if index is out of range
203  if (index < 0 || index >= size()) {
204  throw GException::out_of_range(G_AT, "Orbit Aspect Data index",
205  index, size());
206  }
207 
208  // Return reference
209  return m_oads[index];
210 }
211 
212 
213 /***********************************************************************//**
214  * @brief Return reference to Orbit Aspect Data (const version)
215  *
216  * @param[in] index Orbit Aspect Data index [0,...,size()-1].
217  *
218  * @exception GException::out_of_range
219  * Orbit Aspect Data index is out of range.
220  *
221  * Returns a reference to the Orbit Aspect Data with the specified @p index.
222  ***************************************************************************/
223 const GCOMOad& GCOMOads::at(const int& index) const
224 {
225  // Raise exception if index is out of range
226  if (index < 0 || index >= size()) {
227  throw GException::out_of_range(G_AT, "Orbit Aspect Data index",
228  index, size());
229  }
230 
231  // Return reference
232  return m_oads[index];
233 }
234 
235 
236 /***********************************************************************//**
237  * @brief Append Orbit Aspect Data to container
238  *
239  * @param[in] oad Orbit Aspect Data.
240  * @return Reference to appended Orbit Aspect Data.
241  *
242  * Appends Orbit Aspect Data to the container by making a deep copy of the
243  * Orbit Aspect Data.
244  ***************************************************************************/
246 {
247  // Append oad to list
248  m_oads.push_back(oad);
249 
250  // Return reference
251  return m_oads[size()-1];
252 }
253 
254 
255 /***********************************************************************//**
256  * @brief Insert Orbit Aspect Data into container
257  *
258  * @param[in] index Orbit Aspect Data index (0,...,size()-1).
259  * @param[in] oad Orbit Aspect Data.
260  *
261  * @exception GException::out_of_range
262  * Orbit Aspect Data index is out of range.
263  *
264  * Inserts an @p Orbit Aspect Data into the container before the Orbit
265  * Aspect Data with the specified @p index.
266  ***************************************************************************/
267 GCOMOad& GCOMOads::insert(const int& index, const GCOMOad& oad)
268 {
269  // Compile option: raise exception if index is out of range
270  #if defined(G_RANGE_CHECK)
271  if (is_empty()) {
272  if (index > 0) {
273  throw GException::out_of_range(G_INSERT, "Orbit Aspect Data index",
274  index, size());
275  }
276  }
277  else {
278  if (index < 0 || index >= size()) {
279  throw GException::out_of_range(G_INSERT, "Orbit Aspect Data index",
280  index, size());
281  }
282  }
283  #endif
284 
285  // Inserts Orbit Aspect Data
286  m_oads.insert(m_oads.begin()+index, oad);
287 
288  // Return reference
289  return m_oads[index];
290 }
291 
292 
293 /***********************************************************************//**
294  * @brief Remove Orbit Aspect Data from container
295  *
296  * @param[in] index Orbit Aspect Data index (0,...,size()-1).
297  *
298  * @exception GException::out_of_range
299  * Orbit Aspect Data index is out of range.
300  *
301  * Remove Orbit Aspect Data of specified @p index from container.
302  ***************************************************************************/
303 void GCOMOads::remove(const int& index)
304 {
305  // Compile option: raise exception if index is out of range
306  #if defined(G_RANGE_CHECK)
307  if (index < 0 || index >= size()) {
308  throw GException::out_of_range(G_REMOVE, "Orbit Aspect Data index",
309  index, size());
310  }
311  #endif
312 
313  // Erase Orbit Aspect Data from container
314  m_oads.erase(m_oads.begin() + index);
315 
316  // Return
317  return;
318 }
319 
320 
321 /***********************************************************************//**
322  * @brief Append Orbit Aspect Data container
323  *
324  * @param[in] oads COMPTEL Orbit Aspect Data container.
325  *
326  * Append COMPTEL Orbit Aspect Data container to the container.
327  ***************************************************************************/
328 void GCOMOads::extend(const GCOMOads& oads)
329 {
330  // Do nothing if COMPTEL Orbit Aspect Data container is empty
331  if (!oads.is_empty()) {
332 
333  // Get size. Note that we extract the size first to avoid an
334  // endless loop that arises when a container is appended to
335  // itself.
336  int num = oads.size();
337 
338  // Reserve enough space
339  reserve(size() + num);
340 
341  // Loop over all elements and append them to container
342  for (int i = 0; i < num; ++i) {
343  m_oads.push_back(oads[i]);
344  }
345 
346  } // endif: COMPTEL Orbit Aspect Data container was not empty
347 
348  // Return
349  return;
350 }
351 
352 
353 /***********************************************************************//**
354  * @brief Load COMPTEL Orbit Aspect Data FITS file
355  *
356  * @param[in] filename COMPTEL OAD FITS file name.
357  *
358  * Loads an COMPTEL Orbit Aspect Data FITS file in the container.
359  ***************************************************************************/
360 void GCOMOads::load(const GFilename& filename)
361 {
362  // Open FITS file
363  GFits fits(filename);
364 
365  // Get HDU (pointer is always valid)
366  const GFitsTable& hdu = *fits.table(1);
367 
368  // Read Orbit Aspect Data FITS table
369  read(hdu);
370 
371  // Close FITS file
372  fits.close();
373 
374  // Return
375  return;
376 }
377 
378 
379 /***********************************************************************//**
380  * @brief Read COMPTEL Orbit Aspect Data FITS table
381  *
382  * @param[in] table COMPTEL OAD FITS table.
383  *
384  * Reads COMPTEL Orbit Aspect Data FITS table into the container.
385  ***************************************************************************/
386 void GCOMOads::read(const GFitsTable& table)
387 {
388  // Earth radius in km including atmosphere (see EVP routine EVSCDJ Rev.7)
389  const double era = 6451.03;
390 
391  // Set Earth radius in km excluding atmosphere
392  const double erd = 6378.5;
393  const double erd_min = erd + 300.0; // Lower limit on orbit
394  const double erd_max = erd + 530.0; // Upper limit on orbit
395 
396  // Clear
397  clear();
398 
399  // Extract number of records in FITS file
400  int num = table.nrows();
401 
402  // If there are records then load them
403  if (num > 0) {
404 
405  // Reserve data
406  m_oads.reserve(num);
407 
408  // Get column pointers
409  const GFitsTableCol* ptr_tjd = table["TJD"]; // days
410  const GFitsTableCol* ptr_tics = table["TICS"]; // ticks
411  const GFitsTableCol* ptr_gcaz = table["GCAZ"]; // rad
412  const GFitsTableCol* ptr_gcel = table["GCEL"]; // rad
413  const GFitsTableCol* ptr_posx = table["POSX"]; // km
414  const GFitsTableCol* ptr_posy = table["POSY"]; // km
415  const GFitsTableCol* ptr_posz = table["POSZ"]; // km
416  const GFitsTableCol* ptr_zrasc = table["ZRASC"]; // rad
417  const GFitsTableCol* ptr_zdecl = table["ZDECL"]; // rad
418  const GFitsTableCol* ptr_xrasc = table["XRASC"]; // rad
419  const GFitsTableCol* ptr_xdecl = table["XDECL"]; // rad
420 
421  // Initialise Earth radius angle
422  double georad = 73.5;
423 
424  // Copy data from columns into records
425  for (int i = 0; i < num; ++i) {
426 
427  // Allocate OAD record
428  GCOMOad oad;
429 
430  // Get time of record
431  int tjd = ptr_tjd->integer(i);
432  int tics = ptr_tics->integer(i);
433 
434  // Store time information. The stop time is defined as the start
435  // time plus 131071 tics, since the length of one superpacket is
436  // 16.384 secs, i.e. 16.384 * 8000 = 131072 ticks
437  oad.tjd(tjd);
438  oad.tics(tics);
439  oad.tstart(gammalib::com_time(tjd, tics));
440  oad.tstop(gammalib::com_time(tjd, tics + 131071));
441 
442  // Set geocentre azimuth and zenith angle in deg
443  oad.gcaz(ptr_gcaz->real(i) * gammalib::rad2deg);
444  oad.gcel(ptr_gcel->real(i) * gammalib::rad2deg);
445 
446  // Set telescope Z- and X-axes
447  GSkyDir zaxis;
448  GSkyDir xaxis;
449  zaxis.radec(ptr_zrasc->real(i), ptr_zdecl->real(i));
450  xaxis.radec(ptr_xrasc->real(i), ptr_xdecl->real(i));
451  oad.zaxis(zaxis);
452  oad.xaxis(xaxis);
453 
454  // Compute apparent radius of Earth
455  double radius = std::sqrt(ptr_posx->real(i) * ptr_posx->real(i) +
456  ptr_posy->real(i) * ptr_posy->real(i) +
457  ptr_posz->real(i) * ptr_posz->real(i));
458  if ((radius > erd_min) && (radius < erd_max)) {
459  georad = std::asin(era/radius) * gammalib::rad2deg;
460  }
461  #if defined(G_GEORAD_WARNING)
462  else {
463  std::string msg = "Error in CGRO position. Distance from "
464  "geocentre is "+gammalib::str(radius)+ " km "
465  "while it should be in the interval ["+
466  gammalib::str(erd_min)+","+
467  gammalib::str(erd_max)+"] km. Use previous "
468  "spacecraft altitude.";
470  }
471  #endif
472  oad.georad(georad);
473 
474  // Append record
475  m_oads.push_back(oad);
476 
477  } // endfor: looped over OAD records
478 
479  } // endif: there were records to load
480 
481  // Return
482  return;
483 }
484 
485 
486 /***********************************************************************//**
487  * @brief Print COMPTEL Orbit Aspect Data container
488  *
489  * @param[in] chatter Chattiness.
490  * @return String containing COMPTEL Orbit Aspect Data container information.
491  ***************************************************************************/
492 std::string GCOMOads::print(const GChatter& chatter) const
493 {
494  // Initialise result string
495  std::string result;
496 
497  // Continue only if chatter is not silent
498  if (chatter != SILENT) {
499 
500  // Append header
501  result.append("=== GCOMOads ===");
502 
503  // Append container information
504  result.append("\n"+gammalib::parformat("Superpackets"));
505  result.append(gammalib::str(size()));
506  if (size() > 0) {
507 
508  // Append time range
509  result.append("\n"+gammalib::parformat("TJD range"));
510  result.append(gammalib::str(m_oads[0].tjd()));
511  result.append(":");
512  result.append(gammalib::str(m_oads[0].tics()));
513  result.append(" - ");
514  result.append(gammalib::str(m_oads[size()-1].tjd()));
515  result.append(":");
516  result.append(gammalib::str(m_oads[size()-1].tics()));
517  result.append("\n"+gammalib::parformat("MJD range"));
518  result.append(gammalib::str(m_oads[0].tstart().mjd()));
519  result.append(" - ");
520  result.append(gammalib::str(m_oads[size()-1].tstop().mjd()));
521  result.append(" days");
522  result.append("\n"+gammalib::parformat("UTC range"));
523  result.append(m_oads[0].tstart().utc());
524  result.append(" - ");
525  result.append(m_oads[size()-1].tstop().utc());
526 
527  // Append detailed information
528  GChatter reduced_chatter = gammalib::reduce(chatter);
529  if (reduced_chatter > SILENT) {
530 
531  // Append TJDs
532  int tjd = 0;
533  int num = 0;
534  for (int i = 0; i < size(); ++i) {
535  if (m_oads[i].tjd() != tjd) {
536  if (num > 0) {
537  std::string key = "TJD "+gammalib::str(tjd);
538  result.append("\n"+gammalib::parformat(key));
539  result.append(gammalib::str(num)+" superpackets");
540  }
541  tjd = m_oads[i].tjd();
542  num = 1;
543  }
544  else {
545  num++;
546  }
547  }
548  std::string key = "TJD "+gammalib::str(tjd);
549  result.append("\n"+gammalib::parformat(key));
550  result.append(gammalib::str(num)+" superpackets");
551 
552  } // endif: detailed information requested
553 
554  } // endif: there were records
555 
556  } // endif: chatter was not silent
557 
558  // Return result
559  return result;
560 }
561 
562 
563 /*==========================================================================
564  = =
565  = Private methods =
566  = =
567  ==========================================================================*/
568 
569 /***********************************************************************//**
570  * @brief Initialise class members
571  ***************************************************************************/
573 {
574  // Initialise members
575  m_oads.clear();
576 
577  // Return
578  return;
579 }
580 
581 
582 /***********************************************************************//**
583  * @brief Copy class members
584  *
585  * @param[in] oads COMPTEL Orbit Aspect Data container.
586  ***************************************************************************/
588 {
589  // Copy members
590  m_oads = oads.m_oads;
591 
592  // Return
593  return;
594 }
595 
596 
597 /***********************************************************************//**
598  * @brief Delete class members
599  ***************************************************************************/
601 {
602  // Return
603  return;
604 }
GCOMOads * clone(void) const
Clone COMPTEL Orbit Aspect Data container.
Definition: GCOMOads.cpp:184
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition: GFits.cpp:472
std::vector< GCOMOad > m_oads
Orbit Aspect Data records.
Definition: GCOMOads.hpp:89
void warning(const std::string &origin, const std::string &message)
Emits warning.
Definition: GTools.cpp:1265
GCOMOad & insert(const int &index, const GCOMOad &oad)
Insert Orbit Aspect Data into container.
Definition: GCOMOads.cpp:267
#define G_REMOVE
Definition: GCOMOads.cpp:41
const float & gcel(void) const
Return Geocentre zenith angle.
Definition: GCOMOad.hpp:276
bool is_empty(void) const
Signals if there are no Orbit Aspect Data in container.
Definition: GCOMOads.hpp:156
const GSkyDir & xaxis(void) const
Return telescope X-axis.
Definition: GCOMOad.hpp:364
std::string print(const GChatter &chatter=NORMAL) const
Print COMPTEL Orbit Aspect Data container.
Definition: GCOMOads.cpp:492
const GTime & tstop(void) const
Return stop time of superpacket.
Definition: GCOMOad.hpp:160
int size(void) const
Return number of Orbit Aspect Data in container.
Definition: GCOMOads.hpp:141
FITS file class.
Definition: GFits.hpp:63
void clear(void)
Clear COMPTEL Orbit Aspect Data container.
Definition: GCOMOads.cpp:166
FITS file class interface definition.
GCOMOad & at(const int &index)
Return reference to Orbit Aspect Data.
Definition: GCOMOads.cpp:200
#define G_INSERT
Definition: GCOMOads.cpp:40
GTime com_time(const int &tjd, const int &tics)
Convert TJD and COMPTEL ticks in GTime object.
const float & georad(void) const
Return apparent radius of Earth.
Definition: GCOMOad.hpp:305
COMPTEL Orbit Aspect Data class.
Definition: GCOMOad.hpp:48
Implementation of support function used by COMPTEL classes.
const GTime & tstart(void) const
Return start time of superpacket.
Definition: GCOMOad.hpp:129
void read(const GFitsTable &table)
Read COMPTEL Orbit Aspect Data FITS table.
Definition: GCOMOads.cpp:386
void radec(const double &ra, const double &dec)
Set equatorial sky direction (radians)
Definition: GSkyDir.cpp:175
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1268
void load(const GFilename &filename)
Load COMPTEL Orbit Aspect Data FITS file.
Definition: GCOMOads.cpp:360
#define G_AT
Definition: GCOMOads.cpp:39
const GSkyDir & zaxis(void) const
Return telescope Z-axis.
Definition: GCOMOad.hpp:334
Filename class.
Definition: GFilename.hpp:62
Abstract interface for FITS table column.
virtual ~GCOMOads(void)
Destructor.
Definition: GCOMOads.cpp:114
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
GChatter
Definition: GTypemaps.hpp:33
GCOMOad & append(const GCOMOad &oad)
Append Orbit Aspect Data to container.
Definition: GCOMOads.cpp:245
void extend(const GCOMOads &oads)
Append Orbit Aspect Data container.
Definition: GCOMOads.cpp:328
GVector asin(const GVector &vector)
Computes arcsin of vector elements.
Definition: GVector.cpp:1016
const int & nrows(void) const
Return number of rows in table.
Definition: GFitsTable.hpp:119
void init_members(void)
Initialise class members.
Definition: GCOMOads.cpp:572
void free_members(void)
Delete class members.
Definition: GCOMOads.cpp:600
#define G_READ
Definition: GCOMOads.cpp:42
virtual int integer(const int &row, const int &inx=0) const =0
virtual double real(const int &row, const int &inx=0) const =0
const int & tjd(void) const
Return Truncated Julian Days of Orbit Aspect Record.
Definition: GCOMOad.hpp:189
GCOMOads & operator=(const GCOMOads &oads)
Assignment operator.
Definition: GCOMOads.cpp:136
void reserve(const int &num)
Reserves space for Orbit Aspect Data in container.
Definition: GCOMOads.hpp:170
GChatter reduce(const GChatter &chatter)
Reduce chattiness by one level.
Definition: GTypemaps.hpp:65
Exception handler interface definition.
COMPTEL Orbit Aspect Data container class.
Definition: GCOMOads.hpp:51
GCOMOads(void)
Void constructor.
Definition: GCOMOads.cpp:63
void remove(const int &index)
Remove Orbit Aspect Data from container.
Definition: GCOMOads.cpp:303
const double rad2deg
Definition: GMath.hpp:44
const int & tics(void) const
Return tics of Orbit Aspect Record.
Definition: GCOMOad.hpp:218
void copy_members(const GCOMOads &oads)
Copy class members.
Definition: GCOMOads.cpp:587
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
Sky direction class.
Definition: GSkyDir.hpp:62
void close(void)
Close FITS file.
Definition: GFits.cpp:1314
COMPTEL Orbit Aspect Data container class definition.
Mathematical function definitions.
const float & gcaz(void) const
Return Geocentre azimuth angle.
Definition: GCOMOad.hpp:247
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413
FITS table abstract base class interface definition.