Jürgen Knödlseder / gammalib

Commit 445fd8275cc806376980661ecd70f288ea05c623

Authored by Jürgen Knödlseder
1 parent 2fa08b9f

Add GEphemerides class (#3973)

Inline Side-by-side
Showing 12 changed files with 431 additions and 15 deletions
ChangeLog
  View file @445fd82
1   -2022-01-09
  1 +2022-01-14
2 2  
3 3 * Version 2.0.0 released
4 4 ========================
... ...
Makefile.am
  View file @445fd82
1 1 #############################################################################
2 2 # Makefile.am for GammaLib package #
3 3 # ------------------------------------------------------------------------- #
4   -# Copyright (C) 2016-2021 Juergen Knoedlseder #
  4 +# Copyright (C) 2016-2022 Juergen Knoedlseder #
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # This program is free software: you can redistribute it and/or modify #
... ... @@ -28,7 +28,7 @@ INST_DIR = inst
28 28 endif
29 29  
30 30 # Subdirectories
31   -SUBDIRS = include $(INST_DIR) src pyext examples test
  31 +SUBDIRS = include $(INST_DIR) src pyext refdata examples test
32 32  
33 33 # Pkgconfig directory
34 34 pkgconfigdir = $(libdir)/pkgconfig
... ...
1 1 New Features and Important Changes in GammaLib 2.0.0
2 2  
3   -9 January 2022
  3 +14 January 2022
4 4  
5 5  
6 6 1. Introduction
... ... @@ -17,6 +17,7 @@ The following classes have been added:
17 17 - GResponseVectorCache
18 18 - GPulsar
19 19 - GPulsarEphemeris
  20 +- GEphemerides
20 21 - GSkyDirs
21 22 - GSkyRegionRect
22 23 - GModelAssociation
... ... @@ -277,6 +278,10 @@ of functions (#3328).
277 278  
278 279 11. Observation module
279 280 ----------------------
  281 +Add GEphemerides class to deal with JPL ephemerides and specifically the
  282 +DE200 ephemerides. A JPL DE200 ephemerides file was added to a refdata
  283 +direction, the file is installed in $GAMMALIB/share/refdata (#3973).
  284 +
280 285 Add GPhases::length() method that returns the total length of the phase
281 286 intervals (#3973).
282 287  
... ...
README.md
  View file @445fd82
1 1 GammaLib information
2 2 ====================
3   -* Version: 2.0.0.dev (9 January 2022)
  3 +* Version: 2.0.0.dev (14 January 2022)
4 4  
5 5 [![Build Status](https://cta-jenkins.irap.omp.eu/buildStatus/icon?job=gammalib-integrate-os)](https://cta-jenkins.irap.omp.eu/job/gammalib-integrate-os/)
6 6  
... ...
configure.ac
  View file @445fd82
... ... @@ -2,7 +2,7 @@
2 2 # Configuration for GammaLib - a versatile toolbox for high-level analysis #
3 3 # of astronomical gamma-ray data #
4 4 # ------------------------------------------------------------------------- #
5   -# Copyright (C) 2006-2021 Juergen Knoedlseder #
  5 +# Copyright (C) 2006-2022 Juergen Knoedlseder #
6 6 # ------------------------------------------------------------------------- #
7 7 # #
8 8 # This program is free software: you can redistribute it and/or modify #
... ... @@ -1485,6 +1485,7 @@ AC_CONFIG_FILES([Makefile
1485 1485 pyext/setup.py
1486 1486 pyext/gammalib/__init__.py
1487 1487 dev/conda.recipe/meta.yaml
  1488 + refdata/Makefile
1488 1489 examples/Makefile
1489 1490 examples/cpp/Makefile
1490 1491 examples/python/Makefile
... ...
include/GEphemerides.hpp
  View file @445fd82
... ... @@ -29,10 +29,14 @@
29 29  
30 30 /* __ Includes ___________________________________________________________ */
31 31 #include <string>
  32 +#include <vector>
32 33 #include "GBase.hpp"
33 34 #include "GFilename.hpp"
  35 +#include "GTime.hpp"
  36 +#include "GVector.hpp"
34 37  
35 38 /* __ Forward declarations _______________________________________________ */
  39 +class GSkyDir;
36 40  
37 41 /* __ Constants __________________________________________________________ */
38 42  
... ... @@ -62,17 +66,38 @@ public:
62 66 virtual std::string print(const GChatter& chatter = NORMAL) const;
63 67  
64 68 // Other methods
65   - // TODO: Add any further methods that are needed
  69 + int size(void) const;
  70 + bool is_empty(void) const;
  71 + const std::string& name(void) const;
  72 + void name(const std::string& name);
  73 + void load(const GFilename& filename);
  74 + void ephemeris(const GTime& time,
  75 + GVector* rce,
  76 + GVector* rcs,
  77 + GVector* vce,
  78 + double* etut) const;
  79 + double geo2ssb(const GTime& time,
  80 + const GSkyDir& srcdir) const;
66 81  
67 82 protected:
68 83 // Protected methods
69 84 void init_members(void);
70 85 void copy_members(const GEphemerides& ephemerides);
71 86 void free_members(void);
  87 + void fetch_data(void);
72 88  
73 89 // Protected members
74   - std::string m_name; //!< Ephemerides (e.g. DE200)
75   - GFilename m_filename; //!< Ephemerides filename
  90 + std::string m_name; //!< Ephemerides (e.g. DE200)
  91 + GFilename m_filename; //!< Ephemerides filename
  92 + GTime m_tstart; //!< Ephemerides validity start time
  93 + GTime m_tstop; //!< Ephemerides validity stop time
  94 + std::vector<GTime> m_times; //!< Times of vectors
  95 + std::vector<GVector> m_earth; //!< Earth vectors
  96 + std::vector<GVector> m_earth_dt; //!< First derivative of Earth vectors
  97 + std::vector<GVector> m_earth_d2t; //!< Second derivative of Earth vectors
  98 + std::vector<GVector> m_earth_d3t; //!< Third derivative of Earth vectors
  99 + std::vector<GVector> m_sun; //!< Sun vectors
  100 + std::vector<double> m_tdb2tt; //!< TBD to TT conversion term
76 101 };
77 102  
78 103  
... ... @@ -87,4 +112,57 @@ std::string GEphemerides::classname(void) const
87 112 return ("GEphemerides");
88 113 }
89 114  
  115 +
  116 +/***********************************************************************//**
  117 + * @brief Return number of ephemerides
  118 + *
  119 + * @return Number of ephemerides.
  120 + *
  121 + * Returns the number of ephemerides.
  122 + ***************************************************************************/
  123 +inline
  124 +int GEphemerides::size(void) const
  125 +{
  126 + return (int)m_times.size();
  127 +}
  128 +
  129 +
  130 +/***********************************************************************//**
  131 + * @brief Signals if there are no ephemerides
  132 + *
  133 + * @return True if there are no ephemerides, false otherwise.
  134 + *
  135 + * Signals if there are no ephemerides.
  136 + ***************************************************************************/
  137 +inline
  138 +bool GEphemerides::is_empty(void) const
  139 +{
  140 + return (m_times.empty());
  141 +}
  142 +
  143 +
  144 +/***********************************************************************//**
  145 + * @brief Return ephemerides name
  146 + *
  147 + * @return Ephemerides name
  148 + ***************************************************************************/
  149 +inline
  150 +const std::string& GEphemerides::name(void) const
  151 +{
  152 + return (m_name);
  153 +}
  154 +
  155 +
  156 +/***********************************************************************//**
  157 + * @brief Set ephemerides name
  158 + *
  159 + * @param[in] name Ephemerides name.
  160 + ***************************************************************************/
  161 +inline
  162 +void GEphemerides::name(const std::string& name)
  163 +{
  164 + m_name = name;
  165 + return;
  166 +}
  167 +
90 168 #endif /* GEPHEMERIDES_HPP */
... ...
include/Makefile.am
  View file @445fd82
... ... @@ -21,6 +21,7 @@
21 21  
22 22 # Define package included headers (to be installed)
23 23 pkginclude_HEADERS = GammaLib.hpp \
  24 + GEphemerides.hpp \
24 25 GTypemaps.hpp \
25 26 GBase.hpp \
26 27 GContainer.hpp \
... ...
pyext/GEphemerides.i
  View file @445fd82
... ... @@ -48,7 +48,18 @@ public:
48 48 virtual std::string classname(void) const;
49 49  
50 50 // Other methods
51   - // TODO: Copy methods from GEphemerides.hpp file
  51 + int size(void) const;
  52 + bool is_empty(void) const;
  53 + const std::string& name(void) const;
  54 + void name(const std::string& name);
  55 + void load(const GFilename& filename);
  56 + void ephemeris(const GTime& time,
  57 + GVector* rce,
  58 + GVector* rcs,
  59 + GVector* vce,
  60 + double* etut) const;
  61 + double geo2ssb(const GTime& time,
  62 + const GSkyDir& srcdir) const;
52 63 };
53 64  
54 65  
... ...
refdata/Makefile.am 0 → 100644
  View file @445fd82
  1 +#############################################################################
  2 +# Makefile for reference data #
  3 +# ------------------------------------------------------------------------- #
  4 +# Copyright (C) 2022 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 +# Models to install
  23 +refdata = ephem_jpl_de200.fits
  24 +
  25 +# Install models
  26 +refdatadir = $(datadir)/refdata
  27 +refdata_DATA = $(refdata)
  28 +
  29 +# Extra files to be distributed in the tar file
  30 +EXTRA_DIST = $(refdata)
  31 +
... ...
refdata/ephem_jpl_de200.fits 0 → 100644
View file @445fd82
No preview for this file type
src/obs/GEphemerides.cpp
  View file @445fd82
... ... @@ -29,8 +29,14 @@
29 29 #include <config.h>
30 30 #endif
31 31 #include "GEphemerides.hpp"
  32 +#include "GFits.hpp"
  33 +#include "GFitsTable.hpp"
  34 +#include "GFitsTableCol.hpp"
  35 +#include "GSkyDir.hpp"
32 36  
33 37 /* __ Method name definitions ____________________________________________ */
  38 +#define G_EPHEMERIS "GEphemerides::ephemeris(GTime&, GVector*, GVector*, "\
  39 + "GVector*, double*)"
34 40  
35 41 /* __ Macros _____________________________________________________________ */
36 42  
... ... @@ -157,6 +163,244 @@ GEphemerides* GEphemerides::clone(void) const
157 163  
158 164  
159 165 /***********************************************************************//**
  166 + * @brief Load Ephemerides
  167 + *
  168 + * @param[in] filename Ephemerides file name.
  169 + *
  170 + * Load ephemerides from FITS file.
  171 + ***************************************************************************/
  172 +void GEphemerides::load(const GFilename& filename)
  173 +{
  174 + // Clear members
  175 + clear();
  176 +
  177 + // Open JPL ephemerides FITS file
  178 + GFits fits(filename);
  179 +
  180 + // Get ephemerides table
  181 + const GFitsTable* ephem = fits.table("EPHEM");
  182 +
  183 + // Extract number of events in FITS file
  184 + int num = ephem->nrows();
  185 +
  186 + // Continue if there are ephemerides
  187 + if (num > 0) {
  188 +
  189 + // Reserve data
  190 + m_times.reserve(num);
  191 + m_earth.reserve(num);
  192 + m_earth_dt.reserve(num);
  193 + m_earth_d2t.reserve(num);
  194 + m_earth_d3t.reserve(num);
  195 + m_sun.reserve(num);
  196 + m_tdb2tt.reserve(num);
  197 +
  198 + // Get validity interval
  199 + double tstart = ephem->real("TSTART");
  200 + m_tstart.jd(tstart);
  201 + m_tstop.jd(ephem->real("TSTOP"));
  202 +
  203 + // Get column pointers
  204 + const GFitsTableCol* ptr_earth = (*ephem)["EARTH"]; // light seconds
  205 + const GFitsTableCol* ptr_sun = (*ephem)["SUN"]; // light seconds
  206 + const GFitsTableCol* ptr_tbd2tt = (*ephem)["TIMEDIFF"]; // sec
  207 +
  208 + // Extract data
  209 + for (int i = 0; i < num; ++i) {
  210 +
  211 + // Set time
  212 + GTime time;
  213 + time.jd(tstart + double(i));
  214 +
  215 + // Set vectors
  216 + GVector earth(ptr_earth->real(i,0),
  217 + ptr_earth->real(i,1),
  218 + ptr_earth->real(i,2));
  219 + GVector earth_dt(ptr_earth->real(i,3),
  220 + ptr_earth->real(i,4),
  221 + ptr_earth->real(i,5));
  222 + GVector earth_d2t(ptr_earth->real(i,6),
  223 + ptr_earth->real(i,7),
  224 + ptr_earth->real(i,8));
  225 + GVector earth_d3t(ptr_earth->real(i,9),
  226 + ptr_earth->real(i,10),
  227 + ptr_earth->real(i,11));
  228 + GVector sun(ptr_sun->real(i,0),
  229 + ptr_sun->real(i,1),
  230 + ptr_sun->real(i,2));
  231 +
  232 + // Set tbd2tt
  233 + double tbd2tt = ptr_tbd2tt->real(i);
  234 +
  235 + // Push back data
  236 + m_times.push_back(time);
  237 + m_earth.push_back(earth);
  238 + m_earth_dt.push_back(earth_dt);
  239 + m_earth_d2t.push_back(earth_d2t);
  240 + m_earth_d3t.push_back(earth_d3t);
  241 + m_sun.push_back(sun);
  242 + m_tdb2tt.push_back(tbd2tt);
  243 +
  244 + } // endfor: looped over rows
  245 +
  246 + // Store attributes
  247 + m_name = ephem->string("NAME");
  248 + m_filename = filename;
  249 +
  250 + } // endif: there were ephemerides to load
  251 +
  252 + // Return
  253 + return;
  254 +}
  255 +
  256 +
  257 +/***********************************************************************//**
  258 + * @brief Get ephemeris vector and TBD->TT value for a given time
  259 + *
  260 + * @param[in] time Time.
  261 + * @param[out] rce Pointer to vector from SSBC to Earth (light-s).
  262 + * @param[out] rcs Pointer to vector from SSBC to Sun (light-s).
  263 + * @param[out] vce Pointer to first time derivative of vector from SSBC to
  264 + * Earth (light-s/s).
  265 + * @param[out] etut Time difference TDB-TT (s)
  266 + *
  267 + * Get ephemeris vector and TBD->TT value for a given time. Information is
  268 + * only returned for pointers that are not NULL.
  269 + *
  270 + * @todo Check whether the conversion to UTC is needed and whether utc_to_tdt
  271 + * needs to be added to etut since times are normally in TT, not UTC.
  272 + ***************************************************************************/
  273 +void GEphemerides::ephemeris(const GTime& time,
  274 + GVector* rce,
  275 + GVector* rcs,
  276 + GVector* vce,
  277 + double* etut) const
  278 +{
  279 + // Fetch ephemerides
  280 + const_cast<GEphemerides*>(this)->fetch_data();
  281 +
  282 + // Compute jd_utc as int and dt as day fraction, between -0.5 and 0.5
  283 + double jd = time.jd("UTC");
  284 + int jd_utc = int(jd);
  285 + double dt = jd - double(jd_utc);
  286 + if (dt > 0.5) {
  287 + jd_utc += 1;
  288 + dt -= 1.0;
  289 + }
  290 +
  291 + // Get Julian Day validity interval as int
  292 + int jd0 = int(m_tstart.jd());
  293 + int jd1 = int(m_tstop.jd());
  294 +
  295 + // If date is outside range of ephemeris then throw an exception
  296 + if (jd_utc < jd0) {
  297 + std::string msg = "Time JD "+gammalib::str(jd_utc)+ " is before "
  298 + "validity start JD "+gammalib::str(jd0)+" of "
  299 + "ephemerides. Please specify a time later than "
  300 + "the validity start.";
  301 + throw GException::invalid_argument(G_EPHEMERIS, msg);
  302 + }
  303 + else if (jd_utc > jd1) {
  304 + std::string msg = "Time JD "+gammalib::str(jd_utc)+ " is after "
  305 + "validity end JD "+gammalib::str(jd1)+" of "
  306 + "ephemerides. Please specify a time earlier than "
  307 + "the validity end.";
  308 + throw GException::invalid_argument(G_EPHEMERIS, msg);
  309 + }
  310 +
  311 + // Compute ephemeris index
  312 + int index = jd_utc - jd0;
  313 +
  314 + // Make rough computation of time derivative of TDB-TDT
  315 + double tbd2tt_dot = (index+1 < size()) ? m_tdb2tt[index+1] - m_tdb2tt[index]
  316 + : m_tdb2tt[index] - m_tdb2tt[index-1];
  317 +
  318 + // Get Earth position and velocity vector at the desired time using
  319 + // Taylor expansion
  320 + const double c1 = 1.0/2.0;
  321 + const double c2 = 1.0/6.0;
  322 + if (rce != NULL) {
  323 + *rce = m_earth[index] + (m_earth_dt[index] + (m_earth_d2t[index] *
  324 + c1 + m_earth_d3t[index] * dt * c2) * dt) * dt;
  325 + }
  326 + if (vce != NULL) {
  327 + *vce = (m_earth_dt[index] + (m_earth_d2t[index] + m_earth_d3t[index] *
  328 + dt * c1)) * dt * gammalib::sec2day;
  329 + }
  330 +
  331 + // Get Sun vector
  332 + if (rce != NULL) {
  333 + *rce = m_sun[index];
  334 + }
  335 +
  336 + // Get TBD-TT (seconds)
  337 + if (etut != NULL) {
  338 + double utc_to_tdt = time.secs("TT") - time.secs("UTC");
  339 + *etut = m_tdb2tt[index] + dt * tbd2tt_dot + utc_to_tdt;
  340 + }
  341 +
  342 + // Return
  343 + return;
  344 +}
  345 +
  346 +
  347 +/***********************************************************************//**
  348 + * @brief Get time difference between geocentric and SSB (seconds)
  349 + *
  350 + * @param[in] time Time.
  351 + * @param[in] srcdir Source direction.
  352 + * @return Time difference in seconds.
  353 + ***************************************************************************/
  354 +double GEphemerides::geo2ssb(const GTime& time, const GSkyDir& srcdir) const
  355 +{
  356 + // Set constants
  357 + const double sunrad = 2.315;
  358 + const double aultsc = 499.004782; // What ???
  359 + const double gauss = 0.01720209895; // What ???
  360 + const double schwarzschild_radius = (gauss*gauss) *
  361 + (aultsc*aultsc*aultsc) *
  362 + (gammalib::sec2day*gammalib::sec2day);
  363 +
  364 + // Get ephemerides
  365 + GVector rce(3);
  366 + GVector rcs(3);
  367 + GVector vce(3);
  368 + double etut;
  369 + ephemeris(time, &rce, &rcs, &vce, &etut);
  370 +
  371 + // Get sky direction as celestial vector
  372 + GVector dir = srcdir.celvector();
  373 +
  374 + // Calculate light-travel time to barycenter in Euclidean space
  375 + double light_travel_time = dir * rce;
  376 +
  377 + // Now calculate the time delay due to the gravitational field of the
  378 + // Sun (I.I. Shapiro, Phys. Rev. Lett. 13, 789 (1964))
  379 + GVector rsa = rce - rcs;
  380 + double sundis = norm(rsa);
  381 + double sunsiz = sunrad/sundis;
  382 + double cos_theta = (dir * rsa) / sundis;
  383 +
  384 + // Special handling if sky direction is inside the Sun
  385 + /*
  386 + if (cth + 1.0 < 0.5 * sunsiz * sunsiz) {
  387 +
  388 + }
  389 + */
  390 +
  391 + // Compute Shapiro delay
  392 + double shapiro_delay = -2.0 * schwarzschild_radius * std::log(1.0 + cos_theta);
  393 +
  394 + // Compute time difference between geocentric and Solar System Barycentric
  395 + // frame
  396 + double geo2ssb = light_travel_time - shapiro_delay + etut;
  397 +
  398 + // Return barycentric correction (seconds)
  399 + return geo2ssb;
  400 +}
  401 +
  402 +
  403 +/***********************************************************************//**
160 404 * @brief Print Ephemerides
161 405 *
162 406 * @param[in] chatter Chattiness.
... ... @@ -178,6 +422,16 @@ std::string GEphemerides::print(const GChatter&amp; chatter) const
178 422 result.append(m_name);
179 423 result.append("\n"+gammalib::parformat("File name"));
180 424 result.append(m_filename.url());
  425 + if (!is_empty()) {
  426 + result.append("\n"+gammalib::parformat("Validity MJD range"));
  427 + result.append(gammalib::str(m_tstart.mjd()));
  428 + result.append(" - ");
  429 + result.append(gammalib::str(m_tstop.mjd()));
  430 + result.append("\n"+gammalib::parformat("Validity UTC range"));
  431 + result.append(m_tstart.utc());
  432 + result.append(" - ");
  433 + result.append(m_tstop.utc());
  434 + }
181 435  
182 436 } // endif: chatter was not silent
183 437  
... ... @@ -200,7 +454,16 @@ void GEphemerides::init_members(void)
200 454 // Initialise members
201 455 m_name.clear();
202 456 m_filename.clear();
203   -
  457 + m_tstart.clear();
  458 + m_tstop.clear();
  459 + m_times.clear();
  460 + m_earth.clear();
  461 + m_earth_dt.clear();
  462 + m_earth_d2t.clear();
  463 + m_earth_d3t.clear();
  464 + m_sun.clear();
  465 + m_tdb2tt.clear();
  466 +
204 467 // Return
205 468 return;
206 469 }
... ... @@ -214,8 +477,17 @@ void GEphemerides::init_members(void)
214 477 void GEphemerides::copy_members(const GEphemerides& ephemerides)
215 478 {
216 479 // Copy members
217   - m_name = ephemerides.m_name;
218   - m_filename = ephemerides.m_filename;
  480 + m_name = ephemerides.m_name;
  481 + m_filename = ephemerides.m_filename;
  482 + m_tstart = ephemerides.m_tstart;
  483 + m_tstop = ephemerides.m_tstop;
  484 + m_times = ephemerides.m_times;
  485 + m_earth = ephemerides.m_earth;
  486 + m_earth_dt = ephemerides.m_earth_dt;
  487 + m_earth_d2t = ephemerides.m_earth_d2t;
  488 + m_earth_d3t = ephemerides.m_earth_d3t;
  489 + m_sun = ephemerides.m_sun;
  490 + m_tdb2tt = ephemerides.m_tdb2tt;
219 491  
220 492 // Return
221 493 return;
... ... @@ -230,3 +502,21 @@ void GEphemerides::free_members(void)
230 502 // Return
231 503 return;
232 504 }
  505 +
  506 +
  507 +/***********************************************************************//**
  508 + * @brief Fetch ephemerides data
  509 + *
  510 + * This method loads the DE200 ephemerides data from the file provided in
  511 + * the reference data repository if no ephemerides are loaded.
  512 + ***************************************************************************/
  513 +void GEphemerides::fetch_data(void)
  514 +{
  515 + // If no ephemerides data are loaded then load the JPL DE200 ephemerides
  516 + if (is_empty()) {
  517 + load("$GAMMALIB/share/refdata/ephem_jpl_de200.fits");
  518 + }
  519 +
  520 + // Return
  521 + return;
  522 +}
... ...
src/obs/Makefile.am
  View file @445fd82
... ... @@ -54,8 +54,7 @@ sources = GEnergies.cpp \
54 54 GSource.cpp \
55 55 GPulsar.cpp \
56 56 GPulsarEphemeris.cpp \
57   - GEphemerides.cpp \
58   -
  57 + GEphemerides.cpp
59 58  
60 59 # Build libtool library
61 60 noinst_LTLIBRARIES = libobs.la
... ...