Jürgen Knödlseder / gammalib

Commit df5043453de66d1fb424a2a2b0a62aff410aa24f

Authored by Jürgen Knödlseder
1 parent c16adc29

Support Galactic coordinates in radial spatial models

Inline Side-by-side
Showing 25 changed files with 571 additions and 372 deletions
... ... @@ -81,6 +81,7 @@ The following methods have been added:
81 81 - GModelSpatialPointSource::coordsys()
82 82 - GModelSpatialRadial::is_energy_dependent()
83 83 - GModelSpatialRadial::is_time_dependent()
  84 +- GModelSpatialRadial::coordsys()
84 85 - GModelSpectralFunc(GModelSpectral&, GEnergies&)
85 86 - GModelSpectralFunc::nodes()
86 87 - GModelSpectralFunc::is_empty()
... ... @@ -151,6 +152,8 @@ The following methods have been removed:
151 152 - GFitsTableCol::scale()
152 153 - GModelSpatialPointSource::ra()
153 154 - GModelSpatialPointSource::dec()
  155 +- GModelSpatialRadial::ra()
  156 +- GModelSpatialRadial::dec()
154 157 - GCOMObservation::drm(GModelSky&)
155 158 - GCOMResponse::irf_spatial()
156 159 - GCOMSelection::zeta_min()
... ...
include/GModelSpatialRadial.hpp
  View file @df50434
1 1 /***************************************************************************
2 2 * GModelSpatialRadial.hpp - Abstract radial spatial model base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2020 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2022 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -33,11 +33,14 @@
33 33 #include "GModelSpatial.hpp"
34 34 #include "GModelPar.hpp"
35 35 #include "GSkyDir.hpp"
36   -#include "GEnergy.hpp"
37   -#include "GTime.hpp"
38   -#include "GPhoton.hpp"
39   -#include "GXmlElement.hpp"
40   -#include "GRan.hpp"
  36 +
  37 +/* __ Forward declarations _______________________________________________ */
  38 +class GEnergy;
  39 +class GTime;
  40 +class GPhoton;
  41 +class GRan;
  42 +class GSkyRegion;
  43 +class GXmlElement;
41 44  
42 45  
43 46 /***********************************************************************//**
... ... @@ -88,23 +91,26 @@ public:
88 91 virtual void write(GXmlElement& xml) const;
89 92  
90 93 // Other methods
91   - double ra(void) const;
92   - double dec(void) const;
93   - void ra(const double& ra);
94   - void dec(const double& dec);
95   - GSkyDir dir(void) const;
96   - void dir(const GSkyDir& dir);
  94 + std::string coordsys(void) const;
  95 + const GSkyDir& dir(void) const;
  96 + void dir(const GSkyDir& dir);
97 97  
98 98 protected:
99 99 // Protected methods
100 100 void init_members(void);
101 101 void copy_members(const GModelSpatialRadial& model);
102 102 void free_members(void);
  103 + bool is_celestial(void) const;
103 104 virtual void set_region(void) const = 0;
104 105  
105   - // Proteced members
106   - GModelPar m_ra; //!< Right Ascension (deg)
107   - GModelPar m_dec; //!< Declination (deg)
  106 + // Protected members
  107 + GModelPar m_lon; //!< Right Ascension or Galactic longitude (deg)
  108 + GModelPar m_lat; //!< Declination or Galactic latitude (deg)
  109 +
  110 + // Cached members for sky direction handling
  111 + mutable GSkyDir m_dir; //!< Sky direction representing parameters
  112 + mutable double m_last_lon; //!< Last longitude
  113 + mutable double m_last_lat; //!< Last latitude
108 114 };
109 115  
110 116  
... ... @@ -153,80 +159,49 @@ bool GModelSpatialRadial::is_time_dependent(void) const
153 159  
154 160  
155 161 /***********************************************************************//**
156   - * @brief Return Right Ascencion of model centre
157   - *
158   - * @return Right Ascencion of model centre (degrees).
159   - *
160   - * Returns the Right Ascension of the model centre in degrees.
161   - ***************************************************************************/
162   -inline
163   -double GModelSpatialRadial::ra(void) const
164   -{
165   - return (m_ra.value());
166   -}
167   -
168   -
169   -/***********************************************************************//**
170   - * @brief Set Right Ascencion of model centre
171   - *
172   - * @param[in] ra Right Ascension (degrees).
173   - *
174   - * Sets the Right Ascension of the model centre in degrees.
175   - ***************************************************************************/
176   -inline
177   -void GModelSpatialRadial::ra(const double& ra)
178   -{
179   - m_ra.value(ra);
180   - return;
181   -}
182   -
183   -
184   -/***********************************************************************//**
185   - * @brief Return Declination of model centre
  162 + * @brief Return normalization of radial source for Monte Carlo simulations
186 163 *
187   - * @return Declination of model centre (degrees).
  164 + * @param[in] dir Centre of simulation cone.
  165 + * @param[in] radius Radius of simulation cone (degrees).
  166 + * @return Normalization.
188 167 *
189   - * Returns the Declination of the model centre in degrees.
  168 + * Returns the normalization for a radial source within a circular region.
  169 + * The normalization is 1 if the radial source falls within the circle
  170 + * defined by @p dir and @p radius, 0 otherwise.
190 171 ***************************************************************************/
191 172 inline
192   -double GModelSpatialRadial::dec(void) const
  173 +double GModelSpatialRadial::mc_norm(const GSkyDir& dir,
  174 + const double& radius) const
193 175 {
194   - return (m_dec.value());
  176 + double norm = (dir.dist_deg(this->dir()) <= radius+theta_max()) ? 1.0 : 0.0;
  177 + return (norm);
195 178 }
196 179  
197 180  
198 181 /***********************************************************************//**
199   - * @brief Set Declination of model centre
  182 + * @brief Return coordinate system
200 183 *
201   - * @param[in] dec Declination (degrees).
  184 + * @return Coordinate system of point source model.
202 185 *
203   - * Sets the Declination of the model centre in degrees.
  186 + * Returns "CEL" for a celestial coordinate system and "GAL" for a Galactic
  187 + * coordinate system.
204 188 ***************************************************************************/
205 189 inline
206   -void GModelSpatialRadial::dec(const double& dec)
  190 +std::string GModelSpatialRadial::coordsys(void) const
207 191 {
208   - m_dec.value(dec);
209   - return;
  192 + return (is_celestial() ? "CEL" : "GAL");
210 193 }
211 194  
212 195  
213 196 /***********************************************************************//**
214   - * @brief Return normalization of radial source for Monte Carlo simulations
  197 + * @brief Check if model holds celestial coordinates
215 198 *
216   - * @param[in] dir Centre of simulation cone.
217   - * @param[in] radius Radius of simulation cone (degrees).
218   - * @return Normalization.
219   - *
220   - * Returns the normalization for a radial source within a circular region.
221   - * The normalization is 1 if the radial source falls within the circle
222   - * defined by @p dir and @p radius, 0 otherwise.
  199 + * @return True if model holds celestial coordinates, false otherwise.
223 200 ***************************************************************************/
224 201 inline
225   -double GModelSpatialRadial::mc_norm(const GSkyDir& dir,
226   - const double& radius) const
  202 +bool GModelSpatialRadial::is_celestial(void) const
227 203 {
228   - double norm = (dir.dist_deg(this->dir()) <= radius+theta_max()) ? 1.0 : 0.0;
229   - return (norm);
  204 + return (m_lon.name() == "RA");
230 205 }
231 206  
232 207 #endif /* GMODELSPATIALRADIAL_HPP */
... ...
include/GModelSpatialRadialDisk.hpp
  View file @df50434
1 1 /***************************************************************************
2 2 * GModelSpatialRadialDisk.hpp - Radial disk source model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2020 by Christoph Deil *
  4 + * copyright (C) 2011-2022 by Christoph Deil *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -31,9 +31,15 @@
31 31 #include <string>
32 32 #include "GModelSpatialRadial.hpp"
33 33 #include "GModelPar.hpp"
34   -#include "GSkyDir.hpp"
35   -#include "GSkyRegionCircle.hpp"
36   -#include "GXmlElement.hpp"
  34 +
  35 +/* __ Forward declarations _______________________________________________ */
  36 +class GEnergy;
  37 +class GTime;
  38 +class GPhoton;
  39 +class GRan;
  40 +class GSkyDir;
  41 +class GSkyRegion;
  42 +class GXmlElement;
37 43  
38 44  
39 45 /**************************************************************************
... ... @@ -51,7 +57,9 @@ public:
51 57 // Constructors and destructors
52 58 GModelSpatialRadialDisk(void);
53 59 GModelSpatialRadialDisk(const bool& dummy, const std::string& type);
54   - GModelSpatialRadialDisk(const GSkyDir& dir, const double& radius);
  60 + GModelSpatialRadialDisk(const GSkyDir& dir,
  61 + const double& radius,
  62 + const std::string& coordsys = "CEL");
55 63 explicit GModelSpatialRadialDisk(const GXmlElement& xml);
56 64 GModelSpatialRadialDisk(const GModelSpatialRadialDisk& model);
57 65 virtual ~GModelSpatialRadialDisk(void);
... ...
include/GModelSpatialRadialGauss.hpp
  View file @df50434
1 1 /***************************************************************************
2 2 * GModelSpatialRadialGauss.hpp - Radial Gaussian source model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2020 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2022 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -31,11 +31,15 @@
31 31 #include <string>
32 32 #include "GModelSpatialRadial.hpp"
33 33 #include "GModelPar.hpp"
34   -#include "GSkyDir.hpp"
35   -#include "GSkyRegionCircle.hpp"
36   -#include "GEnergy.hpp"
37   -#include "GTime.hpp"
38   -#include "GXmlElement.hpp"
  34 +
  35 +/* __ Forward declarations _______________________________________________ */
  36 +class GEnergy;
  37 +class GTime;
  38 +class GPhoton;
  39 +class GRan;
  40 +class GSkyDir;
  41 +class GSkyRegion;
  42 +class GXmlElement;
39 43  
40 44  
41 45 /***********************************************************************//**
... ... @@ -52,7 +56,9 @@ public:
52 56 // Constructors and destructors
53 57 GModelSpatialRadialGauss(void);
54 58 GModelSpatialRadialGauss(const bool& dummy, const std::string& type);
55   - GModelSpatialRadialGauss(const GSkyDir& dir, const double& sigma);
  59 + GModelSpatialRadialGauss(const GSkyDir& dir,
  60 + const double& sigma,
  61 + const std::string& coordsys = "CEL");
56 62 explicit GModelSpatialRadialGauss(const GXmlElement& xml);
57 63 GModelSpatialRadialGauss(const GModelSpatialRadialGauss& model);
58 64 virtual ~GModelSpatialRadialGauss(void);
... ...
include/GModelSpatialRadialGeneralGauss.hpp
  View file @df50434
... ... @@ -32,11 +32,15 @@
32 32 #include <string>
33 33 #include "GModelSpatialRadial.hpp"
34 34 #include "GModelPar.hpp"
35   -#include "GSkyDir.hpp"
36   -#include "GSkyRegionCircle.hpp"
37   -#include "GEnergy.hpp"
38   -#include "GTime.hpp"
39   -#include "GXmlElement.hpp"
  35 +
  36 +/* __ Forward declarations _______________________________________________ */
  37 +class GEnergy;
  38 +class GTime;
  39 +class GPhoton;
  40 +class GRan;
  41 +class GSkyDir;
  42 +class GSkyRegion;
  43 +class GXmlElement;
