Jürgen Knödlseder / gammalib

Commit c16adc29e79dbd2a735d640bcb8f939702ca9987

Authored by Jürgen Knödlseder
1 parent fb633b7c

Support Galactic coordinates in GModelSpatialPointSource (#3988)

Inline Side-by-side
Showing 8 changed files with 304 additions and 208 deletions
ChangeLog
  View file @c16adc2
1   -2022-01-19
  1 +2022-02-11
2 2  
3 3 * Version 2.0.0 released
4 4 ========================
5 5  
  6 + Support Galactic coordinates in spatial models (#3988)
6 7 Implement CGRO clock correction (#3973)
7 8 Add support for COMPTEL pulsar analysis (#3973)
8 9 Add GCOMObservation::tim() and GCOMObservation::oads() setter methods
... ...
1 1 New Features and Important Changes in GammaLib 2.0.0
2 2  
3   -19 January 2022
  3 +11 February 2022
4 4  
5 5  
6 6 1. Introduction
... ... @@ -78,6 +78,7 @@ The following methods have been added:
78 78 - GModelSpatialDiffuseCube::flux(GSkyRegion&, GEnergy&, GTime&)
79 79 - GModelSpatialDiffuseMap::flux(GSkyRegion&, GEnergy&, GTime&)
80 80 - GModelSpatialPointSource::flux(GSkyRegion&, GEnergy&, GTime&)
  81 +- GModelSpatialPointSource::coordsys()
81 82 - GModelSpatialRadial::is_energy_dependent()
82 83 - GModelSpatialRadial::is_time_dependent()
83 84 - GModelSpectralFunc(GModelSpectral&, GEnergies&)
... ... @@ -148,6 +149,8 @@ The following methods have been added:
148 149  
149 150 The following methods have been removed:
150 151 - GFitsTableCol::scale()
  152 +- GModelSpatialPointSource::ra()
  153 +- GModelSpatialPointSource::dec()
151 154 - GCOMObservation::drm(GModelSky&)
152 155 - GCOMResponse::irf_spatial()
153 156 - GCOMSelection::zeta_min()
... ... @@ -160,6 +163,10 @@ The following methods have been renamed:
160 163 The arguments for the following methods have been changed:
161 164 - GModelSpatial::flux(GSkyRegion&, GEnergy&, GTime&) was
162 165 GModelSpatial::flux(GSkyRegion*, GEnergy&, GTime&)
  166 +- GModelSpatialPointSource::GModelSpatialPointSource(GSkyDir&, std::string&) was
  167 + GModelSpatialPointSource::GModelSpatialPointSource(GSkyDir&)
  168 +- GModelSpatialPointSource::GModelSpatialPointSource(double&, double&, std::string&) was
  169 + GModelSpatialPointSource::GModelSpatialPointSource(double&, double&)
163 170 - GCOMObservation::drm(GSkyModel&) was GCOMObservation::drm(GSource&)
164 171  
165 172 The return value of the following methods has been changed:
... ... @@ -228,6 +235,12 @@ Implemented GMatrixSparse::row() and GMatrixSparse::add_to_row() methods.
228 235  
229 236 9. Model module
230 237 ---------------
  238 +Added support for Galactic coordinates in spatial models. The ra() and dec()
  239 +methods have been removed for spatial models, access to the coordinates is now
  240 +only supported via the dir() method. Models with celestial coordinates have
  241 +the parameters "RA" and "DEC", models with Galactic coordinates have the
  242 +parameters "GLON" and "GLAT" (#3988).
  243 +
231 244 Added GModelSpatialRadialGeneralGauss and GModelSpatialEllipticalGeneralGauss
232 245 classes to model general Gaussian functions as spatial components (#3924).
233 246  
... ...
README.md
  View file @c16adc2
1 1 GammaLib information
2 2 ====================
3   -* Version: 2.0.0.dev (19 January 2022)
  3 +* Version: 2.0.0.dev (11 February 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  
... ...
doc/source/admin/release_history/2.0.rst
  View file @c16adc2
1 1 .. _2.0:
2 2  
3   -GammaLib 2.0 release (19 January 2022)
4   -======================================
  3 +GammaLib 2.0 release (11 February 2022)
  4 +=======================================
5 5  
6 6 Overview
7 7 --------
... ... @@ -14,6 +14,7 @@ In particular, this release provides:
14 14 * full support of all spatial models for COMPTEL
15 15 * support for pulsar analysis
16 16 * vectorised response computation
  17 +* support Galactic coordinates in spatial models
17 18  
18 19  
19 20 Bug fixes
... ... @@ -40,6 +41,8 @@ Bug fixes
40 41 Improvements
41 42 ------------
42 43  
  44 +* [`3988 <https://cta-redmine.irap.omp.eu/issues/3988>`_] -
  45 + Support Galactic coordinates in spatial models
43 46 * Add ``GCOMObservation::tim()`` and ``GCOMObservation::oads()`` setter methods
44 47 * [`3917 <https://cta-redmine.irap.omp.eu/issues/3913>`_] -
45 48 Add ``GCOMSelection::fpmtflag()``, ``GCOMSelection::use_d1()`` and ``GCOMSelection::use_d2()`` methods
... ...
include/GModelSpatialPointSource.hpp
  View file @c16adc2
1 1 /***************************************************************************
2 2 * GModelSpatialPointSource.hpp - Spatial point source model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2021 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2022 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -32,8 +32,14 @@
32 32 #include "GModelSpatial.hpp"
33 33 #include "GModelPar.hpp"
34 34 #include "GSkyDir.hpp"
35   -#include "GSkyRegionCircle.hpp"
36   -#include "GXmlElement.hpp"
  35 +
  36 +/* __ Forward declarations _______________________________________________ */
  37 +class GEnergy;
  38 +class GTime;
  39 +class GPhoton;
  40 +class GRan;
  41 +class GSkyRegion;
  42 +class GXmlElement;
37 43  
38 44  
39 45 /***********************************************************************//**
... ... @@ -42,8 +48,11 @@
42 48 * @brief Point source spatial model
43 49 *
44 50 * This class implements a point source as the spatial component of the
45   - * factorised source model. The point source has two parameters: the Right
46   - * Ascension and Declination of the point source location.
  51 + * factorised source model. The point source has two parameters: the
  52 + * longitude and latitude of the point source location. The parameters may
  53 + * either be specified in celestial or Galactic coordinates. For celestial
  54 + * coordinates the parameter names are "RA" and "DEC", for Galactic
  55 + * coordinates the parameter names are "GLON" and "GLAT".
47 56 ***************************************************************************/
48 57 class GModelSpatialPointSource : public GModelSpatial {
49 58  
... ... @@ -51,8 +60,11 @@ public:
51 60 // Constructors and destructors
52 61 GModelSpatialPointSource(void);
53 62 GModelSpatialPointSource(const bool& dummy, const std::string& type);
54   - explicit GModelSpatialPointSource(const GSkyDir& dir);
55   - GModelSpatialPointSource(const double& ra, const double& dec);
  63 + GModelSpatialPointSource(const GSkyDir& dir,
  64 + const std::string& coordsys = "CEL");
  65 + GModelSpatialPointSource(const double& lon,
  66 + const double& lat,
  67 + const std::string& coordsys = "CEL");
56 68 explicit GModelSpatialPointSource(const GXmlElement& xml);
57 69 GModelSpatialPointSource(const GModelSpatialPointSource& model);
58 70 virtual ~GModelSpatialPointSource(void);
... ... @@ -84,23 +96,26 @@ public:
84 96 const GTime& srcTime = GTime()) const;
85 97  
86 98 // Other methods
87   - double ra(void) const;
88   - double dec(void) const;
89   - void ra(const double& ra);
90   - void dec(const double& dec);
91   - GSkyDir dir(void) const;
92   - void dir(const GSkyDir& dir);
  99 + std::string coordsys(void) const;
  100 + const GSkyDir& dir(void) const;
  101 + void dir(const GSkyDir& dir);
93 102  
94 103 protected:
95 104 // Protected methods
96 105 void init_members(void);
97 106 void copy_members(const GModelSpatialPointSource& model);
98 107 void free_members(void);
  108 + bool is_celestial(void) const;
99 109 virtual void set_region(void) const;
100 110  
101 111 // Protected members
102   - GModelPar m_ra; //!< Right Ascension (deg)
103   - GModelPar m_dec; //!< Declination (deg)
  112 + GModelPar m_lon; //!< Right Ascension or Galactic longitude (deg)
  113 + GModelPar m_lat; //!< Declination or Galactic latitude (deg)
  114 +
  115 + // Cached members for sky direction handling
  116 + mutable GSkyDir m_dir; //!< Sky direction representing parameters
  117 + mutable double m_last_lon; //!< Last longitude
  118 + mutable double m_last_lat; //!< Last latitude
104 119 };
105 120  
106 121  
... ... @@ -131,80 +146,49 @@ GClassCode GModelSpatialPointSource::code(void) const
131 146  
132 147  
133 148 /***********************************************************************//**
134   - * @brief Return Right Ascencion of model centre
135   - *
136   - * @return Right Ascencion of model centre (degrees).
137   - *
138   - * Returns the Right Ascension of the model centre in degrees.
139   - ***************************************************************************/
140   -inline
141   -double GModelSpatialPointSource::ra(void) const
142   -{
143   - return (m_ra.value());
144   -}
145   -
146   -
147   -/***********************************************************************//**
148   - * @brief Set Right Ascencion of model centre
149   - *
150   - * @param[in] ra Right Ascencion of model centre.
151   - *
152   - * Sets the Right Ascencion of model centre.
153   - ***************************************************************************/
154   -inline
155   -void GModelSpatialPointSource::ra(const double& ra)
156   -{
157   - m_ra.value(ra);
158   - return;
159   -}
160   -
161   -
162   -/***********************************************************************//**
163   - * @brief Return Declination of model centre
  149 + * @brief Return normalization of point source for Monte Carlo simulations
164 150 *
165   - * @return Declination of model centre (degrees).
  151 + * @param[in] dir Centre of simulation cone.
  152 + * @param[in] radius Radius of simulation cone (degrees).
  153 + * @return Normalization.
166 154 *
167   - * Returns the Declination of the model centre in degrees.
  155 + * Returns the normalization for a point source within a circular region.
  156 + * The normalization is 1 if the point source falls within the circle
  157 + * defined by @p dir and @p radius, 0 otherwise.
168 158 ***************************************************************************/
169 159 inline
170   -double GModelSpatialPointSource::dec(void) const
  160 +double GModelSpatialPointSource::mc_norm(const GSkyDir& dir,
  161 + const double& radius) const
171 162 {
172   - return (m_dec.value());
  163 + double norm = (dir.dist_deg(this->dir()) <= radius) ? 1.0 : 0.0;
  164 + return (norm);
173 165 }
174 166  
175 167  
176 168 /***********************************************************************//**
177   - * @brief Set Declination of model centre
  169 + * @brief Return coordinate system
178 170 *
179   - * @param[in] dec Declination of model centre.
  171 + * @return Coordinate system of point source model.
180 172 *
181   - * Sets the Declination of model centre.
  173 + * Returns "CEL" for a celestial coordinate system and "GAL" for a Galactic
  174 + * coordinate system.
182 175 ***************************************************************************/
183 176 inline
184   -void GModelSpatialPointSource::dec(const double& dec)
  177 +std::string GModelSpatialPointSource::coordsys(void) const
185 178 {
186   - m_dec.value(dec);
187   - return;
  179 + return (is_celestial() ? "CEL" : "GAL");
188 180 }
189 181  
190 182  
191 183 /***********************************************************************//**
192   - * @brief Return normalization of point source for Monte Carlo simulations
  184 + * @brief Check if model holds celestial coordinates
193 185 *
194   - * @param[in] dir Centre of simulation cone.
195   - * @param[in] radius Radius of simulation cone (degrees).
196   - * @return Normalization.
197   - *
198   - * Returns the normalization for a point source within a circular region.
199   - * The normalization is 1 if the point source falls within the circle
200   - * defined by @p dir and @p radius, 0 otherwise.
  186 + * @return True if model holds celestial coordinates, false otherwise.
201 187 ***************************************************************************/
202 188 inline
203   -double GModelSpatialPointSource::mc_norm(const GSkyDir& dir,
204   - const double& radius) const
  189 +bool GModelSpatialPointSource::is_celestial(void) const
205 190 {
206   - double norm = (dir.dist_deg(this->dir()) <= radius) ? 1.0 : 0.0;
207   - return (norm);
  191 + return (m_lon.name() == "RA");
208 192 }
209 193  
210 194 #endif /* GMODELSPATIALPOINTSOURCE_HPP */
... ...
pyext/GModelSpatialPointSource.i
  View file @c16adc2
1 1 /***************************************************************************
2 2 * GModelSpatialPointSource.i - Spatial point source model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2021 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2022 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -39,8 +39,11 @@ class GModelSpatialPointSource : public GModelSpatial {
39 39 public:
40 40 // Constructors and destructors
41 41 GModelSpatialPointSource(void);
42   - explicit GModelSpatialPointSource(const GSkyDir& dir);
43   - GModelSpatialPointSource(const double& ra, const double& dec);
  42 + GModelSpatialPointSource(const GSkyDir& dir,
  43 + const std::string& coordsys = "CEL");
  44 + GModelSpatialPointSource(const double& lon,
  45 + const double& lat,
  46 + const std::string& coordsys = "CEL");
44 47 explicit GModelSpatialPointSource(const GXmlElement& xml);
45 48 GModelSpatialPointSource(const GModelSpatialPointSource& model);
46 49 virtual ~GModelSpatialPointSource(void);
... ... @@ -67,12 +70,9 @@ public:
67 70 const GTime& srcTime = GTime()) const;
68 71  
69 72 // Other methods
70   - double ra(void) const;
71   - double dec(void) const;
72   - void ra(const double& ra);
73   - void dec(const double& dec);
74   - GSkyDir dir(void) const;
75   - void dir(const GSkyDir& dir);
  73 + std::string coordsys(void) const;
  74 + const GSkyDir& dir(void) const;
  75 + void dir(const GSkyDir& dir);
76 76 };
77 77  
78 78  
... ... @@ -85,12 +85,10 @@ public:
85 85 }
86 86 %pythoncode {
87 87 def __getstate__(self):
88   - state = self.type(), self[0], self[1]
  88 + state = {'type': self.type(), 'dir': self.dir(), 'coordsys': self.coordsys()}
89 89 return state
90 90 def __setstate__(self, state):
91   - self.__init__()
92   - self.type(state[0])
93   - self[0] = state[1]
94   - self[1] = state[2]
  91 + self.__init__(state['dir'],state['coordsys'])
  92 + self.type(state['type'])
95 93 }
96 94 };
... ...
src/model/GModelSpatialPointSource.cpp
  View file @c16adc2
1 1 /***************************************************************************
2 2 * GModelSpatialPointSource.cpp - Spatial point source model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2021 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2022 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -47,6 +47,10 @@ const GModelSpatialRegistry g_spatial_ptsrc_legacy_registry(&amp;g_spatial_ptsrc_
47 47 #endif
48 48  
49 49 /* __ Method name definitions ____________________________________________ */
  50 +#define G_CONSTRUCTOR1 "GModelSpatialPointSource::GModelSpatialPointSource("\
  51 + "GSkyDir&, std::string&)"
  52 +#define G_CONSTRUCTOR2 "GModelSpatialPointSource::GModelSpatialPointSource("\
  53 + "double&, double&, std::string&)"
50 54 #define G_READ "GModelSpatialPointSource::read(GXmlElement&)"
51 55 #define G_WRITE "GModelSpatialPointSource::write(GXmlElement&)"
52 56  
... ... @@ -105,15 +109,40 @@ GModelSpatialPointSource::GModelSpatialPointSource(const bool&amp; dummy,
105 109 * @brief Sky direction constructor
106 110 *
107 111 * @param[in] dir Sky direction.
  112 + * @param[in] coordsys Coordinate system (either "CEL" or "GAL")
108 113 *
109   - * Construct a point source spatial model from a sky direction.
  114 + * @exception GException::invalid_argument
  115 + * Invalid @p coordsys argument specified.
  116 + *
  117 + * Construct a point source spatial model from a sky direction. The
  118 + * @p coordsys parameter specified whether the sky direction should be
  119 + * interpreted in the celestial or Galactic coordinata system.
110 120 ***************************************************************************/
111   -GModelSpatialPointSource::GModelSpatialPointSource(const GSkyDir& dir) :
  121 +GModelSpatialPointSource::GModelSpatialPointSource(const GSkyDir& dir,
  122 + const std::string& coordsys) :
112 123 GModelSpatial()
113 124 {
  125 + // Throw an exception if the coordinate system is invalid
  126 + if ((coordsys != "CEL") && (coordsys != "GAL")) {
  127 + std::string msg = "Invalid coordinate system \""+coordsys+"\" "
  128 + "specified. Please specify either \"CEL\" or "
  129 + "\"GAL\".";
  130 + throw GException::invalid_argument(G_CONSTRUCTOR1, msg);
  131 + }
  132 +
114 133 // Initialise members
115 134 init_members();
116 135  
  136 + // Set parameter names
  137 + if (coordsys == "CEL") {
  138 + m_lon.name("RA");
  139 + m_lat.name("DEC");
  140 + }
  141 + else {
  142 + m_lon.name("GLON");
  143 + m_lat.name("GLAT");
  144 + }
  145 +
117 146 // Assign direction
118 147 this->dir(dir);
119 148  
... ... @@ -125,22 +154,49 @@ GModelSpatialPointSource::GModelSpatialPointSource(const GSkyDir&amp; dir) :
125 154 /***********************************************************************//**
126 155 * @brief Value constructor
127 156 *
128   - * @param[in] ra Right Ascencion of model centre.
129   - * @param[in] dec Declination of model centre.
  157 + * @param[in] lon Longitude of point source (deg).
  158 + * @param[in] lat Latitude of point source (deg).
  159 + * @param[in] coordsys Coordinate system (either "CEL" or "GAL")
130 160 *
131   - * Construct a point source spatial model from the Right Ascension and
132   - * Declination of the model centre.
  161 + * @exception GException::invalid_argument
  162 + * Invalid @p coordsys argument specified.
  163 + *
  164 + * Construct a point source spatial model from the longitude and latitude
  165 + * of a point source. The @p coordsys parameter specifies the coordinate
  166 + * system in which @p lon and @p lat are provided. By default a celestial
  167 + * coordinate system is assumed, which means that @p lon and @p lat are
  168 + * taken as Right Ascension and Declination. If "GAL" is specified, @p lon
  169 + * and @p lat are taken as Galactic longitude and latitude.
133 170 ***************************************************************************/
134   -GModelSpatialPointSource::GModelSpatialPointSource(const double& ra,
135   - const double& dec) :
  171 +GModelSpatialPointSource::GModelSpatialPointSource(const double& lon,
  172 + const double& lat,
  173 + const std::string& coordsys) :
136 174 GModelSpatial()
137 175 {
  176 + // Throw an exception if the coordinate system is invalid
  177 + if ((coordsys != "CEL") && (coordsys != "GAL")) {
  178 + std::string msg = "Invalid coordinate system \""+coordsys+"\" "
  179 + "specified. Please specify either \"CEL\" or "
  180 + "\"GAL\".";
  181 + throw GException::invalid_argument(G_CONSTRUCTOR2, msg);
  182 + }
  183 +
138 184 // Initialise members
139 185 init_members();
140 186  
  187 + // Set parameter names
  188 + if (coordsys == "CEL") {
  189 + m_lon.name("RA");
  190 + m_lat.name("DEC");
  191 + }
  192 + else {
  193 + m_lon.name("GLON");
  194 + m_lat.name("GLAT");
  195 + }
  196 +
141 197 // Set values
142   - m_ra.value(ra);
143   - m_dec.value(dec);
  198 + m_lon.value(lon);
  199 + m_lat.value(lat);
144 200  
145 201 // Return
146 202 return;
... ... @@ -352,15 +408,15 @@ bool GModelSpatialPointSource::contains(const GSkyDir&amp; dir,
352 408 * format
353 409 *
354 410 * <spatialModel type="PointSource">
355   - * <parameter free="0" max="360" min="-360" name="RA" scale="1" value="83.6331" />
356   - * <parameter free="0" max="90" min="-90" name="DEC" scale="1" value="22.0145" />
  411 + * <parameter name="RA" scale="1" value="83.6331" min="-360" max="360" free="0" />
  412 + * <parameter name="DEC" scale="1" value="22.0145" min="-90" max="90" free="0" />
357 413 * </spatialModel>
358 414 *
359 415 * or
360 416 *
361 417 * <spatialModel type="PointSource">
362   - * <parameter free="0" max="360" min="-360" name="GLON" scale="1" value="83.6331" />
363   - * <parameter free="0" max="90" min="-90" name="GLAT" scale="1" value="22.0145" />
  418 + * <parameter name="GLON" scale="1" value="184.5575" min="-360" max="360" free="0" />
  419 + * <parameter name="GLAT" scale="1" value="-5.7843" min="-90" max="90" free="0" />
364 420 * </spatialModel>
365 421 *
366 422 ***************************************************************************/
... ... @@ -377,8 +433,8 @@ void GModelSpatialPointSource::read(const GXmlElement&amp; xml)
377 433 const GXmlElement* dec = gammalib::xml_get_par(G_READ, xml, "DEC");
378 434  
379 435 // Read parameters
380   - m_ra.read(*ra);
381   - m_dec.read(*dec);
  436 + m_lon.read(*ra);
  437 + m_lat.read(*dec);
382 438  
383 439 }
384 440  
... ... @@ -390,18 +446,8 @@ void GModelSpatialPointSource::read(const GXmlElement&amp; xml)
390 446 const GXmlElement* glat = gammalib::xml_get_par(G_READ, xml, "GLAT");
391 447  
392 448 // Read parameters
393   - m_ra.read(*glon);
394   - m_dec.read(*glat);
395   -
396   - // Convert into RA/DEC
397   - GSkyDir dir;
398   - dir.lb_deg(ra(), dec()),
399   - m_ra.value(dir.ra_deg());
400   - m_dec.value(dir.dec_deg());
401   -
402   - // Set names to RA/DEC
403   - m_ra.name("RA");
404   - m_dec.name("DEC");
  449 + m_lon.read(*glon);
  450 + m_lat.read(*glat);
405 451  
406 452 }
407 453  
... ... @@ -415,16 +461,21 @@ void GModelSpatialPointSource::read(const GXmlElement&amp; xml)
415 461 *
416 462 * @param[in] xml XML element into which model information is written.
417 463 *
418   - * Write the point source information into an XML element with the following
419   - * format
  464 + * Depending on the coordinate system, write the point source information
  465 + * into an XML element with the following format
420 466 *
421 467 * <spatialModel type="PointSource">
422   - * <parameter free="0" max="360" min="-360" name="RA" scale="1" value="83.6331" />
423   - * <parameter free="0" max="90" min="-90" name="DEC" scale="1" value="22.0145" />
  468 + * <parameter name="RA" scale="1" value="83.6331" min="-360" max="360" free="0" />
  469 + * <parameter name="DEC" scale="1" value="22.0145" min="-90" max="90" free="0" />
  470 + * </spatialModel>
  471 + *
  472 + * or
  473 + *
  474 + * <spatialModel type="PointSource">
  475 + * <parameter name="GLON" scale="1" value="184.5575" min="-360" max="360" free="0" />
  476 + * <parameter name="GLAT" scale="1" value="-5.7843" min="-90" max="90" free="0" />
424 477 * </spatialModel>
425 478 *
426   - * @todo The case that an existing spatial XML element with "GLON" and "GLAT"
427   - * as coordinates is not supported.
428 479 ***************************************************************************/
429 480 void GModelSpatialPointSource::write(GXmlElement& xml) const
430 481 {
... ... @@ -432,12 +483,12 @@ void GModelSpatialPointSource::write(GXmlElement&amp; xml) const
432 483 gammalib::xml_check_type(G_WRITE, xml, type());
433 484  
434 485 // Get or create parameters
435   - GXmlElement* ra = gammalib::xml_need_par(G_WRITE, xml, m_ra.name());
436   - GXmlElement* dec = gammalib::xml_need_par(G_WRITE, xml, m_dec.name());
  486 + GXmlElement* lon = gammalib::xml_need_par(G_WRITE, xml, m_lon.name());
  487 + GXmlElement* lat = gammalib::xml_need_par(G_WRITE, xml, m_lat.name());
437 488  
438 489 // Write parameters
439   - m_ra.write(*ra);
440   - m_dec.write(*dec);
  490 + m_lon.write(*lon);
  491 + m_lat.write(*lat);
441 492  
442 493 // Return
443 494 return;
... ... @@ -511,29 +562,53 @@ std::string GModelSpatialPointSource::print(const GChatter&amp; chatter) const
511 562 *
512 563 * Returns the sky direction of the point source.
513 564 ***************************************************************************/
514   -GSkyDir GModelSpatialPointSource::dir(void) const
  565 +const GSkyDir& GModelSpatialPointSource::dir(void) const
515 566 {
516   - // Allocate sky direction
517   - GSkyDir srcDir;
  567 + // Get longitude and latitude values
  568 + double lon = m_lon.value();
  569 + double lat = m_lat.value();
518 570  
519   - // Set sky direction
520   - srcDir.radec_deg(ra(), dec());
  571 + // If longitude or latitude values have changed then update sky
  572 + // direction cache
  573 + if ((lon != m_last_lon) || (lat != m_last_lat)) {
  574 +
  575 + // Update last values
  576 + m_last_lon = lon;
  577 + m_last_lat = lat;
  578 +
  579 + // Update sky direction dependent on model coordinate system
  580 + if (is_celestial()) {
  581 + m_dir.radec_deg(m_last_lon, m_last_lat);
  582 + }
  583 + else {
  584 + m_dir.lb_deg(m_last_lon, m_last_lat);
  585 + }
521 586  
522   - // Return direction
523   - return srcDir;
  587 + } // endif: update of sky direction cache required
  588 +
  589 + // Return sky direction
  590 + return (m_dir);
524 591 }
525 592  
526 593  
527 594 /***********************************************************************//**
528 595 * @brief Set position of point source
529 596 *
  597 + * @param[in] dir Sky direction of point source.
  598 + *
530 599 * Sets the sky direction of the point source.
531 600 ***************************************************************************/
532 601 void GModelSpatialPointSource::dir(const GSkyDir& dir)
533 602 {
534   - // Assign Right Ascension and Declination
535   - m_ra.value(dir.ra_deg());
536   - m_dec.value(dir.dec_deg());
  603 + // Assign sky direction depending on the model coordinate system
  604 + if (is_celestial()) {
  605 + m_lon.value(dir.ra_deg());
  606 + m_lat.value(dir.dec_deg());
  607 + }
  608 + else {
  609 + m_lon.value(dir.l_deg());
  610 + m_lat.value(dir.b_deg());
  611 + }
537 612  
538 613 // Return
539 614 return;
... ... @@ -555,27 +630,32 @@ void GModelSpatialPointSource::init_members(void)
555 630 m_type = "PointSource";
556 631  
557 632 // Initialise Right Ascension
558   - m_ra.clear();
559   - m_ra.name("RA");
560   - m_ra.unit("deg");
561   - m_ra.fix();
562   - m_ra.scale(1.0);
563   - m_ra.gradient(0.0);
564   - m_ra.has_grad(false);
  633 + m_lon.clear();
  634 + m_lon.name("RA");
  635 + m_lon.unit("deg");
  636 + m_lon.fix();
  637 + m_lon.scale(1.0);
  638 + m_lon.gradient(0.0);
  639 + m_lon.has_grad(false);
565 640  
566 641 // Initialise Declination
567   - m_dec.clear();
568   - m_dec.name("DEC");
569   - m_dec.unit("deg");
570   - m_dec.fix();
571   - m_dec.scale(1.0);
572   - m_dec.gradient(0.0);
573   - m_dec.has_grad(false);
  642 + m_lat.clear();
  643 + m_lat.name("DEC");
  644 + m_lat.unit("deg");
  645 + m_lat.fix();
  646 + m_lat.scale(1.0);
  647 + m_lat.gradient(0.0);
  648 + m_lat.has_grad(false);
574 649  
575 650 // Set parameter pointer(s)
576 651 m_pars.clear();
577   - m_pars.push_back(&m_ra);
578   - m_pars.push_back(&m_dec);
  652 + m_pars.push_back(&m_lon);
  653 + m_pars.push_back(&m_lat);
  654 +
  655 + // Initialise cache
  656 + m_dir.clear();
  657 + m_last_lon = -9999.0;
  658 + m_last_lat = -9999.0;
579 659  
580 660 // Return
581 661 return;
... ... @@ -591,13 +671,18 @@ void GModelSpatialPointSource::copy_members(const GModelSpatialPointSource&amp; mode
591 671 {
592 672 // Copy members
593 673 m_type = model.m_type; // Needed to conserve model type
594   - m_ra = model.m_ra;
595   - m_dec = model.m_dec;
  674 + m_lon = model.m_lon;
  675 + m_lat = model.m_lat;
596 676  
597 677 // Set parameter pointer(s)
598 678 m_pars.clear();
599   - m_pars.push_back(&m_ra);
600   - m_pars.push_back(&m_dec);
  679 + m_pars.push_back(&m_lon);
  680 + m_pars.push_back(&m_lat);
  681 +
  682 + // Copy cache
  683 + m_dir = model.m_dir;
  684 + m_last_lon = model.m_last_lon;
  685 + m_last_lat = model.m_last_lat;
601 686  
602 687 // Return
603 688 return;
... ... @@ -620,7 +705,7 @@ void GModelSpatialPointSource::free_members(void)
620 705 void GModelSpatialPointSource::set_region(void) const
621 706 {
622 707 // Set sky region circle
623   - GSkyRegionCircle region(m_ra.value(), m_dec.value(), 0.0);
  708 + GSkyRegionCircle region(dir(), 0.0);
624 709  
625 710 // Set region (circumvent const correctness)
626 711 const_cast<GModelSpatialPointSource*>(this)->m_region = region;
... ...
test/test_GModel.cpp
  View file @c16adc2
... ... @@ -882,78 +882,90 @@ void TestGModel::test_point_source(void)
882 882 GModelSpatialPointSource model1;
883 883 test_value(model1.type(), "PointSource");
884 884  
885   - // Test sky direction constructor
  885 + // Test celestial sky direction constructor
886 886 GSkyDir dir1;
887 887 dir1.radec_deg(83.6331, +22.0145);
888 888 GModelSpatialPointSource model2(dir1);
889   - test_value(model2.ra(), 83.6331);
890   - test_value(model2.dec(), +22.0145);
891   -
892   - // Test value constructor
893   - GModelSpatialPointSource model3(83.6331, +22.0145);
894   - test_value(model3.ra(), 83.6331);
895   - test_value(model3.dec(), +22.0145);
  889 + test_value(model2.dir().ra_deg(), 83.6331, "Test Right Ascension for celestial sky direction constructor");
  890 + test_value(model2.dir().dec_deg(), +22.0145, "Test Declination for celestial sky direction constructor");
  891 + test_value(model2.dir().l_deg(), 184.557460, "Test Galactic longitude for celestial sky direction constructor");
  892 + test_value(model2.dir().b_deg(), -5.784346, "Test Galactic latitude for celestial sky direction constructor");
  893 +
  894 + // Test Galactic sky direction constructor
  895 + GModelSpatialPointSource model3(dir1, "GAL");
  896 + test_value(model3.dir().ra_deg(), 83.6331, "Test Right Ascension for Galactic sky direction constructor");
  897 + test_value(model3.dir().dec_deg(), +22.0145, "Test Declination for Galactic sky direction constructor");
  898 + test_value(model3.dir().l_deg(), 184.557460, "Test Galactic longitude for Galactic sky direction constructor");
  899 + test_value(model3.dir().b_deg(), -5.784346, "Test Galactic latitude for Galactic sky direction constructor");
  900 +
  901 + // Test celestial value constructor
  902 + GModelSpatialPointSource model4(83.6331, +22.0145);
  903 + test_value(model4.dir().ra_deg(), 83.6331, "Test Right Ascension for celestial value constructor");
  904 + test_value(model4.dir().dec_deg(), +22.0145, "Test Declination for celestial value constructor");
  905 + test_value(model4.dir().l_deg(), 184.557460, "Test Galactic longitude for celestial value constructor");
  906 + test_value(model4.dir().b_deg(), -5.784346, "Test Galactic latitude for celestial value constructor");
  907 +
  908 + // Test Galactic value constructor
  909 + GModelSpatialPointSource model5(184.5575, -5.7843, "GAL");
  910 + test_value(model5.dir().ra_deg(), 83.633165, "Test Right Ascension for Galactic value constructor");
  911 + test_value(model5.dir().dec_deg(), +22.014491, "Test Declination for Galactic value constructor");
  912 + test_value(model5.dir().l_deg(), 184.5575, "Test Galactic longitude for Galactic value constructor");
  913 + test_value(model5.dir().b_deg(), -5.7843, "Test Galactic latitude for Galactic value constructor");
896 914  
897 915 // Test XML constructor
898 916 GXml xml(m_xml_model_point_plaw);
899 917 GXmlElement* element = xml.element(0)->element(0)->element("spatialModel", 0);
900   - GModelSpatialPointSource model4(*element);
901   - test_value(model4.size(), 2);
902   - test_value(model4.type(), "PointSource");
903   - test_value(model4.ra(), 83.6331);
904   - test_value(model4.dec(), +22.0145);
905   -
906   - // Test ra method
907   - model4.ra(3.9);
908   - test_value(model4.ra(), 3.9);
909   -
910   - // Test dec method
911   - model4.dec(3.9);
912   - test_value(model4.dec(), 3.9);
  918 + GModelSpatialPointSource model6(*element);
  919 + test_value(model6.size(), 2);
  920 + test_value(model6.type(), "PointSource");
  921 + test_value(model6.dir().ra_deg(), 83.6331, "Test Right Ascension for celestial XML constructor");
  922 + test_value(model6.dir().dec_deg(), +22.0145, "Test Declination for celestial XML constructor");
  923 + test_value(model6.dir().l_deg(), 184.557460, "Test Galactic longitude for celestial XML constructor");
  924 + test_value(model6.dir().b_deg(), -5.784346, "Test Galactic latitude for celestial XML constructor");
913 925  
914 926 // Test region method
915 927 GSkyDir dir2;
916   - dir2.radec_deg(3.9, 3.9);
917   - test_assert(model4.region()->contains(dir2), "Test region() method (inside)");
918   - dir2.radec_deg(4.9, 3.9);
919   - test_assert(!model4.region()->contains(dir2), "Test region() method (outside)");
  928 + dir2.radec_deg(83.6331, +22.0145);
  929 + test_assert(model6.region()->contains(dir2), "Test region() method (inside)");
  930 + dir2.radec_deg(84.6331, +22.0145);
  931 + test_assert(!model6.region()->contains(dir2), "Test region() method (outside)");
920 932  
921 933 // Test operator access
922 934 const char* strarray[] = {"RA", "DEC"};
923 935 for (int i = 0; i < 2; ++i) {
924 936 std::string keyname(strarray[i]);
925   - model4[keyname].value(2.1);
926   - model4[keyname].error(1.9);
927   - model4[keyname].gradient(0.8);
928   - test_value(model4[keyname].value(), 2.1);
929   - test_value(model4[keyname].error(), 1.9);
930   - test_value(model4[keyname].gradient(), 0.8);
  937 + model6[keyname].value(2.1);
  938 + model6[keyname].error(1.9);
  939 + model6[keyname].gradient(0.8);
  940 + test_value(model6[keyname].value(), 2.1);
  941 + test_value(model6[keyname].error(), 1.9);
  942 + test_value(model6[keyname].gradient(), 0.8);
931 943 }
932 944  
933 945 // Test size() method
934   - GModelSpatialPointSource model5;
935   - test_value(model5.size(), 2, "Test size() method.");
  946 + GModelSpatialPointSource model7;
  947 + test_value(model7.size(), 2, "Test size() method.");
936 948  
937 949 // Test has_par() method
938   - test_assert(model5.has_par("RA"), "Test has_par() method for \"RA\"");
939   - test_assert(model5.has_par("DEC"), "Test has_par() method for \"DEC\"");
940   - test_assert(!model5.has_par("XXX"), "Test has_par() method for \"XXX\"");
  950 + test_assert(model7.has_par("RA"), "Test has_par() method for \"RA\"");
  951 + test_assert(model7.has_par("DEC"), "Test has_par() method for \"DEC\"");
  952 + test_assert(!model7.has_par("XXX"), "Test has_par() method for \"XXX\"");
941 953  
942 954 // Test has_free_pars() method
943   - model5["RA"].fix();
944   - model5["DEC"].fix();
945   - test_assert(!model5.has_free_pars(), "Test has_free_pars() method for fixed parameters");
946   - model5["DEC"].free();
947   - test_assert(model5.has_free_pars(), "Test has_free_pars() method for free parameters");
  955 + model7["RA"].fix();
  956 + model7["DEC"].fix();
  957 + test_assert(!model7.has_free_pars(), "Test has_free_pars() method for fixed parameters");
  958 + model7["DEC"].free();
  959 + test_assert(model7.has_free_pars(), "Test has_free_pars() method for free parameters");
948 960  
949 961 // Test autoscale() method
950   - model5["RA"].value(180.0);
951   - model5["DEC"].value(-30.0);
952   - model5.autoscale();
953   - test_value(model5["RA"].scale(), 180.0, "Test autoscale() method (\"RA\" scale).");
954   - test_value(model5["DEC"].scale(), -30.0, "Test autoscale() method (\"DEC\" scale).");
955   - test_value(model5["RA"].factor_value(), 1.0, "Test autoscale() method (\"RA\" factor).");
956   - test_value(model5["DEC"].factor_value(), 1.0, "Test autoscale() method (\"DEC\" factor).");
  962 + model7["RA"].value(180.0);
  963 + model7["DEC"].value(-30.0);
  964 + model7.autoscale();
  965 + test_value(model7["RA"].scale(), 180.0, "Test autoscale() method (\"RA\" scale).");
  966 + test_value(model7["DEC"].scale(), -30.0, "Test autoscale() method (\"DEC\" scale).");
  967 + test_value(model7["RA"].factor_value(), 1.0, "Test autoscale() method (\"RA\" factor).");
  968 + test_value(model7["DEC"].factor_value(), 1.0, "Test autoscale() method (\"DEC\" factor).");
957 969  
958 970 // Exit test
959 971 return;
... ... @@ -4173,8 +4185,8 @@ void TestGModel::test_legacy_model_point_const(void)
4173 4185  
4174 4186 // Test model values
4175 4187 test_value(spatial->type(), "SkyDirFunction", "Check spatial model type");
4176   - test_value(spatial->ra(), 83.6331, 1.0e-7, "Check Right Ascension");
4177   - test_value(spatial->dec(), 22.0145, 1.0e-7, "Check Declination");
  4188 + test_value(spatial->dir().ra_deg(), 83.6331, 1.0e-7, "Check Right Ascension");
  4189 + test_value(spatial->dir().dec_deg(), 22.0145, 1.0e-7, "Check Declination");
4178 4190 test_value(spectral->type(), "ConstantValue", "Check spectral model type");
4179 4191 test_value(spectral->value(), 5.7e-16, 1.0e-7, "Check constant value");
4180 4192  
... ... @@ -4193,8 +4205,8 @@ void TestGModel::test_legacy_model_point_const(void)
4193 4205  
4194 4206 // Re-test model values
4195 4207 test_value(spatial->type(), "SkyDirFunction", "Check spatial model type");
4196   - test_value(spatial->ra(), 83.6331, 1.0e-7, "Check Right Ascension");
4197   - test_value(spatial->dec(), 22.0145, 1.0e-7, "Check Declination");
  4208 + test_value(spatial->dir().ra_deg(), 83.6331, 1.0e-7, "Check Right Ascension");
  4209 + test_value(spatial->dir().dec_deg(), 22.0145, 1.0e-7, "Check Declination");
4198 4210 test_value(spectral->type(), "ConstantValue", "Check spectral model type");
4199 4211 test_value(spectral->value(), 5.7e-16, 1.0e-7, "Check constant value");
4200 4212  
... ...