GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCOMSupport.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCOMSupport.cpp - COMPTEL support functions *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2012-2019 by Juergen Knoedlseder *
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 GCOMSupport.hpp
23  * @brief Implementation of support function used by COMPTEL classes
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GTools.hpp"
32 #include "GMath.hpp"
33 #include "GWcs.hpp"
34 #include "GWcsCAR.hpp"
35 #include "GTime.hpp"
36 #include "GSkyMap.hpp"
37 #include "GCOMSupport.hpp"
38 
39 /* __ Coding definitions _________________________________________________ */
40 
41 /* __ Debug definitions __________________________________________________ */
42 
43 
44 /***********************************************************************//**
45  * @brief Changes Mercator's projection to cartesian projection
46  *
47  * @param[in,out] map Skymap.
48  *
49  * Changes the World Coordinate System of a sky map from Mercartor's
50  * projection to a cartesian projection. This transformation is needed to
51  * correct the wrong WCS headers that are found in the COMPTEL files
52  * distributed through the HEASARC web site.
53  *
54  * The method operates only on sky maps that use Mercator's projection.
55  * Nothing is done for all other sky map. This method can thus be
56  * transparently be applied to the sky maps read from COMPTEL data files.
57  ***************************************************************************/
59 {
60  // Get WCS poiunter
61  const GWcs* wcs = dynamic_cast<const GWcs*>(map.projection());
62 
63  // Apply kluge only if the skymap is in Mercator's projection
64  if (wcs != NULL && wcs->code() == "MER") {
65 
66  // Extract original WCS definition
67  double crval1 = wcs->crval(0);
68  double crval2 = wcs->crval(1);
69  double crpix1 = wcs->crpix(0);
70  double crpix2 = wcs->crpix(1);
71  double cdelt1 = wcs->cdelt(0);
72  double cdelt2 = wcs->cdelt(1);
73  int naxis1 = map.nx();
74  int naxis2 = map.ny();
75 
76  // Compute adjustment of WCS definition to map centre
77  double adjust1 = (double(naxis1) + 1.0) / 2.0 - crpix1;
78  double adjust2 = (double(naxis2) + 1.0) / 2.0 - crpix2;
79 
80  // Adjust WCS definition to map centre
81  crval1 += adjust1;
82  crpix1 += adjust1;
83  crval2 += adjust2;
84  crpix2 += adjust2;
85 
86  // Allocate WCS
87  GWcsCAR car(wcs->coordsys(), crval1, crval2,
88  crpix1, crpix2,
89  cdelt1, cdelt2);
90 
91  // Set new WCS
92  map.projection(car);
93 
94  } // endif: skymap was in Mercator projection
95 
96  // Return
97  return;
98 }
99 
100 
101 /***********************************************************************//**
102  * @brief Return D1 energy deposit
103  *
104  * @param[in] energy Input photon energy (MeV).
105  * @param[in] phigeo Geometrical scatter angle (deg).
106  * @return D1 energy deposit in MeV.
107  ***************************************************************************/
108 double gammalib::com_energy1(const double& energy, const double& phigeo)
109 {
110  // Compute D2 energy deposit
111  double e2 = com_energy2(energy, phigeo);
112 
113  // Return D1 energy deposit
114  return (energy - e2);
115 }
116 
117 
118 /***********************************************************************//**
119  * @brief Return D2 energy deposit
120  *
121  * @param[in] energy Input photon energy (MeV).
122  * @param[in] phigeo Geometrical scatter angle (deg).
123  * @return D2 energy deposit in MeV.
124  ***************************************************************************/
125 double gammalib::com_energy2(const double& energy, const double& phigeo)
126 {
127  // Compute 1-cos(phigeo)
128  double one_minus_cos = 1.0 - std::cos(phigeo * gammalib::deg2rad);
129 
130  // Compute D2 energy deposit
131  double e2 = energy / (one_minus_cos * energy / gammalib::mec2 + 1.0);
132 
133  // Return D2 energy deposit
134  return e2;
135 }
136 
137 
138 /***********************************************************************//**
139  * @brief Convert TJD and COMPTEL ticks in GTime object
140  *
141  * @param[in] tjd Truncated Julian Days (days).
142  * @param[in] tics COMPTEL ticks (1/8 ms).
143  * @return Time.
144  *
145  * Converts TJD and COMPTEL ticks into a GTime object. COMPTEL times are
146  * given in UTC, i.e. 8393:0 converts into 1991-05-17T00:00:00 UT
147  * (see COM-RP-UNH-DRG-037).
148  ***************************************************************************/
149 GTime gammalib::com_time(const int& tjd, const int& tics)
150 {
151  // Compute MJD
152  double mjd = double(tjd) + 40000.0;
153 
154  // Set time and retrieve result in native seconds
155  GTime time;
156  time.mjd(mjd, "UTC");
157  double secs = time.secs();
158 
159  // Add ticks and set time in native seconds
160  secs += double(tics) * 0.000125;
161 
162  // Set time
163  time.secs(secs);
164 
165  // Return time
166  return time;
167 }
168 
169 
170 /***********************************************************************//**
171  * @brief Convert GTime in COMPTEL TJD
172  *
173  * @param[in] time Time.
174  * @return Truncated Julian Days (days).
175  *
176  * Converts GTime object in TJD.
177  ***************************************************************************/
178 int gammalib::com_tjd(const GTime& time)
179 {
180  // Compute TJD
181  int tjd = int(time.mjd("UTC") - 40000.0);
182 
183  // Return TJD
184  return tjd;
185 }
186 
187 
188 /***********************************************************************//**
189  * @brief Convert GTime in COMPTEL tics
190  *
191  * @param[in] time Time.
192  * @return COMPTEL ticks (1/8 ms).
193  *
194  * Converts GTime object in COMPTEL ticks (1/8 ms).
195  ***************************************************************************/
196 int gammalib::com_tics(const GTime& time)
197 {
198  // Compute COMPTEL time at 0 tics
199  GTime tjd = com_time(com_tjd(time), 0);
200 
201  // Compute time difference in seconds
202  int tics = int((time - tjd) * 8000.0 + 0.5); // rounding to nearest int
203 
204  // Return tics
205  return tics;
206 }
const double mec2
Definition: GTools.hpp:52
Sky map class.
Definition: GSkyMap.hpp:89
const int & ny(void) const
Returns number of pixels in y coordinate.
Definition: GSkyMap.hpp:366
const GSkyProjection * projection(void) const
Returns pointer to sky projection.
Definition: GSkyMap.hpp:414
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1100
double com_energy2(const double &energy, const double &phigeo)
Return D2 energy deposit.
double crval(const int &inx) const
Return value of reference pixel.
Definition: GWcs.cpp:877
Time class.
Definition: GTime.hpp:54
Gammalib tools definition.
GTime com_time(const int &tjd, const int &tics)
Convert TJD and COMPTEL ticks in GTime object.
Sky map class definition.
virtual std::string code(void) const =0
Implementation of support function used by COMPTEL classes.
const double deg2rad
Definition: GMath.hpp:43
double crpix(const int &inx) const
Return reference pixel.
Definition: GWcs.cpp:900
void com_wcs_mer2car(GSkyMap &map)
Changes Mercator&#39;s projection to cartesian projection.
Definition: GCOMSupport.cpp:58
int com_tjd(const GTime &time)
Convert GTime in COMPTEL TJD.
double cdelt(const int &inx) const
Return pixel size.
Definition: GWcs.cpp:923
Plate carree (CAR) projection class definition.
virtual std::string coordsys(void) const
Returns coordinate system.
Plate carree (CAR) projection class definition.
Definition: GWcsCAR.hpp:43
const double & secs(void) const
Return time in seconds in native reference (TT)
Definition: GTime.hpp:153
Abstract world coordinate system base class.
Definition: GWcs.hpp:51
int com_tics(const GTime &time)
Convert GTime in COMPTEL tics.
const int & nx(void) const
Returns number of pixels in x coordinate.
Definition: GSkyMap.hpp:350
double mjd(void) const
Return time in Modified Julian Days (TT)
Definition: GTime.cpp:319
double com_energy1(const double &energy, const double &phigeo)
Return D1 energy deposit.
Time class interface definition.
Abstract world coordinate system base class definition.
Mathematical function definitions.