40 44  
41 45  
42 46 /***********************************************************************//**
... ... @@ -52,9 +56,10 @@ class GModelSpatialRadialGeneralGauss : public GModelSpatialRadial {
52 56 public:
53 57 // Constructors and destructors
54 58 GModelSpatialRadialGeneralGauss(void);
55   - GModelSpatialRadialGeneralGauss(const GSkyDir& dir,
56   - const double& radius,
57   - const double& ridx);
  59 + GModelSpatialRadialGeneralGauss(const GSkyDir& dir,
  60 + const double& radius,
  61 + const double& ridx,
  62 + const std::string& coordsys = "CEL");
58 63 explicit GModelSpatialRadialGeneralGauss(const GXmlElement& xml);
59 64 GModelSpatialRadialGeneralGauss(const GModelSpatialRadialGeneralGauss& model);
60 65 virtual ~GModelSpatialRadialGeneralGauss(void);
... ...
include/GModelSpatialRadialRing.hpp
  View file @df50434
1 1 /***************************************************************************
2 2 * GModelSpatialRadialRing.hpp - Radial ring source model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2020 by Pierrick Martin *
  4 + * copyright (C) 2020-2022 by Pierrick Martin *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -31,9 +31,15 @@
31 31 #include <string>
32 32 #include "GModelSpatialRadial.hpp"
33 33 #include "GModelPar.hpp"
34   -#include "GSkyDir.hpp"
35   -#include "GSkyRegionCircle.hpp"
36   -#include "GXmlElement.hpp"
  34 +
  35 +/* __ Forward declarations _______________________________________________ */
  36 +class GEnergy;
  37 +class GTime;
  38 +class GPhoton;
  39 +class GRan;
  40 +class GSkyDir;
  41 +class GSkyRegion;
  42 +class GXmlElement;
37 43  
38 44  
39 45 /**************************************************************************
... ... @@ -50,7 +56,10 @@ class GModelSpatialRadialRing : public GModelSpatialRadial {
50 56 public:
51 57 // Constructors and destructors
52 58 GModelSpatialRadialRing(void);
53   - GModelSpatialRadialRing(const GSkyDir& dir, const double& radius, const double& width);
  59 + GModelSpatialRadialRing(const GSkyDir& dir,
  60 + const double& radius,
  61 + const double& width,
  62 + const std::string& coordsys = "CEL");
54 63 explicit GModelSpatialRadialRing(const GXmlElement& xml);
55 64 GModelSpatialRadialRing(const GModelSpatialRadialRing& model);
56 65 virtual ~GModelSpatialRadialRing(void);
... ...
include/GModelSpatialRadialShell.hpp
  View file @df50434
1 1 /***************************************************************************
2 2 * GModelSpatialRadialShell.hpp - Radial shell source model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2020 by Christoph Deil *
  4 + * copyright (C) 2011-2022 by Christoph Deil *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -31,9 +31,15 @@
31 31 #include <string>
32 32 #include "GModelSpatialRadial.hpp"
33 33 #include "GModelPar.hpp"
34   -#include "GSkyDir.hpp"
35   -#include "GSkyRegionCircle.hpp"
36   -#include "GXmlElement.hpp"
  34 +
  35 +/* __ Forward declarations _______________________________________________ */
  36 +class GEnergy;
  37 +class GTime;
  38 +class GPhoton;
  39 +class GRan;
  40 +class GSkyDir;
  41 +class GSkyRegion;
  42 +class GXmlElement;
37 43  
38 44  
39 45 /**************************************************************************
... ... @@ -61,9 +67,10 @@ public:
61 67 // Constructors and destructors
62 68 GModelSpatialRadialShell(void);
63 69 GModelSpatialRadialShell(const bool& dummy, const std::string& type);
64   - GModelSpatialRadialShell(const GSkyDir& dir,
65   - const double& radius,
66   - const double& width);
  70 + GModelSpatialRadialShell(const GSkyDir& dir,
  71 + const double& radius,
  72 + const double& width,
  73 + const std::string& coordsys = "CEL");
67 74 explicit GModelSpatialRadialShell(const GXmlElement& xml);
68 75 GModelSpatialRadialShell(const GModelSpatialRadialShell& model);
69 76 virtual ~GModelSpatialRadialShell(void);
... ... @@ -110,13 +117,13 @@ protected:
110 117 GModelPar m_width; //!< Shell thickness (deg)
111 118  
112 119 // Cached members used for pre-computations
113   - mutable double m_last_radius; //!< Last shell radius (deg)
114   - mutable double m_last_width; //!< Last shell width (deg)
115   - mutable double m_theta_in; //!< Inner shell radius (rad)
116   - mutable double m_x_in; //!< sin(m_theta_in)^2
117   - mutable double m_theta_out; //!< Outer shell radius (rad)
118   - mutable double m_x_out; //!< sin(m_theta_out)^2
119   - mutable double m_norm; //!< Shell normalization
  120 + mutable double m_last_radius; //!< Last shell radius (deg)
  121 + mutable double m_last_width; //!< Last shell width (deg)
  122 + mutable double m_theta_in; //!< Inner shell radius (rad)
  123 + mutable double m_x_in; //!< sin(m_theta_in)^2
  124 + mutable double m_theta_out; //!< Outer shell radius (rad)
  125 + mutable double m_x_out; //!< sin(m_theta_out)^2
  126 + mutable double m_norm; //!< Shell normalization
120 127 };
121 128  
122 129  
... ...
inst/cta/src/GCTAResponseIrf.cpp
  View file @df50434
1 1 /***************************************************************************
2 2 * GCTAResponseIrf.cpp - CTA instrument response function class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2021 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2022 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -3042,8 +3042,8 @@ double GCTAResponseIrf::nroi_radial(const GModelSky&amp; model,
3042 3042 // given by (rho,omega), into celestial coordinates.
3043 3043 GMatrix ry;
3044 3044 GMatrix rz;
3045   - ry.eulery(spatial->dec() - 90.0);
3046   - rz.eulerz(-spatial->ra());
  3045 + ry.eulery(spatial->dir().dec_deg() - 90.0);
  3046 + rz.eulerz(-spatial->dir().ra_deg());
3047 3047 GMatrix rot = (ry * rz).transpose();
3048 3048  
3049 3049 // Compute position angle of ROI centre with respect to model
... ... @@ -3266,8 +3266,8 @@ double GCTAResponseIrf::nroi_elliptical(const GModelSky&amp; model,
3266 3266 // given by (rho,omega), into celestial coordinates.
3267 3267 GMatrix ry;
3268 3268 GMatrix rz;
3269   - ry.eulery(spatial->dec() - 90.0);
3270   - rz.eulerz(-spatial->ra());
  3269 + ry.eulery(spatial->dir().dec_deg() - 90.0);
  3270 + rz.eulerz(-spatial->dir().ra_deg());
3271 3271 GMatrix rot = (ry * rz).transpose();
3272 3272  
3273 3273 // Compute position angle of ROI centre with respect to model
... ...
inst/cta/src/cta_helpers_response_stacked_vector.cpp
  View file @df50434
1 1 /***************************************************************************
2 2 * CTA helper classes for stacked vector response *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2020 by Juergen Knoedlseder *
  4 + * copyright (C) 2020-2022 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -113,9 +113,18 @@ cta_psf_radial_kerns_delta::cta_psf_radial_kerns_delta(const GCTAResponseCube*
113 113 m_sin_zeta = std::sin(zeta);
114 114 m_cos_theta_max = std::cos(theta_max);
115 115  
116   - // Store references to Right Ascension and Declination parameters
117   - m_par_ra = &((*(const_cast<GModelSpatialRadial*>(model)))["RA"]);
118   - m_par_dec = &((*(const_cast<GModelSpatialRadial*>(model)))["DEC"]);
  116 + // Set coordinate system
  117 + m_par_cel = (model->coordsys() == "CEL");
  118 +
  119 + // Store references to longitude and latitude parameters
  120 + if (m_par_cel) {
  121 + m_par_lon = &((*(const_cast<GModelSpatialRadial*>(model)))["RA"]);
  122 + m_par_lat = &((*(const_cast<GModelSpatialRadial*>(model)))["DEC"]);
  123 + }
  124 + else {
  125 + m_par_lon = &((*(const_cast<GModelSpatialRadial*>(model)))["GLON"]);
  126 + m_par_lat = &((*(const_cast<GModelSpatialRadial*>(model)))["GLAT"]);
  127 + }
119 128  
120 129 // Pre-compute constants needed for gradient computation. The following
121 130 // exceptions need to be handled
... ... @@ -123,10 +132,10 @@ cta_psf_radial_kerns_delta::cta_psf_radial_kerns_delta(const GCTAResponseCube*
123 132 // * m_sin_zeta = 0
124 133 // * denom = 0
125 134 if (grad) {
126   - double alpha_0 = model->ra() * gammalib::deg2rad;
127   - double beta_0 = model->dec() * gammalib::deg2rad;
128   - double alpha_reco = obsDir.ra();
129   - double beta_reco = obsDir.dec();
  135 + double alpha_0 = m_par_cel ? model->dir().ra() : model->dir().l();
  136 + double beta_0 = m_par_cel ? model->dir().dec() : model->dir().b();
  137 + double alpha_reco = m_par_cel ? obsDir.ra() : obsDir.l();
  138 + double beta_reco = m_par_cel ? obsDir.dec() : obsDir.b();
130 139 double sin_beta_0 = std::sin(beta_0);
131 140 double cos_beta_0 = std::cos(beta_0);
132 141 double tan_beta_0 = std::tan(beta_0); // Exception: beta_0 = 90 deg
... ... @@ -396,8 +405,8 @@ GVector cta_psf_radial_kerns_phi::eval(const double&amp; phi)
396 405 // into the respective factor gradients so that the
397 406 // GModelSpatialRadial::eval() method can use the information.
398 407 if (m_outer->m_grad) {
399   - double g_ra = 0.0;
400   - double g_dec = 0.0;
  408 + double g_lon = 0.0;
  409 + double g_lat = 0.0;
401 410 if (theta > 0.0) {
402 411  
403 412 // Compute sin(theta) by avoiding to use a call to sine
... ... @@ -413,19 +422,19 @@ GVector cta_psf_radial_kerns_phi::eval(const double&amp; phi)
413 422 double dtheta_dphi = (m_sin_delta_sin_zeta * sin_phi) *
414 423 norm;
415 424 if (m_outer->m_zeta != 0.0) {
416   - g_ra = dtheta_dzeta * m_outer->m_dzeta_dalpha_0 +
  425 + g_lon = dtheta_dzeta * m_outer->m_dzeta_dalpha_0 +
417 426 dtheta_dphi * m_outer->m_dphi_dalpha_0;
418   - g_dec = dtheta_dzeta * m_outer->m_dzeta_dbeta_0 +
  427 + g_lat = dtheta_dzeta * m_outer->m_dzeta_dbeta_0 +
419 428 dtheta_dphi * m_outer->m_dphi_dbeta_0;
420 429 }
421 430 else {
422   - g_ra = 0.0; //TODO: Set correct value
423   - g_dec = -cos_phi;
  431 + g_lon = 0.0; //TODO: Set correct value
  432 + g_lat = -cos_phi;
424 433 }
425 434 }
426 435 }
427   - m_outer->m_par_ra->factor_gradient(g_ra);
428   - m_outer->m_par_dec->factor_gradient(g_dec);
  436 + m_outer->m_par_lon->factor_gradient(g_lon);
  437 + m_outer->m_par_lat->factor_gradient(g_lat);
429 438 }
430 439  
431 440 // Compute model value and optionally model parameter gradients
... ...
inst/cta/src/cta_helpers_response_stacked_vector.hpp
  View file @df50434
1 1 /***************************************************************************
2 2 * CTA helper classes for stacked vector response *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2020 by Juergen Knoedlseder *
  4 + * copyright (C) 2020-2022 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -63,8 +63,9 @@ public:
63 63 protected:
64 64 const GCTAResponseCube* m_rsp; //!< Response cube
65 65 const GModelSpatialRadial* m_model; //!< Radial model
66   - GModelPar* m_par_ra; //!< Right Ascension parameter
67   - GModelPar* m_par_dec; //!< Declination parameter
  66 + GModelPar* m_par_lon; //!< Longitude parameter
  67 + GModelPar* m_par_lat; //!< Latitude parameter
  68 + bool m_par_cel; //!< Celestial or galactic coordinates
68 69 GSkyDir m_obsDir; //!< Reconstructed event direction
69 70 GEnergies m_srcEngs; //!< True photon energies
70 71 double m_zeta; //!< Distance of model from Psf
... ...
pyext/GModelSpatialRadial.i
  View file @df50434
1 1 /***************************************************************************
2 2 * GModelSpatialRadial.i - Abstract radial spatial model base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2020 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2022 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -68,12 +68,9 @@ public:
68 68 virtual void write(GXmlElement& xml) const;
69 69  
70 70 // Other methods
71   - double ra(void) const;
72   - double dec(void) const;
73   - void ra(const double& ra);
74   - void dec(const double& dec);
75   - GSkyDir dir(void) const;
76   - void dir(const GSkyDir& dir);
  71 + std::string coordsys(void) const;
  72 + const GSkyDir& dir(void) const;
  73 + void dir(const GSkyDir& dir);
77 74 };
78 75  
79 76  
... ...
pyext/GModelSpatialRadialDisk.i
  View file @df50434
1 1 /***************************************************************************
2 2 * GModelSpatialRadialDisk.i - Radial disk source model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2020 by Christoph Deil *
  4 + * copyright (C) 2011-2022 by Christoph Deil *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -42,7 +42,9 @@ class GModelSpatialRadialDisk : public GModelSpatialRadial {
42 42 public:
43 43 // Constructors and destructors
44 44 GModelSpatialRadialDisk(void);
45   - GModelSpatialRadialDisk(const GSkyDir& dir, const double& radius);
  45 + GModelSpatialRadialDisk(const GSkyDir& dir,
  46 + const double& radius,
  47 + const std::string& coordsys = "CEL");
46 48 explicit GModelSpatialRadialDisk(const GXmlElement& xml);
47 49 GModelSpatialRadialDisk(const GModelSpatialRadialDisk& model);
48 50 virtual ~GModelSpatialRadialDisk(void);
... ...
pyext/GModelSpatialRadialGauss.i
  View file @df50434
1 1 /***************************************************************************
2 2 * GModelSpatialRadialGauss.i - Radial Gaussian source model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2020 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2022 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -39,7 +39,9 @@ class GModelSpatialRadialGauss : public GModelSpatialRadial {
39 39 public:
40 40 // Constructors and destructors
41 41 GModelSpatialRadialGauss(void);
42   - GModelSpatialRadialGauss(const GSkyDir& dir, const double& sigma);
  42 + GModelSpatialRadialGauss(const GSkyDir& dir,
  43 + const double& sigma,
  44 + const std::string& coordsys = "CEL");
43 45 explicit GModelSpatialRadialGauss(const GXmlElement& xml);
44 46 GModelSpatialRadialGauss(const GModelSpatialRadialGauss& model);
45 47 virtual ~GModelSpatialRadialGauss(void);
... ...
pyext/GModelSpatialRadialGeneralGauss.i
  View file @df50434
... ... @@ -40,9 +40,10 @@ class GModelSpatialRadialGeneralGauss : public GModelSpatialRadial {
40 40 public:
41 41 // Constructors and destructors
42 42 GModelSpatialRadialGeneralGauss(void);
43   - GModelSpatialRadialGeneralGauss(const GSkyDir& dir,
44   - const double& radius,
45   - const double& ridx);
  43 + GModelSpatialRadialGeneralGauss(const GSkyDir& dir,
  44 + const double& radius,
  45 + const double& ridx,
  46 + const std::string& coordsys = "CEL");
46 47 explicit GModelSpatialRadialGeneralGauss(const GXmlElement& xml);
47 48 GModelSpatialRadialGeneralGauss(const GModelSpatialRadialGeneralGauss& model);
48 49 virtual ~GModelSpatialRadialGeneralGauss(void);
... ...
pyext/GModelSpatialRadialRing.i
  View file @df50434
1 1 /***************************************************************************
2 2 * GModelSpatialRadialRing.i - Radial disk source model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2020 by Pierrick Martin *
  4 + * copyright (C) 2020-2022 by Pierrick Martin *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -42,7 +42,10 @@ class GModelSpatialRadialRing : public GModelSpatialRadial {
42 42 public:
43 43 // Constructors and destructors
44 44 GModelSpatialRadialRing(void);
45   - GModelSpatialRadialRing(const GSkyDir& dir, const double& radius, const double& width);
  45 + GModelSpatialRadialRing(const GSkyDir& dir,
  46 + const double& radius,
  47 + const double& width,
  48 + const std::string& coordsys = "CEL");
46 49 explicit GModelSpatialRadialRing(const GXmlElement& xml);
47 50 GModelSpatialRadialRing(const GModelSpatialRadialRing& model);
48 51 virtual ~GModelSpatialRadialRing(void);
... ...
pyext/GModelSpatialRadialShell.i
  View file @df50434
1 1 /***************************************************************************
2 2 * GModelSpatialRadialShell.i - Radial spatial shell source model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2020 by Christoph Deil *
  4 + * copyright (C) 2011-2022 by Christoph Deil *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -38,9 +38,10 @@ class GModelSpatialRadialShell : public GModelSpatialRadial {
38 38 public:
39 39 // Constructors and destructors
40 40 GModelSpatialRadialShell(void);
41   - GModelSpatialRadialShell(const GSkyDir& dir,
42   - const double& radius,
43   - const double& width);
  41 + GModelSpatialRadialShell(const GSkyDir& dir,
  42 + const double& radius,
  43 + const double& width,
  44 + const std::string& coordsys = "CEL");
44 45 explicit GModelSpatialRadialShell(const GXmlElement& xml);
45 46 GModelSpatialRadialShell(const GModelSpatialRadialShell& model);
46 47 virtual ~GModelSpatialRadialShell(void);
... ...
src/model/GModelSpatialPointSource.cpp
  View file @df50434
... ... @@ -115,8 +115,8 @@ GModelSpatialPointSource::GModelSpatialPointSource(const bool&amp; dummy,
115 115 * Invalid @p coordsys argument specified.
116 116 *
117 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.
  118 + * @p coordsys parameter specifies whether the sky direction should be
  119 + * interpreted in the celestial or Galactic coordinate system.
120 120 ***************************************************************************/
121 121 GModelSpatialPointSource::GModelSpatialPointSource(const GSkyDir& dir,
122 122 const std::string& coordsys) :
... ...
src/model/GModelSpatialRadial.cpp
  View file @df50434
1 1 /***************************************************************************
2 2 * GModelSpatialRadial.cpp - Abstract radial spatial model base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2018 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2022 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -199,16 +199,16 @@ double GModelSpatialRadial::eval(const GPhoton&amp; photon,
199 199 * XML element in the following format
200 200 *
201 201 * <spatialModel type="...">
202   - * <parameter name="RA" scale="1" value="83.63" min="-360" max="360" free="1"/>
203   - * <parameter name="DEC" scale="1" value="22.01" min="-90" max="90" free="1"/>
  202 + * <parameter name="RA" scale="1" value="83.6331" min="-360" max="360" free="0" />
  203 + * <parameter name="DEC" scale="1" value="22.0145" min="-90" max="90" free="0" />
204 204 * ...
205 205 * </spatialModel>
206 206 *
207 207 * or
208 208 *
209 209 * <spatialModel type="...">
210   - * <parameter name="GLON" scale="1" value="83.63" min="-360" max="360" free="1"/>
211   - * <parameter name="GLAT" scale="1" value="22.01" min="-90" max="90" free="1"/>
  210 + * <parameter name="GLON" scale="1" value="184.5575" min="-360" max="360" free="0" />
  211 + * <parameter name="GLAT" scale="1" value="-5.7843" min="-90" max="90" free="0" />
212 212 * ...
213 213 * </spatialModel>
214 214 *
... ... @@ -223,8 +223,8 @@ void GModelSpatialRadial::read(const GXmlElement&amp; xml)
223 223 const GXmlElement* dec = gammalib::xml_get_par(G_READ, xml, "DEC");
224 224  
225 225 // Read parameters
226   - m_ra.read(*ra);
227   - m_dec.read(*dec);
  226 + m_lon.read(*ra);
  227 + m_lat.read(*dec);
228 228  
229 229 }
230 230  
... ... @@ -236,18 +236,8 @@ void GModelSpatialRadial::read(const GXmlElement&amp; xml)
236 236 const GXmlElement* glat = gammalib::xml_get_par(G_READ, xml, "GLAT");
237 237  
238 238 // Read parameters
239   - m_ra.read(*glon);
240   - m_dec.read(*glat);
241   -
242   - // Convert into RA/DEC
243   - GSkyDir dir;
244   - dir.lb_deg(ra(), dec()),
245   - m_ra.value(dir.ra_deg());
246   - m_dec.value(dir.dec_deg());
247   -
248   - // Set names to RA/DEC
249   - m_ra.name("RA");
250   - m_dec.name("DEC");
  239 + m_lon.read(*glon);
  240 + m_lat.read(*glat);
251 241  
252 242 }
253 243  
... ... @@ -261,24 +251,23 @@ void GModelSpatialRadial::read(const GXmlElement&amp; xml)
261 251 *
262 252 * @param[in] xml XML element into which model information is written.
263 253 *
264   - * @exception GException::model_invalid_spatial
265   - * Existing XML element is not of type 'GaussFunction'
266   - * @exception GException::model_invalid_parnum
267   - * Invalid number of model parameters found in XML element.
268   - * @exception GException::model_invalid_parnames
269   - * Invalid model parameter names found in XML element.
270   - *
271   - * Writes the radial source location and position angle information into an
272   - * XML element in the following format
  254 + * Depending on the coordinate system, write the radial source information
  255 + * into an XML element with the following format
273 256 *
274 257 * <spatialModel type="...">
275   - * <parameter name="RA" scale="1" value="83.63" min="-360" max="360" free="1"/>
276   - * <parameter name="DEC" scale="1" value="22.01" min="-90" max="90" free="1"/>
  258 + * <parameter name="RA" scale="1" value="83.6331" min="-360" max="360" free="0" />
  259 + * <parameter name="DEC" scale="1" value="22.0145" min="-90" max="90" free="0" />
  260 + * ...
  261 + * </spatialModel>
  262 + *
  263 + * or
  264 + *
  265 + * <spatialModel type="PointSource">
  266 + * <parameter name="GLON" scale="1" value="184.5575" min="-360" max="360" free="0" />
  267 + * <parameter name="GLAT" scale="1" value="-5.7843" min="-90" max="90" free="0" />
277 268 * ...
278 269 * </spatialModel>
279 270 *
280   - * @todo The case that an existing spatial XML element with "GLON" and "GLAT"
281   - * as coordinates is not supported.
282 271 ***************************************************************************/
283 272 void GModelSpatialRadial::write(GXmlElement& xml) const
284 273 {
... ... @@ -286,12 +275,12 @@ void GModelSpatialRadial::write(GXmlElement&amp; xml) const
286 275 gammalib::xml_check_type(G_WRITE, xml, type());
287 276  
288 277 // Get or create parameters
289   - GXmlElement* ra = gammalib::xml_need_par(G_WRITE, xml, m_ra.name());
290   - GXmlElement* dec = gammalib::xml_need_par(G_WRITE, xml, m_dec.name());
  278 + GXmlElement* lon = gammalib::xml_need_par(G_WRITE, xml, m_lon.name());
  279 + GXmlElement* lat = gammalib::xml_need_par(G_WRITE, xml, m_lat.name());
291 280  
292 281 // Write parameters
293   - m_ra.write(*ra);
294   - m_dec.write(*dec);
  282 + m_lon.write(*lon);
  283 + m_lat.write(*lat);
295 284  
296 285 // Return
297 286 return;
... ... @@ -300,28 +289,58 @@ void GModelSpatialRadial::write(GXmlElement&amp; xml) const
300 289  
301 290 /***********************************************************************//**
302 291 * @brief Return position of radial spatial model
  292 + *
  293 + * @return Radial spatial model sky direction.
  294 + *
  295 + * Returns the sky direction of the radial spatial model.
303 296 ***************************************************************************/
304   -GSkyDir GModelSpatialRadial::dir(void) const
  297 +const GSkyDir& GModelSpatialRadial::dir(void) const
305 298 {
306   - // Allocate sky direction
307   - GSkyDir srcDir;
308   -
309   - // Set sky direction
310   - srcDir.radec_deg(ra(), dec());
311   -
312   - // Return direction
313   - return srcDir;
  299 + // Get longitude and latitude values
  300 + double lon = m_lon.value();
  301 + double lat = m_lat.value();
  302 +
  303 + // If longitude or latitude values have changed then update sky
  304 + // direction cache
  305 + if ((lon != m_last_lon) || (lat != m_last_lat)) {
  306 +
  307 + // Update last values
  308 + m_last_lon = lon;
  309 + m_last_lat = lat;
  310 +
  311 + // Update sky direction dependent on model coordinate system
  312 + if (is_celestial()) {
  313 + m_dir.radec_deg(m_last_lon, m_last_lat);
  314 + }
  315 + else {
  316 + m_dir.lb_deg(m_last_lon, m_last_lat);
  317 + }
  318 +
  319 + } // endif: update of sky direction cache required
  320 +
  321 + // Return sky direction
  322 + return (m_dir);
314 323 }
315 324  
316 325  
317 326 /***********************************************************************//**
318 327 * @brief Set position of radial spatial model
  328 + *
  329 + * @param[in] dir Sky direction of radial spatial model.
  330 + *
  331 + * Sets the sky direction of the radial spatial model.
319 332 ***************************************************************************/
320 333 void GModelSpatialRadial::dir(const GSkyDir& dir)
321 334 {
322   - // Assign Right Ascension and Declination
323   - m_ra.value(dir.ra_deg());
324   - m_dec.value(dir.dec_deg());
  335 + // Assign sky direction depending on the model coordinate system
  336 + if (is_celestial()) {
  337 + m_lon.value(dir.ra_deg());
  338 + m_lat.value(dir.dec_deg());
  339 + }
  340 + else {
  341 + m_lon.value(dir.l_deg());
  342 + m_lat.value(dir.b_deg());
  343 + }
325 344  
326 345 // Return
327 346 return;
... ... @@ -340,27 +359,32 @@ void GModelSpatialRadial::dir(const GSkyDir&amp; dir)
340 359 void GModelSpatialRadial::init_members(void)
341 360 {
342 361 // Initialise Right Ascension
343   - m_ra.clear();
344   - m_ra.name("RA");
345   - m_ra.unit("deg");
346   - m_ra.fix();
347   - m_ra.scale(1.0);
348   - m_ra.gradient(0.0);
349   - m_ra.has_grad(false);
  362 + m_lon.clear();
  363 + m_lon.name("RA");
  364 + m_lon.unit("deg");
  365 + m_lon.fix();
  366 + m_lon.scale(1.0);
  367 + m_lon.gradient(0.0);
  368 + m_lon.has_grad(false);
350 369  
351 370 // Initialise Declination
352   - m_dec.clear();
353   - m_dec.name("DEC");
354   - m_dec.unit("deg");
355   - m_dec.fix();
356   - m_dec.scale(1.0);
357   - m_dec.gradient(0.0);
358   - m_dec.has_grad(false);
  371 + m_lat.clear();
  372 + m_lat.name("DEC");
  373 + m_lat.unit("deg");
  374 + m_lat.fix();
  375 + m_lat.scale(1.0);
  376 + m_lat.gradient(0.0);
  377 + m_lat.has_grad(false);
359 378  
360 379 // Set parameter pointer(s)
361 380 m_pars.clear();
362   - m_pars.push_back(&m_ra);
363   - m_pars.push_back(&m_dec);
  381 + m_pars.push_back(&m_lon);
  382 + m_pars.push_back(&m_lat);
  383 +
  384 + // Initialise cache
  385 + m_dir.clear();
  386 + m_last_lon = -9999.0;
  387 + m_last_lat = -9999.0;
364 388  
365 389 // Return
366 390 return;
... ... @@ -375,13 +399,18 @@ void GModelSpatialRadial::init_members(void)
375 399 void GModelSpatialRadial::copy_members(const GModelSpatialRadial& model)
376 400 {
377 401 // Copy members
378   - m_ra = model.m_ra;
379   - m_dec = model.m_dec;
  402 + m_lon = model.m_lon;
  403 + m_lat = model.m_lat;
380 404  
381 405 // Set parameter pointer(s)
382 406 m_pars.clear();
383   - m_pars.push_back(&m_ra);
384   - m_pars.push_back(&m_dec);
  407 + m_pars.push_back(&m_lon);
  408 + m_pars.push_back(&m_lat);
  409 +
  410 + // Copy cache
  411 + m_dir = model.m_dir;
  412 + m_last_lon = model.m_last_lon;
  413 + m_last_lat = model.m_last_lat;
385 414  
386 415 // Return
387 416 return;
... ...
src/model/GModelSpatialRadialDisk.cpp
  View file @df50434
1 1 /***************************************************************************
2 2 * GModelSpatialRadialDisk.cpp - Radial disk source model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2021 by Christoph Deil *
  4 + * copyright (C) 2011-2022 by Christoph Deil *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -45,6 +45,8 @@ const GModelSpatialRegistry g_radial_disk_legacy_registry(&amp;g_radial_disk_legac
45 45 #endif
46 46  
47 47 /* __ Method name definitions ____________________________________________ */
  48 +#define G_CONSTRUCTOR "GModelSpatialRadialDisk::GModelSpatialRadialDisk("\
  49 + "GSkyDir&, double&, std::string&)"
48 50 #define G_READ "GModelSpatialRadialDisk::read(GXmlElement&)"
49 51 #define G_WRITE "GModelSpatialRadialDisk::write(GXmlElement&)"
50 52  
... ... @@ -104,17 +106,42 @@ GModelSpatialRadialDisk::GModelSpatialRadialDisk(const bool&amp; dummy,
104 106 *
105 107 * @param[in] dir Sky position of disk centre.
106 108 * @param[in] radius Disk radius (degrees).
  109 + * @param[in] coordsys Coordinate system (either "CEL" or "GAL")
  110 + *
  111 + * @exception GException::invalid_argument
  112 + * Invalid @p coordsys argument specified.
107 113 *
108 114 * Constructs radial disk model from the sky position of the disk centre
109   - * (@p dir) and the disk @p radius in degrees.
  115 + * (@p dir) and the disk @p radius in degrees. The @p coordsys parameter
  116 + * specifies whether the sky direction should be interpreted in the
  117 + * celestial or Galactic coordinate system.
110 118 ***************************************************************************/
111   -GModelSpatialRadialDisk::GModelSpatialRadialDisk(const GSkyDir& dir,
112   - const double& radius) :
  119 +GModelSpatialRadialDisk::GModelSpatialRadialDisk(const GSkyDir& dir,
  120 + const double& radius,
  121 + const std::string& coordsys) :
113 122 GModelSpatialRadial()
114 123 {
  124 + // Throw an exception if the coordinate system is invalid
  125 + if ((coordsys != "CEL") && (coordsys != "GAL")) {
  126 + std::string msg = "Invalid coordinate system \""+coordsys+"\" "
  127 + "specified. Please specify either \"CEL\" or "
  128 + "\"GAL\".";
  129 + throw GException::invalid_argument(G_CONSTRUCTOR, msg);
  130 + }
  131 +
115 132 // Initialise members
116 133 init_members();
117 134  
  135 + // Set parameter names
  136 + if (coordsys == "CEL") {
  137 + m_lon.name("RA");
  138 + m_lat.name("DEC");
  139 + }
  140 + else {
  141 + m_lon.name("GLON");
  142 + m_lat.name("GLAT");
  143 + }
  144 +
118 145 // Assign parameters
119 146 this->dir(dir);
120 147 this->radius(radius);
... ... @@ -382,11 +409,6 @@ double GModelSpatialRadialDisk::theta_max(void) const
382 409 *
383 410 * @param[in] xml XML element.
384 411 *
385   - * @exception GException::model_invalid_parnum
386   - * Invalid number of model parameters found in XML element.
387   - * @exception GException::model_invalid_parnames
388   - * Invalid model parameter names found in XML element.
389   - *
390 412 * Reads the radial disk model information from an XML element. The XML
391 413 * element shall have either the format
392 414 *
... ... @@ -611,7 +633,7 @@ void GModelSpatialRadialDisk::update() const
611 633 void GModelSpatialRadialDisk::set_region(void) const
612 634 {
613 635 // Set sky region circle
614   - GSkyRegionCircle region(m_ra.value(), m_dec.value(), m_radius.value());
  636 + GSkyRegionCircle region(dir(), m_radius.value());
615 637  
616 638 // Set region (circumvent const correctness)
617 639 const_cast<GModelSpatialRadialDisk*>(this)->m_region = region;
... ...
src/model/GModelSpatialRadialGauss.cpp
  View file @df50434
1 1 /***************************************************************************
2 2 * GModelSpatialRadialGauss.cpp - Radial Gaussian source model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2020 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2022 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -45,6 +45,8 @@ const GModelSpatialRegistry g_radial_gauss_legacy_registry(&amp;g_radial_gauss_le
45 45 #endif
46 46  
47 47 /* __ Method name definitions ____________________________________________ */
  48 +#define G_CONSTRUCTOR "GModelSpatialRadialGauss::GModelSpatialRadialGauss("\
  49 + "GSkyDir&, double&, std::string&)"
48 50 #define G_EVAL "GModelSpatialRadialGauss::eval(double&, GEnergy&, GTime&, "\
49 51 "bool&)"
50 52 #define G_READ "GModelSpatialRadialGauss::read(GXmlElement&)"
... ... @@ -106,17 +108,42 @@ GModelSpatialRadialGauss::GModelSpatialRadialGauss(const bool&amp; dummy,
106 108 *
107 109 * @param[in] dir Sky position of Gaussian.
108 110 * @param[in] sigma Width of Gaussian (in degrees).
  111 + * @param[in] coordsys Coordinate system (either "CEL" or "GAL")
  112 + *
  113 + * @exception GException::invalid_argument
  114 + * Invalid @p coordsys argument specified.
109 115 *
110 116 * Constructs a Gaussian spatial model using a sky direction (@p dir) and
111   - * a Gaussian width parameter @p sigma in degrees.
  117 + * a Gaussian width parameter @p sigma in degrees. The @p coordsys parameter
  118 + * specifies whether the sky direction should be interpreted in the
  119 + * celestial or Galactic coordinate system.
112 120 ***************************************************************************/
113   -GModelSpatialRadialGauss::GModelSpatialRadialGauss(const GSkyDir& dir,
114   - const double& sigma) :
  121 +GModelSpatialRadialGauss::GModelSpatialRadialGauss(const GSkyDir& dir,
  122 + const double& sigma,
  123 + const std::string& coordsys) :
115 124 GModelSpatialRadial()
116 125 {
  126 + // Throw an exception if the coordinate system is invalid
  127 + if ((coordsys != "CEL") && (coordsys != "GAL")) {
  128 + std::string msg = "Invalid coordinate system \""+coordsys+"\" "
  129 + "specified. Please specify either \"CEL\" or "
  130 + "\"GAL\".";
  131 + throw GException::invalid_argument(G_CONSTRUCTOR, msg);
  132 + }
  133 +
117 134 // Initialise members
118 135 init_members();
119 136  
  137 + // Set parameter names
  138 + if (coordsys == "CEL") {
  139 + m_lon.name("RA");
  140 + m_lat.name("DEC");
  141 + }
  142 + else {
  143 + m_lon.name("GLON");
  144 + m_lat.name("GLAT");
  145 + }
  146 +
120 147 // Assign direction and sigma
121 148 this->dir(dir);
122 149 this->sigma(sigma);
... ... @@ -312,7 +339,7 @@ GModelSpatialRadialGauss* GModelSpatialRadialGauss::clone(void) const
312 339 * centre relies on the computation of the partial derivatives
313 340 * \f$\frac{\partial \theta}{\partial \alpha_0}\f$ and
314 341 * \f$\frac{\partial \theta}{\partial \beta_0}\f$ which are assumed to
315   - * be stored in `m_ra.factor_gradient()` and `m_dec.factor_gradient()`
  342 + * be stored in `m_lon.factor_gradient()` and `m_lat.factor_gradient()`
316 343 * when the method is entered. No check is performed whether these values
317 344 * are actually available and valid.
318 345 *
... ... @@ -336,22 +363,23 @@ double GModelSpatialRadialGauss::eval(const double&amp; theta,
336 363 if (gradients) {
337 364  
338 365 // Evaluate position gradients. Note that this requires that the
339   - // partial derivatives of theta with respect to Right Ascension
340   - // and Declination are preset in the m_ra.factor_gradient() and
341   - // m_dec.factor_gradient() members.
342   - double g_ra = 0.0;
343   - double g_dec = 0.0;
344   - if (m_ra.is_free() || m_dec.is_free()) {
  366 + // partial derivatives of theta with respect to longitude and
  367 + // latitude are preset in the m_lon.factor_gradient() and
  368 + // m_lat.factor_gradient() members. An example where these values
  369 + // are set is the cta_psf_radial_kerns_phi::eval() method.
  370 + double g_lon = 0.0;
  371 + double g_lat = 0.0;
  372 + if (m_lon.is_free() || m_lat.is_free()) {
345 373 double g_theta = -theta * m_g_theta_norm * value;
346   - if (m_ra.is_free()) {
347   - g_ra = g_theta * m_ra.factor_gradient() * m_ra.scale();
  374 + if (m_lon.is_free()) {
  375 + g_lon = g_theta * m_lon.factor_gradient() * m_lon.scale();
348 376 }
349   - if (m_dec.is_free()) {
350   - g_dec = g_theta * m_dec.factor_gradient() * m_dec.scale();
  377 + if (m_lat.is_free()) {
  378 + g_lat = g_theta * m_lat.factor_gradient() * m_lat.scale();
351 379 }
352 380 }
353   - m_ra.factor_gradient(g_ra);
354   - m_dec.factor_gradient(g_dec);
  381 + m_lon.factor_gradient(g_lon);
  382 + m_lat.factor_gradient(g_lat);
355 383  
356 384 // Evaluate sigma gradient
357 385 double g_sigma = 0.0;
... ... @@ -578,8 +606,8 @@ void GModelSpatialRadialGauss::init_members(void)
578 606 m_sigma.gradient(0.0);
579 607  
580 608 // Signal that this model provides analytical parameter gradients
581   - m_ra.has_grad(true);
582   - m_dec.has_grad(true);
  609 + m_lon.has_grad(true);
  610 + m_lat.has_grad(true);
583 611 m_sigma.has_grad(true);
584 612  
585 613 // Set parameter pointer(s)
... ... @@ -678,7 +706,7 @@ void GModelSpatialRadialGauss::update(const bool&amp; gradients) const
678 706 void GModelSpatialRadialGauss::set_region(void) const
679 707 {
680 708 // Set sky region circle (5 times Gaussian sigma)
681   - GSkyRegionCircle region(m_ra.value(), m_dec.value(), m_sigma.value() * 5.0);
  709 + GSkyRegionCircle region(dir(), m_sigma.value() * 5.0);
682 710  
683 711 // Set region (circumvent const correctness)
684 712 const_cast<GModelSpatialRadialGauss*>(this)->m_region = region;
... ...
src/model/GModelSpatialRadialGeneralGauss.cpp
  View file @df50434
... ... @@ -42,6 +42,8 @@ const GModelSpatialRadialGeneralGauss g_radial_general_gauss_seed;
42 42 const GModelSpatialRegistry g_radial_general_gauss_registry(&g_radial_general_gauss_seed);
43 43  
44 44 /* __ Method name definitions ____________________________________________ */
  45 +#define G_CONSTRUCTOR "GModelSpatialRadialGeneralGauss::"\
  46 + "GModelSpatialRadialGeneralGauss(GSkyDir&, double&, double&, std::string&)"
45 47 #define G_EVAL "GModelSpatialRadialGeneralGauss::eval(double&, GEnergy&, "\
46 48 "GTime&, bool&)"
47 49 #define G_READ "GModelSpatialRadialGeneralGauss::read(GXmlElement&)"
... ... @@ -83,19 +85,36 @@ GModelSpatialRadialGeneralGauss::GModelSpatialRadialGeneralGauss(void) :
83 85 * @param[in] dir Sky position of Gaussian.
84 86 * @param[in] radius Width of generalised Gaussian (deg).
85 87 * @param[in] ridx Reciprocal of exponential index of radial profile.
  88 + * @param[in] coordsys Coordinate system (either "CEL" or "GAL")
  89 + *
  90 + * @exception GException::invalid_argument
  91 + * Invalid @p coordsys argument specified.
86 92 *
87 93 * Constructs a generalised Gaussian spatial model using a sky direction
88 94 * (@p dir), a Gaussian width parameter @p radius in degrees and a reciprocal
89   - * of the exponential index @p ridx.
  95 + * of the exponential index @p ridx. The @p coordsys parameter specifies
  96 + * whether the sky direction should be interpreted in the celestial or
  97 + * Galactic coordinate system.
90 98 ***************************************************************************/
91   -GModelSpatialRadialGeneralGauss::GModelSpatialRadialGeneralGauss(const GSkyDir& dir,
92   - const double& radius,
93   - const double& ridx) :
  99 +GModelSpatialRadialGeneralGauss::GModelSpatialRadialGeneralGauss(const GSkyDir& dir,
  100 + const double& radius,
  101 + const double& ridx,
  102 + const std::string& coordsys) :
94 103 GModelSpatialRadial()
95 104 {
96 105 // Initialise members
97 106 init_members();
98 107  
  108 + // Set parameter names
  109 + if (coordsys == "CEL") {
  110 + m_lon.name("RA");
  111 + m_lat.name("DEC");
  112 + }
  113 + else {
  114 + m_lon.name("GLON");
  115 + m_lat.name("GLAT");
  116 + }
  117 +
99 118 // Assign direction and radius
100 119 this->dir(dir);
101 120 this->radius(radius);
... ... @@ -634,7 +653,7 @@ void GModelSpatialRadialGeneralGauss::update() const
634 653 void GModelSpatialRadialGeneralGauss::set_region(void) const
635 654 {
636 655 // Set sky region circle (5 times radius)
637   - GSkyRegionCircle region(m_ra.value(), m_dec.value(), m_radius.value() * 5.0);
  656 + GSkyRegionCircle region(dir(), m_radius.value() * 5.0);
638 657  
639 658 // Set region (circumvent const correctness)
640 659 const_cast<GModelSpatialRadialGeneralGauss*>(this)->m_region = region;
... ...
src/model/GModelSpatialRadialProfile.cpp
  View file @df50434
1 1 /***************************************************************************
2 2 * GModelSpatialRadialProfile.cpp - Radial profile source model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2016-2020 by Juergen Knoedlseder *
  4 + * copyright (C) 2016-2022 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -441,7 +441,7 @@ int GModelSpatialRadialProfile::cache_index(void) const
441 441 void GModelSpatialRadialProfile::set_region(void) const
442 442 {
443 443 // Set sky region circle
444   - GSkyRegionCircle region(m_ra.value(), m_dec.value(), theta_max()*gammalib::rad2deg);
  444 + GSkyRegionCircle region(dir(), theta_max()*gammalib::rad2deg);
445 445  
446 446 // Set region (circumvent const correctness)
447 447 const_cast<GModelSpatialRadialProfile*>(this)->m_region = region;
... ...
src/model/GModelSpatialRadialRing.cpp
  View file @df50434
1 1 /***************************************************************************
2 2 * GModelSpatialRadialRing.cpp - Radial ring source model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2020-2021 by Pierrick Martin *
  4 + * copyright (C) 2020-2022 by Pierrick Martin *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -41,6 +41,8 @@ const GModelSpatialRadialRing g_radial_ring_seed;
41 41 const GModelSpatialRegistry g_radial_ring_registry(&g_radial_ring_seed);
42 42  
43 43 /* __ Method name definitions ____________________________________________ */
  44 +#define G_CONSTRUCTOR "GModelSpatialRadialRing::GModelSpatialRadialRing("\
  45 + "GSkyDir&, double&, double&, std::string&)"
44 46 #define G_READ "GModelSpatialRadialRing::read(GXmlElement&)"
45 47 #define G_WRITE "GModelSpatialRadialRing::write(GXmlElement&)"
46 48  
... ... @@ -78,19 +80,43 @@ GModelSpatialRadialRing::GModelSpatialRadialRing(void) : GModelSpatialRadial()
78 80 * @param[in] dir Sky position of ring centre.
79 81 * @param[in] radius Ring inner radius (degrees).
80 82 * @param[in] width Ring width (degrees).
  83 + * @param[in] coordsys Coordinate system (either "CEL" or "GAL")
  84 + *
  85 + * @exception GException::invalid_argument
  86 + * Invalid @p coordsys argument specified.
81 87 *
82 88 * Constructs radial ring model from the sky position of the ring centre
83 89 * (@p dir) and the ring inner radius (@p radius) and width (@p width) in
84   - * degrees.
  90 + * degrees. The @p coordsys parameter specifies whether the sky direction
  91 + * should be interpreted in the celestial or Galactic coordinate system.
85 92 ***************************************************************************/
86   -GModelSpatialRadialRing::GModelSpatialRadialRing(const GSkyDir& dir,
87   - const double& radius,
88   - const double& width) :
  93 +GModelSpatialRadialRing::GModelSpatialRadialRing(const GSkyDir& dir,
  94 + const double& radius,
  95 + const double& width,
  96 + const std::string& coordsys) :
89 97 GModelSpatialRadial()
90 98 {
  99 + // Throw an exception if the coordinate system is invalid
  100 + if ((coordsys != "CEL") && (coordsys != "GAL")) {
  101 + std::string msg = "Invalid coordinate system \""+coordsys+"\" "
  102 + "specified. Please specify either \"CEL\" or "
  103 + "\"GAL\".";
  104 + throw GException::invalid_argument(G_CONSTRUCTOR, msg);
  105 + }
  106 +
91 107 // Initialise members
92 108 init_members();
93 109  
  110 + // Set parameter names
  111 + if (coordsys == "CEL") {
  112 + m_lon.name("RA");
  113 + m_lat.name("DEC");
  114 + }
  115 + else {
  116 + m_lon.name("GLON");
  117 + m_lat.name("GLAT");
  118 + }
  119 +
94 120 // Assign center location, radius and width
95 121 this->dir(dir);
96 122 this->radius(radius);
... ... @@ -635,7 +661,7 @@ void GModelSpatialRadialRing::set_region(void) const
635 661 {
636 662 // Set sky region circle (maximum Gaussian sigma times a scaling
637 663 // factor (actually 3))
638   - GSkyRegionCircle region(m_ra.value(), m_dec.value(), radius()+width());
  664 + GSkyRegionCircle region(dir(), radius()+width());
639 665  
640 666 // Set region (circumvent const correctness)
641 667 const_cast<GModelSpatialRadialRing*>(this)->m_region = region;
... ...
src/model/GModelSpatialRadialShell.cpp
  View file @df50434
1 1 /***************************************************************************
2 2 * GModelSpatialRadialShell.cpp - Radial shell source model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2021 by Christoph Deil *
  4 + * copyright (C) 2011-2022 by Christoph Deil *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -45,6 +45,8 @@ const GModelSpatialRegistry g_radial_shell_legacy_registry(&amp;g_radial_shell_le
45 45 #endif
46 46  
47 47 /* __ Method name definitions ____________________________________________ */
  48 +#define G_CONSTRUCTOR "GModelSpatialRadialShell::GModelSpatialRadialShell("\
  49 + "GSkyDir&, double&, double&, std::string&)"
48 50 #define G_MC "GModelSpatialRadialShell::mc(GEnergy&, GTime&, GRan&)"
49 51 #define G_READ "GModelSpatialRadialShell::read(GXmlElement&)"
50 52 #define G_WRITE "GModelSpatialRadialShell::write(GXmlElement&)"
... ... @@ -107,18 +109,43 @@ GModelSpatialRadialShell::GModelSpatialRadialShell(const bool&amp; dummy,
107 109 * @param[in] dir Sky position of shell centre.
108 110 * @param[in] radius Inner shell radius (degrees).
109 111 * @param[in] width Shell width (degrees).
  112 + * @param[in] coordsys Coordinate system (either "CEL" or "GAL")
  113 + *
  114 + * @exception GException::invalid_argument
  115 + * Invalid @p coordsys argument specified.
110 116 *
111 117 * Constructs the shell model from the shell centre (@p dir), the inner
112   - * shell @p radius, and the shell @p width.
  118 + * shell @p radius, and the shell @p width. The @p coordsys parameter
  119 + * specifies whether the sky direction should be interpreted in the celestial
  120 + * or Galactic coordinate system.
113 121 ***************************************************************************/
114   -GModelSpatialRadialShell::GModelSpatialRadialShell(const GSkyDir& dir,
115   - const double& radius,
116   - const double& width) :
  122 +GModelSpatialRadialShell::GModelSpatialRadialShell(const GSkyDir& dir,
  123 + const double& radius,
  124 + const double& width,
  125 + const std::string& coordsys) :
117 126 GModelSpatialRadial()
118 127 {
  128 + // Throw an exception if the coordinate system is invalid
  129 + if ((coordsys != "CEL") && (coordsys != "GAL")) {
  130 + std::string msg = "Invalid coordinate system \""+coordsys+"\" "
  131 + "specified. Please specify either \"CEL\" or "
  132 + "\"GAL\".";
  133 + throw GException::invalid_argument(G_CONSTRUCTOR, msg);
  134 + }
  135 +
119 136 // Initialise members
120 137 init_members();
121 138  
  139 + // Set parameter names
  140 + if (coordsys == "CEL") {
  141 + m_lon.name("RA");
  142 + m_lat.name("DEC");
  143 + }
  144 + else {
  145 + m_lon.name("GLON");
  146 + m_lat.name("GLAT");
  147 + }
  148 +
122 149 // Assign parameters
123 150 this->dir(dir);
124 151 this->radius(radius);
... ... @@ -758,7 +785,7 @@ void GModelSpatialRadialShell::set_region(void) const
758 785 {
759 786 // Set sky region circle (maximum Gaussian sigma times a scaling
760 787 // factor (actually 3))
761   - GSkyRegionCircle region(m_ra.value(), m_dec.value(), m_radius.value() + m_width.value());
  788 + GSkyRegionCircle region(dir(), m_radius.value() + m_width.value());
762 789  
763 790 // Set region (circumvent const correctness)
764 791 const_cast<GModelSpatialRadialShell*>(this)->m_region = region;
... ...
test/test_GModel.cpp
  View file @df50434
... ... @@ -1210,9 +1210,19 @@ void TestGModel::test_radial_disk(void)
1210 1210 GSkyDir dir1;
1211 1211 dir1.radec_deg(83.6331, +22.0145);
1212 1212 GModelSpatialRadialDisk model2(dir1, 3.0);
1213   - test_value(model2.ra(), 83.6331);
1214   - test_value(model2.dec(), 22.0145);
1215   - test_value(model2.radius(), 3.0);
  1213 + test_value(model2.dir().ra_deg(), 83.6331, "Test Right Ascension for celestial sky direction constructor");
  1214 + test_value(model2.dir().dec_deg(), 22.0145, "Test Declination for celestial sky direction constructor");
  1215 + test_value(model2.dir().l_deg(), 184.557460, "Test Galactic longitude for celestial sky direction constructor");
  1216 + test_value(model2.dir().b_deg(), -5.784346, "Test Galactic latitude for celestial sky direction constructor");
  1217 + test_value(model2.radius(), 3.0, "Test radius for celestial sky direction constructor");
  1218 +
  1219 + // Test Galactic sky direction constructor
  1220 + GModelSpatialRadialDisk model2g(dir1, 3.0, "GAL");
  1221 + test_value(model2g.dir().ra_deg(), 83.6331, "Test Right Ascension for Galactic sky direction constructor");
  1222 + test_value(model2g.dir().dec_deg(), +22.0145, "Test Declination for Galactic sky direction constructor");
  1223 + test_value(model2g.dir().l_deg(), 184.557460, "Test Galactic longitude for Galactic sky direction constructor");
  1224 + test_value(model2g.dir().b_deg(), -5.784346, "Test Galactic latitude for Galactic sky direction constructor");
  1225 + test_value(model2g.radius(), 3.0, "Test radius for Galactic sky direction constructor");
1216 1226  
1217 1227 // Test XML constructor
1218 1228 GXml xml(m_xml_model_radial_disk);
... ... @@ -1220,18 +1230,10 @@ void TestGModel::test_radial_disk(void)
1220 1230 GModelSpatialRadialDisk model3(*element);
1221 1231 test_value(model3.size(), 3);
1222 1232 test_value(model3.type(), "RadialDisk");
1223   - test_value(model3.ra(), 83.6331);
1224   - test_value(model3.dec(), 22.0145);
  1233 + test_value(model3.dir().ra_deg(), 83.6331, "Test Right Ascension for celestial XML constructor");
  1234 + test_value(model3.dir().dec_deg(), 22.0145, "Test Declination for celestial XML constructor");
1225 1235 test_value(model3.radius(), 0.45);
1226 1236  
1227   - // Test ra method
1228   - model3.ra(100.0);
1229   - test_value(model3.ra(), 100.0);
1230   -
1231   - // Test dec method
1232   - model3.dec(10.0);
1233   - test_value(model3.dec(), 10.0);
1234   -
1235 1237 // Test dir method
1236 1238 GSkyDir dir2;
1237 1239 dir2.radec_deg(83.6331, +22.0145);
... ... @@ -1278,11 +1280,22 @@ void TestGModel::test_radial_ring(void)
1278 1280 // Test value constructor
1279 1281 GSkyDir dir1;
1280 1282 dir1.radec_deg(83.6331, +22.0145);
1281   - GModelSpatialRadialRing model2(dir1, 3.0,1.0);
1282   - test_value(model2.ra(), 83.6331);
1283   - test_value(model2.dec(), 22.0145);
1284   - test_value(model2.radius(), 3.0);
1285   - test_value(model2.width(), 1.0);
  1283 + GModelSpatialRadialRing model2(dir1, 3.0, 1.0);
  1284 + test_value(model2.dir().ra_deg(), 83.6331, "Test Right Ascension for celestial sky direction constructor");
  1285 + test_value(model2.dir().dec_deg(), 22.0145, "Test Declination for celestial sky direction constructor");
  1286 + test_value(model2.dir().l_deg(), 184.557460, "Test Galactic longitude for celestial sky direction constructor");
  1287 + test_value(model2.dir().b_deg(), -5.784346, "Test Galactic latitude for celestial sky direction constructor");
  1288 + test_value(model2.radius(), 3.0, "Test radius for celestial sky direction constructor");
  1289 + test_value(model2.width(), 1.0, "Test width for celestial sky direction constructor");
  1290 +
  1291 + // Test Galactic sky direction constructor
  1292 + GModelSpatialRadialRing model2g(dir1, 3.0, 1.0, "GAL");
  1293 + test_value(model2g.dir().ra_deg(), 83.6331, "Test Right Ascension for Galactic sky direction constructor");
  1294 + test_value(model2g.dir().dec_deg(), +22.0145, "Test Declination for Galactic sky direction constructor");
  1295 + test_value(model2g.dir().l_deg(), 184.557460, "Test Galactic longitude for Galactic sky direction constructor");
  1296 + test_value(model2g.dir().b_deg(), -5.784346, "Test Galactic latitude for Galactic sky direction constructor");
  1297 + test_value(model2g.radius(), 3.0, "Test radius for Galactic sky direction constructor");
  1298 + test_value(model2g.width(), 1.0, "Test width for Galactic sky direction constructor");
1286 1299  
1287 1300 // Test XML constructor
1288 1301 GXml xml(m_xml_model_radial_ring);
... ... @@ -1290,19 +1303,11 @@ void TestGModel::test_radial_ring(void)
1290 1303 GModelSpatialRadialRing model3(*element);
1291 1304 test_value(model3.size(), 4);
1292 1305 test_value(model3.type(), "RadialRing");
1293   - test_value(model3.ra(), 83.6331);
1294   - test_value(model3.dec(), 22.0145);
  1306 + test_value(model3.dir().ra_deg(), 83.6331, "Test Right Ascension for celestial XML constructor");
  1307 + test_value(model3.dir().dec_deg(), 22.0145, "Test Declination for celestial XML constructor");
1295 1308 test_value(model3.radius(), 0.45);
1296 1309 test_value(model3.width(), 0.15);
1297 1310  
1298   - // Test ra method
1299   - model3.ra(100.0);
1300   - test_value(model3.ra(), 100.0);
1301   -
1302   - // Test dec method
1303   - model3.dec(10.0);
1304   - test_value(model3.dec(), 10.0);
1305   -
1306 1311 // Test dir method
1307 1312 GSkyDir dir2;
1308 1313 dir2.radec_deg(83.6331, +22.0145);
... ... @@ -1361,9 +1366,19 @@ void TestGModel::test_radial_gauss(void)
1361 1366 GSkyDir dir1;
1362 1367 dir1.radec_deg(83.6331, +22.0145);
1363 1368 GModelSpatialRadialGauss model2(dir1, 3.0);
1364   - test_value(model2.ra(), 83.6331);
1365   - test_value(model2.dec(), 22.0145);
1366   - test_value(model2.sigma(), 3.0);
  1369 + test_value(model2.dir().ra_deg(), 83.6331, "Test Right Ascension for celestial sky direction constructor");
  1370 + test_value(model2.dir().dec_deg(), 22.0145, "Test Declination for celestial sky direction constructor");
  1371 + test_value(model2.dir().l_deg(), 184.557460, "Test Galactic longitude for celestial sky direction constructor");
  1372 + test_value(model2.dir().b_deg(), -5.784346, "Test Galactic latitude for celestial sky direction constructor");
  1373 + test_value(model2.sigma(), 3.0, "Test sigma for celestial sky direction constructor");
  1374 +
  1375 + // Test Galactic sky direction constructor
  1376 + GModelSpatialRadialGauss model2g(dir1, 3.0, "GAL");
  1377 + test_value(model2g.dir().ra_deg(), 83.6331, "Test Right Ascension for Galactic sky direction constructor");
  1378 + test_value(model2g.dir().dec_deg(), +22.0145, "Test Declination for Galactic sky direction constructor");
  1379 + test_value(model2g.dir().l_deg(), 184.557460, "Test Galactic longitude for Galactic sky direction constructor");
  1380 + test_value(model2g.dir().b_deg(), -5.784346, "Test Galactic latitude for Galactic sky direction constructor");
  1381 + test_value(model2g.sigma(), 3.0, "Test sigma for Galactic sky direction constructor");
1367 1382  
1368 1383 // Test XML constructor
1369 1384 GXml xml(m_xml_model_radial_gauss);
... ... @@ -1371,18 +1386,10 @@ void TestGModel::test_radial_gauss(void)
1371 1386 GModelSpatialRadialGauss model3(*element);
1372 1387 test_value(model3.size(), 3);
1373 1388 test_value(model3.type(), "RadialGaussian");
1374   - test_value(model3.ra(), 83.6331);
1375   - test_value(model3.dec(), 22.0145);
  1389 + test_value(model3.dir().ra_deg(), 83.6331, "Test Right Ascension for celestial XML constructor");
  1390 + test_value(model3.dir().dec_deg(), 22.0145, "Test Declination for celestial XML constructor");
1376 1391 test_value(model3.sigma(), 0.20);
1377 1392  
1378   - // Test ra method
1379   - model3.ra(100.0);
1380   - test_value(model3.ra(), 100.0);
1381   -
1382   - // Test dec method
1383   - model3.dec(10.0);
1384   - test_value(model3.dec(), 10.0);
1385   -
1386 1393 // Test dir method
1387 1394 GSkyDir dir2;
1388 1395 dir2.radec_deg(83.6331, +22.0145);
... ... @@ -1429,10 +1436,21 @@ void TestGModel::test_radial_general_gauss(void)
1429 1436 GSkyDir dir1;
1430 1437 dir1.radec_deg(83.6331, +22.0145);
1431 1438 GModelSpatialRadialGeneralGauss model2(dir1, 3.0, 0.5);
1432   - test_value(model2.ra(), 83.6331);
1433   - test_value(model2.dec(), 22.0145);
1434   - test_value(model2.radius(), 3.0);
1435   - test_value(model2.ridx(), 0.5);
  1439 + test_value(model2.dir().ra_deg(), 83.6331, "Test Right Ascension for celestial sky direction constructor");
  1440 + test_value(model2.dir().dec_deg(), 22.0145, "Test Declination for celestial sky direction constructor");
  1441 + test_value(model2.dir().l_deg(), 184.557460, "Test Galactic longitude for celestial sky direction constructor");
  1442 + test_value(model2.dir().b_deg(), -5.784346, "Test Galactic latitude for celestial sky direction constructor");
  1443 + test_value(model2.radius(), 3.0, "Test radius for celestial sky direction constructor");
  1444 + test_value(model2.ridx(), 0.5, "Test radial index for celestial sky direction constructor");
  1445 +
  1446 + // Test Galactic sky direction constructor
  1447 + GModelSpatialRadialGeneralGauss model2g(dir1, 3.0, 0.5, "GAL");
  1448 + test_value(model2g.dir().ra_deg(), 83.6331, "Test Right Ascension for Galactic sky direction constructor");
  1449 + test_value(model2g.dir().dec_deg(), +22.0145, "Test Declination for Galactic sky direction constructor");
  1450 + test_value(model2g.dir().l_deg(), 184.557460, "Test Galactic longitude for Galactic sky direction constructor");
  1451 + test_value(model2g.dir().b_deg(), -5.784346, "Test Galactic latitude for Galactic sky direction constructor");
  1452 + test_value(model2g.radius(), 3.0, "Test radius for Galactic sky direction constructor");
  1453 + test_value(model2g.ridx(), 0.5, "Test radial index for Galactic sky direction constructor");
1436 1454  
1437 1455 // Test XML constructor
1438 1456 GXml xml(m_xml_model_radial_general_gauss);
... ... @@ -1440,19 +1458,11 @@ void TestGModel::test_radial_general_gauss(void)
1440 1458 GModelSpatialRadialGeneralGauss model3(*element);
1441 1459 test_value(model3.size(), 4);
1442 1460 test_value(model3.type(), "RadialGeneralGaussian");
1443   - test_value(model3.ra(), 83.6331);
1444   - test_value(model3.dec(), 22.0145);
  1461 + test_value(model3.dir().ra_deg(), 83.6331, "Test Right Ascension for celestial XML constructor");
  1462 + test_value(model3.dir().dec_deg(), 22.0145, "Test Declination for celestial XML constructor");
1445 1463 test_value(model3.radius(), 0.20);
1446 1464 test_value(model3.ridx(), 0.3);
1447 1465  
1448   - // Test ra method
1449   - model3.ra(100.0);
1450   - test_value(model3.ra(), 100.0);
1451   -
1452   - // Test dec method
1453   - model3.dec(10.0);
1454   - test_value(model3.dec(), 10.0);
1455   -
1456 1466 // Test dir method
1457 1467 GSkyDir dir2;
1458 1468 dir2.radec_deg(83.6331, +22.0145);
... ... @@ -1504,10 +1514,21 @@ void TestGModel::test_radial_shell(void)
1504 1514 GSkyDir dir1;
1505 1515 dir1.radec_deg(83.6331, +22.0145);
1506 1516 GModelSpatialRadialShell model2(dir1, 3.0, 1.0);
1507   - test_value(model2.ra(), 83.6331);
1508   - test_value(model2.dec(), 22.0145);
1509   - test_value(model2.radius(), 3.0);
1510   - test_value(model2.width(), 1.0);
  1517 + test_value(model2.dir().ra_deg(), 83.6331, "Test Right Ascension for celestial sky direction constructor");
  1518 + test_value(model2.dir().dec_deg(), 22.0145, "Test Declination for celestial sky direction constructor");
  1519 + test_value(model2.dir().l_deg(), 184.557460, "Test Galactic longitude for celestial sky direction constructor");
  1520 + test_value(model2.dir().b_deg(), -5.784346, "Test Galactic latitude for celestial sky direction constructor");
  1521 + test_value(model2.radius(), 3.0, "Test radius for celestial sky direction constructor");
  1522 + test_value(model2.width(), 1.0, "Test width for celestial sky direction constructor");
  1523 +
  1524 + // Test Galactic sky direction constructor
  1525 + GModelSpatialRadialShell model2g(dir1, 3.0, 1.0, "GAL");
  1526 + test_value(model2g.dir().ra_deg(), 83.6331, "Test Right Ascension for Galactic sky direction constructor");
  1527 + test_value(model2g.dir().dec_deg(), +22.0145, "Test Declination for Galactic sky direction constructor");
  1528 + test_value(model2g.dir().l_deg(), 184.557460, "Test Galactic longitude for Galactic sky direction constructor");
  1529 + test_value(model2g.dir().b_deg(), -5.784346, "Test Galactic latitude for Galactic sky direction constructor");
  1530 + test_value(model2g.radius(), 3.0, "Test radius for Galactic sky direction constructor");
  1531 + test_value(model2g.width(), 1.0, "Test width for Galactic sky direction constructor");
1511 1532  
1512 1533 // Test XML constructor
1513 1534 GXml xml(m_xml_model_radial_shell);
... ... @@ -1515,19 +1536,11 @@ void TestGModel::test_radial_shell(void)
1515 1536 GModelSpatialRadialShell model3(*element);
1516 1537 test_value(model3.size(), 4);
1517 1538 test_value(model3.type(), "RadialShell");
1518   - test_value(model3.ra(), 83.6331);
1519   - test_value(model3.dec(), 22.0145);
  1539 + test_value(model3.dir().ra_deg(), 83.6331, "Test Right Ascension for celestial XML constructor");
  1540 + test_value(model3.dir().dec_deg(), 22.0145, "Test Declination for celestial XML constructor");
1520 1541 test_value(model3.radius(), 0.30);
1521 1542 test_value(model3.width(), 0.10);
1522 1543  
1523   - // Test ra method
1524   - model3.ra(100.0);
1525   - test_value(model3.ra(), 100.0);
1526   -
1527   - // Test dec method
1528   - model3.dec(10.0);
1529   - test_value(model3.dec(), 10.0);
1530   -
1531 1544 // Test dir method
1532 1545 GSkyDir dir2;
1533 1546 dir2.radec_deg(83.6331, +22.0145);
... ... @@ -1581,8 +1594,8 @@ void TestGModel::test_profile_dmburkert(void)
1581 1594 GModelSpatialRadialProfileDMBurkert model3(*element);
1582 1595 test_value(model3.size(), 8);
1583 1596 test_value(model3.type(), "DMBurkertProfile");
1584   - test_value(model3.ra(), 83.6331);
1585   - test_value(model3.dec(), 22.0145);
  1597 + test_value(model3.dir().ra_deg(), 83.6331, "Test Right Ascension for celestial XML constructor");
  1598 + test_value(model3.dir().dec_deg(), 22.0145, "Test Declination for celestial XML constructor");
1586 1599 test_value(model3.scale_radius(), 21.5);
1587 1600 test_value(model3.scale_density(), 0.2);
1588 1601 test_value(model3.halo_distance(), 7.94);
... ... @@ -1607,8 +1620,8 @@ void TestGModel::test_profile_dmeinasto(void)
1607 1620 GModelSpatialRadialProfileDMEinasto model3(*element);
1608 1621 test_value(model3.size(), 9);
1609 1622 test_value(model3.type(), "DMEinastoProfile");
1610   - test_value(model3.ra(), 83.6331);
1611   - test_value(model3.dec(), 22.0145);
  1623 + test_value(model3.dir().ra_deg(), 83.6331, "Test Right Ascension for celestial XML constructor");
  1624 + test_value(model3.dir().dec_deg(), 22.0145, "Test Declination for celestial XML constructor");
1612 1625 test_value(model3.scale_radius(), 21.5);
1613 1626 test_value(model3.scale_density(), 0.2);
1614 1627 test_value(model3.halo_distance(), 7.94);
... ... @@ -1634,8 +1647,8 @@ void TestGModel::test_profile_dmzhao(void)
1634 1647 GModelSpatialRadialProfileDMZhao model3(*element);
1635 1648 test_value(model3.size(), 11);
1636 1649 test_value(model3.type(), "DMZhaoProfile");
1637   - test_value(model3.ra(), 83.6331);
1638   - test_value(model3.dec(), 22.0145);
  1650 + test_value(model3.dir().ra_deg(), 83.6331, "Test Right Ascension for celestial XML constructor");
  1651 + test_value(model3.dir().dec_deg(), 22.0145, "Test Declination for celestial XML constructor");
1639 1652 test_value(model3.scale_radius(), 21.5);
1640 1653 test_value(model3.scale_density(), 0.2);
1641 1654 test_value(model3.halo_distance(), 7.94);
... ... @@ -1661,8 +1674,10 @@ void TestGModel::test_elliptical_disk(void)
1661 1674 GSkyDir dir1;
1662 1675 dir1.radec_deg(83.6331, +22.0145);
1663 1676 GModelSpatialEllipticalDisk model2(dir1, 3.0, 2.0, 45.0);
1664   - test_value(model2.ra(), 83.6331);
1665   - test_value(model2.dec(), 22.0145);
  1677 + test_value(model2.dir().ra_deg(), 83.6331, "Test Right Ascension for celestial sky direction constructor");
  1678 + test_value(model2.dir().dec_deg(), 22.0145, "Test Declination for celestial sky direction constructor");
  1679 + test_value(model2.dir().l_deg(), 184.557460, "Test Galactic longitude for celestial sky direction constructor");
  1680 + test_value(model2.dir().b_deg(), -5.784346, "Test Galactic latitude for celestial sky direction constructor");
1666 1681 test_value(model2.posangle(), 45.0);
1667 1682 test_value(model2.semimajor(), 3.0);
1668 1683 test_value(model2.semiminor(), 2.0);
... ... @@ -1673,8 +1688,8 @@ void TestGModel::test_elliptical_disk(void)
1673 1688 GModelSpatialEllipticalDisk model3(*element);
1674 1689 test_value(model3.size(), 5);
1675 1690 test_value(model3.type(), "EllipticalDisk");
1676   - test_value(model3.ra(), 83.6331);
1677   - test_value(model3.dec(), 22.0145);
  1691 + test_value(model3.dir().ra_deg(), 83.6331, "Test Right Ascension for celestial XML constructor");
  1692 + test_value(model3.dir().dec_deg(), 22.0145, "Test Declination for celestial XML constructor");
1678 1693 test_value(model3.posangle(), 45.0);
1679 1694 test_value(model3.semimajor(), 2.0);
1680 1695 test_value(model3.semiminor(), 0.5);
... ... @@ -1741,8 +1756,10 @@ void TestGModel::test_elliptical_gauss(void)
1741 1756 GSkyDir dir1;
1742 1757 dir1.radec_deg(83.6331, +22.0145);
1743 1758 GModelSpatialEllipticalGauss model2(dir1, 3.0, 2.0, 45.0);
1744   - test_value(model2.ra(), 83.6331);
1745   - test_value(model2.dec(), 22.0145);
  1759 + test_value(model2.dir().ra_deg(), 83.6331, "Test Right Ascension for celestial sky direction constructor");
  1760 + test_value(model2.dir().dec_deg(), 22.0145, "Test Declination for celestial sky direction constructor");
  1761 + test_value(model2.dir().l_deg(), 184.557460, "Test Galactic longitude for celestial sky direction constructor");
  1762 + test_value(model2.dir().b_deg(), -5.784346, "Test Galactic latitude for celestial sky direction constructor");
1746 1763 test_value(model2.posangle(), 45.0);
1747 1764 test_value(model2.semimajor(), 3.0);
1748 1765 test_value(model2.semiminor(), 2.0);
... ... @@ -1753,8 +1770,8 @@ void TestGModel::test_elliptical_gauss(void)
1753 1770 GModelSpatialEllipticalGauss model3(*element);
1754 1771 test_value(model3.size(), 5);
1755 1772 test_value(model3.type(), "EllipticalGaussian");
1756   - test_value(model3.ra(), 83.6331);
1757   - test_value(model3.dec(), 22.0145);
  1773 + test_value(model3.dir().ra_deg(), 83.6331, "Test Right Ascension for celestial XML constructor");
  1774 + test_value(model3.dir().dec_deg(), 22.0145, "Test Declination for celestial XML constructor");
1758 1775 test_value(model3.posangle(), 45.0);
1759 1776 test_value(model3.semimajor(), 0.3);
1760 1777 test_value(model3.semiminor(), 0.1);
... ... @@ -1821,8 +1838,10 @@ void TestGModel::test_elliptical_general_gauss(void)
1821 1838 GSkyDir dir1;
1822 1839 dir1.radec_deg(83.6331, +22.0145);
1823 1840 GModelSpatialEllipticalGeneralGauss model2(dir1, 3.0, 2.0, 45.0,0.5);
1824   - test_value(model2.ra(), 83.6331);
1825   - test_value(model2.dec(), 22.0145);
  1841 + test_value(model2.dir().ra_deg(), 83.6331, "Test Right Ascension for celestial sky direction constructor");
  1842 + test_value(model2.dir().dec_deg(), 22.0145, "Test Declination for celestial sky direction constructor");
  1843 + test_value(model2.dir().l_deg(), 184.557460, "Test Galactic longitude for celestial sky direction constructor");
  1844 + test_value(model2.dir().b_deg(), -5.784346, "Test Galactic latitude for celestial sky direction constructor");
1826 1845 test_value(model2.posangle(), 45.0);
1827 1846 test_value(model2.semimajor(), 3.0);
1828 1847 test_value(model2.semiminor(), 2.0);
... ... @@ -1835,8 +1854,8 @@ void TestGModel::test_elliptical_general_gauss(void)
1835 1854 GModelSpatialEllipticalGeneralGauss model3(*element);
1836 1855 test_value(model3.size(), 6);
1837 1856 test_value(model3.type(), "EllipticalGeneralGaussian");
1838   - test_value(model3.ra(), 83.6331);
1839   - test_value(model3.dec(), 22.0145);
  1857 + test_value(model3.dir().ra_deg(), 83.6331, "Test Right Ascension for celestial XML constructor");
  1858 + test_value(model3.dir().dec_deg(), 22.0145, "Test Declination for celestial XML constructor");
1840 1859 test_value(model3.posangle(), 45.0);
1841 1860 test_value(model3.semimajor(), 0.3);
1842 1861 test_value(model3.semiminor(), 0.1);
... ... @@ -3913,8 +3932,8 @@ void TestGModel::test_legacy_model_radial_disk(void)
3913 3932  
3914 3933 // Test model values
3915 3934 test_value(spatial->type(), "DiskFunction", "Check model type");
3916   - test_value(spatial->ra(), 83.6331, 1.0e-7, "Check Right Ascension");
3917   - test_value(spatial->dec(), 22.0145, 1.0e-7, "Check Declination");
  3935 + test_value(spatial->dir().ra_deg(), 83.6331, 1.0e-7, "Check Right Ascension");
  3936 + test_value(spatial->dir().dec_deg(), 22.0145, 1.0e-7, "Check Declination");
3918 3937 test_value(spatial->radius(), 0.45, 1.0e-7, "Check radius");
3919 3938  
3920 3939 // Save file to disk and reload file (tests the proper saving and loading)
... ... @@ -3928,8 +3947,8 @@ void TestGModel::test_legacy_model_radial_disk(void)
3928 3947  
3929 3948 // Re-test model values
3930 3949 test_value(spatial->type(), "DiskFunction", "Check model type");
3931   - test_value(spatial->ra(), 83.6331, 1.0e-7, "Check Right Ascension");
3932   - test_value(spatial->dec(), 22.0145, 1.0e-7, "Check Declination");
  3950 + test_value(spatial->dir().ra_deg(), 83.6331, 1.0e-7, "Check Right Ascension");
  3951 + test_value(spatial->dir().dec_deg(), 22.0145, 1.0e-7, "Check Declination");
3933 3952 test_value(spatial->radius(), 0.45, 1.0e-7, "Check radius");
3934 3953  
3935 3954 // Exit test
... ... @@ -3952,8 +3971,8 @@ void TestGModel::test_legacy_model_radial_gauss(void)
3952 3971  
3953 3972 // Test model values
3954 3973 test_value(spatial->type(), "GaussFunction", "Check model type");
3955   - test_value(spatial->ra(), 83.6331, 1.0e-7, "Check Right Ascension");
3956   - test_value(spatial->dec(), 22.0145, 1.0e-7, "Check Declination");
  3974 + test_value(spatial->dir().ra_deg(), 83.6331, 1.0e-7, "Check Right Ascension");
  3975 + test_value(spatial->dir().dec_deg(), 22.0145, 1.0e-7, "Check Declination");
3957 3976 test_value(spatial->sigma(), 0.20, 1.0e-7, "Check Sigma");
3958 3977  
3959 3978 // Save file to disk and reload file (tests the proper saving and loading)
... ... @@ -3967,8 +3986,8 @@ void TestGModel::test_legacy_model_radial_gauss(void)
3967 3986  
3968 3987 // Re-test model values
3969 3988 test_value(spatial->type(), "GaussFunction", "Check model type");
3970   - test_value(spatial->ra(), 83.6331, 1.0e-7, "Check Right Ascension");
3971   - test_value(spatial->dec(), 22.0145, 1.0e-7, "Check Declination");
  3989 + test_value(spatial->dir().ra_deg(), 83.6331, 1.0e-7, "Check Right Ascension");
  3990 + test_value(spatial->dir().dec_deg(), 22.0145, 1.0e-7, "Check Declination");
3972 3991 test_value(spatial->sigma(), 0.20, 1.0e-7, "Check Sigma");
3973 3992  
3974 3993 // Exit test
... ... @@ -3991,8 +4010,8 @@ void TestGModel::test_legacy_model_radial_shell(void)
3991 4010  
3992 4011 // Test model values
3993 4012 test_value(spatial->type(), "ShellFunction", "Check model type");
3994   - test_value(spatial->ra(), 83.6331, 1.0e-7, "Check Right Ascension");
3995   - test_value(spatial->dec(), 22.0145, 1.0e-7, "Check Declination");
  4013 + test_value(spatial->dir().ra_deg(), 83.6331, 1.0e-7, "Check Right Ascension");
  4014 + test_value(spatial->dir().dec_deg(), 22.0145, 1.0e-7, "Check Declination");
3996 4015 test_value(spatial->radius(), 0.30, 1.0e-7, "Check radius");
3997 4016 test_value(spatial->width(), 0.10, 1.0e-7, "Check width");
3998 4017  
... ... @@ -4007,8 +4026,8 @@ void TestGModel::test_legacy_model_radial_shell(void)
4007 4026  
4008 4027 // Re-test model values
4009 4028 test_value(spatial->type(), "ShellFunction", "Check model type");
4010   - test_value(spatial->ra(), 83.6331, 1.0e-7, "Check Right Ascension");
4011   - test_value(spatial->dec(), 22.0145, 1.0e-7, "Check Declination");
  4029 + test_value(spatial->dir().ra_deg(), 83.6331, 1.0e-7, "Check Right Ascension");
  4030 + test_value(spatial->dir().dec_deg(), 22.0145, 1.0e-7, "Check Declination");
4012 4031 test_value(spatial->radius(), 0.30, 1.0e-7, "Check radius");
4013 4032 test_value(spatial->width(), 0.10, 1.0e-7, "Check width");
4014 4033  
... ... @@ -4032,8 +4051,8 @@ void TestGModel::test_legacy_model_elliptical_gauss(void)
4032 4051  
4033 4052 // Test model values
4034 4053 test_value(spatial->type(), "EllipticalGauss", "Check model type");
4035   - test_value(spatial->ra(), 83.6331, 1.0e-7, "Check Right Ascension");
4036   - test_value(spatial->dec(), 22.0145, 1.0e-7, "Check Declination");
  4054 + test_value(spatial->dir().ra_deg(), 83.6331, 1.0e-7, "Check Right Ascension");
  4055 + test_value(spatial->dir().dec_deg(), 22.0145, 1.0e-7, "Check Declination");
4037 4056 test_value(spatial->posangle(), 45.0, 1.0e-7, "Check position angle");
4038 4057 test_value(spatial->semimajor(), 0.3, 1.0e-7, "Check semi major axis");
4039 4058 test_value(spatial->semiminor(), 0.1, 1.0e-7, "Check semi minor axis");
... ... @@ -4049,8 +4068,8 @@ void TestGModel::test_legacy_model_elliptical_gauss(void)
4049 4068  
4050 4069 // Re-test model values
4051 4070 test_value(spatial->type(), "EllipticalGauss", "Check model type");
4052   - test_value(spatial->ra(), 83.6331, 1.0e-7, "Check Right Ascension");
4053   - test_value(spatial->dec(), 22.0145, 1.0e-7, "Check Declination");
  4071 + test_value(spatial->dir().ra_deg(), 83.6331, 1.0e-7, "Check Right Ascension");
  4072 + test_value(spatial->dir().dec_deg(), 22.0145, 1.0e-7, "Check Declination");
4054 4073 test_value(spatial->posangle(), 45.0, 1.0e-7, "Check position angle");
4055 4074 test_value(spatial->semimajor(), 0.3, 1.0e-7, "Check semi major axis");
4056 4075 test_value(spatial->semiminor(), 0.1, 1.0e-7, "Check semi minor axis");
... ...