Johan Bregeon / gammalib

Commit 244ac5bb95eedb5205eb65ba778d56efbea4c06e

Authored by Michael Mayer
1 parent d0351202

Implemented GModelSpatialEllipticalDisk and fixed bug in GModelSpatialElliptical

Inline Side-by-side
Showing 9 changed files with 838 additions and 8 deletions
configure.ac
  View file @244ac5b
... ... @@ -35,7 +35,7 @@ AC_CONFIG_MACRO_DIR([m4])
35 35 #AM_INIT_AUTOMAKE([gnu no-dependencies])
36 36 AM_INIT_AUTOMAKE([gnu no-dependencies -Wno-portability])
37 37 #AM_INIT_AUTOMAKE([gnu -Wno-portability])
38   -AM_CONFIG_HEADER([config.h])
  38 +AC_CONFIG_HEADER([config.h])
39 39  
40 40  
41 41 #############################################################################
... ...
include/GModelSpatialEllipticalDisk.hpp 0 → 100644
  View file @244ac5b
  1 +/***************************************************************************
  2 + * GModelSpatialEllipticalDisk.hpp - Elliptical disk source model class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2011-2013 by Christoph Deil *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GModelSpatialEllipticalDisk.hpp
  23 + * @brief Elliptical disk model class interface definition
  24 + * @author Christoph Deil
  25 + */
  26 +
  27 +#ifndef GMODELSPATIALELLIPTICALDISK_HPP
  28 +#define GMODELSPATIALELLIPTICALDISK_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include "GModelSpatialElliptical.hpp"
  33 +#include "GModelPar.hpp"
  34 +#include "GSkyDir.hpp"
  35 +#include "GXmlElement.hpp"
  36 +
  37 +
  38 +/**************************************************************************
  39 + * @class GModelSpatialEllipticalDisk
  40 + *
  41 + * @brief Disk source model class
  42 + *
  43 + * This class implements the spatial component of the factorised source
  44 + * model for a disk source, i.e. constant surface brightness within some
  45 + * radius and no emission outside.
  46 + ***************************************************************************/
  47 +class GModelSpatialEllipticalDisk : public GModelSpatialElliptical {
  48 +
  49 +public:
  50 + // Constructors and destructors
  51 + GModelSpatialEllipticalDisk(void);
  52 + explicit GModelSpatialEllipticalDisk(const GSkyDir& dir, const double& minor,const double& major, const double& posangle);
  53 + explicit GModelSpatialEllipticalDisk(const GXmlElement& xml);
  54 + GModelSpatialEllipticalDisk(const GModelSpatialEllipticalDisk& model);
  55 + virtual ~GModelSpatialEllipticalDisk(void);
  56 +
  57 + // Operators
  58 + virtual GModelSpatialEllipticalDisk& operator=(const GModelSpatialEllipticalDisk& model);
  59 +
  60 + // Implemented pure virtual methods
  61 + virtual void clear(void);
  62 + virtual GModelSpatialEllipticalDisk* clone(void) const;
  63 + virtual std::string type(void) const { return "EllipticalDisk"; }
  64 + virtual double eval(const double& theta, const double &posangle) const;
  65 + virtual double eval_gradients(const double& theta, const double &posangle) const;
  66 + virtual GSkyDir mc(GRan& ran) const;
  67 + virtual double theta_max(void) const;
  68 + virtual void read(const GXmlElement& xml);
  69 + virtual void write(GXmlElement& xml) const;
  70 + virtual std::string print(void) const;
  71 +
  72 +// virtual double eval(const GSkyDir& srcDir) const;
  73 +// virtual double eval_gradients(const GSkyDir& srcDir) const;
  74 +
  75 +
  76 + // Other methods
  77 + double minor(void) const { return m_minor.real_value(); }
  78 + double major(void) const { return m_major.real_value(); }
  79 +
  80 + void minor(const double& minor) { m_minor.real_value(minor); }
  81 + void major(const double& major) { m_major.real_value(major); }
  82 +
  83 +protected:
  84 + // Protected methods
  85 + void init_members(void);
  86 + void copy_members(const GModelSpatialEllipticalDisk& model);
  87 + void free_members(void);
  88 + void update(void) const;
  89 +
  90 + // Protected members
  91 + GModelPar m_minor; //!< minor axis of ellipse (degrees)
  92 + GModelPar m_major; //!< major axis of ellipse (degrees)
  93 +
  94 + // Cached members used for pre-computations
  95 + mutable double m_last_minor; //!< Last minor axis
  96 + mutable double m_last_major; //!< Last major axis
  97 + mutable double m_minor_rad; //!< Radius in radians
  98 + mutable double m_major_rad; //!< Radius in radians
  99 + mutable double m_norm; //!< Normalization
  100 +};
  101 +
  102 +#endif /* GMODELSPATIALELLIPTICALDISK_HPP */
... ...
include/GammaLib.hpp
  View file @244ac5b
... ... @@ -171,6 +171,7 @@
171 171 #include "GModelSpatialRadialGauss.hpp"
172 172 #include "GModelSpatialRadialShell.hpp"
173 173 #include "GModelSpatialElliptical.hpp"
  174 +#include "GModelSpatialEllipticalDisk.hpp"
174 175 #include "GModelSpatialDiffuse.hpp"
175 176 #include "GModelSpatialDiffuseConst.hpp"
176 177 #include "GModelSpatialDiffuseCube.hpp"
... ...
include/Makefile.am
  View file @244ac5b
... ... @@ -116,6 +116,7 @@ pkginclude_HEADERS = GammaLib.hpp \
116 116 GModelSpatialRadialGauss.hpp \
117 117 GModelSpatialRadialShell.hpp \
118 118 GModelSpatialElliptical.hpp \
  119 + GModelSpatialEllipticalDisk.hpp \
119 120 GModelSpatialDiffuse.hpp \
120 121 GModelSpatialDiffuseConst.hpp \
121 122 GModelSpatialDiffuseCube.hpp \
... ...
pyext/GModelSpatialEllipticalDisk.i 0 → 100644
  View file @244ac5b
  1 +/***************************************************************************
  2 + * GModelSpatialEllipticalDisk.i - Abstract elliptical disc model class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Michael Mayer *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GModelSpatialEllipticalDisk.i
  23 + * @brief Elliptical disk model class Python interface
  24 + * @author Michael Mayer
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GModelSpatialEllipticalDisk.hpp"
  29 +#include "GTools.hpp"
  30 +%}
  31 +
  32 +
  33 +/***********************************************************************//**
  34 + * @class GModelSpatialEllipticalDisk
  35 + *
  36 + * @brief elliptical disk model class
  37 + ***************************************************************************/
  38 +class GModelSpatialEllipticalDisk : public GModelSpatialElliptical {
  39 +public:
  40 + // Constructors and destructors
  41 + GModelSpatialEllipticalDisk(void);
  42 + GModelSpatialEllipticalDisk(const GModelSpatialEllipticalDisk& model);
  43 + explicit GModelSpatialEllipticalDisk(const GSkyDir& dir, const double& minor,const double& major, const double& posangle);
  44 + explicit GModelSpatialEllipticalDisk(const GXmlElement& xml);
  45 + virtual ~GModelSpatialEllipticalDisk(void);
  46 +
  47 + // Pure virtual methods
  48 + virtual void clear(void);
  49 + virtual GModelSpatialEllipticalDisk* clone(void) const;
  50 + virtual std::string type(void) const;
  51 + virtual double eval(const double& theta,const double& posangle) const;
  52 + virtual double eval_gradients(const double& theta,const double& posangle)const;
  53 + virtual GSkyDir mc(GRan& ran) const;
  54 + virtual double theta_max(void) const;
  55 +
  56 + // Implemented virtual methods
  57 + virtual void read(const GXmlElement& xml);
  58 + virtual void write(GXmlElement& xml) const;
  59 +
  60 + // Other methods
  61 + double ra(void) const;
  62 + double dec(void) const;
  63 + double posangle(void) const;
  64 + GSkyDir dir(void) const;
  65 + void dir(const GSkyDir& dir);
  66 + void posangle(const double& posangle);
  67 +};
  68 +
  69 +
  70 +/***********************************************************************//**
  71 + * @brief GModelSpatialEllipticalDisk class extension
  72 + ***************************************************************************/
  73 +%extend GModelSpatialEllipticalDisk {
  74 + GModelSpatialEllipticalDisk copy() {
  75 + return (*self);
  76 + }
  77 + double eval(const GSkyDir& srcDir) const {
  78 + return self->GModelSpatialElliptical::eval(srcDir);
  79 + }
  80 + double eval_gradients(const GSkyDir& srcDir) const {
  81 + return self->GModelSpatialElliptical::eval_gradients(srcDir);
  82 + }
  83 +};
... ...
pyext/gammalib/model.i
  View file @244ac5b
... ... @@ -58,6 +58,7 @@
58 58 %include "GModelSpatialRadialGauss.i"
59 59 %include "GModelSpatialRadialShell.i"
60 60 %include "GModelSpatialElliptical.i"
  61 +%include "GModelSpatialEllipticalDisk.i"
61 62 %include "GModelSpatialDiffuse.i"
62 63 %include "GModelSpatialDiffuseConst.i"
63 64 %include "GModelSpatialDiffuseCube.i"
... ...
src/model/GModelSpatialElliptical.cpp
  View file @244ac5b
... ... @@ -456,18 +456,19 @@ void GModelSpatialElliptical::init_members(void)
456 456 m_dec.hasgrad(false);
457 457  
458 458 // Initialise Position Angle
459   - m_dec.clear();
460   - m_dec.name("PA");
461   - m_dec.unit("deg");
462   - m_dec.fix();
463   - m_dec.scale(1.0);
464   - m_dec.gradient(0.0);
465   - m_dec.hasgrad(false);
  459 + m_posangle.clear();
  460 + m_posangle.name("PA");
  461 + m_posangle.unit("deg");
  462 + m_posangle.fix();
  463 + m_posangle.scale(1.0);
  464 + m_posangle.gradient(0.0);
  465 + m_posangle.hasgrad(false);
466 466  
467 467 // Set parameter pointer(s)
468 468 m_pars.clear();
469 469 m_pars.push_back(&m_ra);
470 470 m_pars.push_back(&m_dec);
  471 + m_pars.push_back(&m_posangle);
471 472  
472 473 // Return
473 474 return;
... ...
src/model/GModelSpatialEllipticalDisk.cpp 0 → 100644
  View file @244ac5b
  1 +/***************************************************************************
  2 + * GModelSpatialEllipticalDisk.cpp - Elliptical disk source model class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2011-2013 by Michael Mayer *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GModelSpatialEllipticalDisk.cpp
  23 + * @brief Elliptical disk model class implementation
  24 + * @author Michael Mayer
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "GException.hpp"
  32 +#include "GTools.hpp"
  33 +#include "GModelSpatialEllipticalDisk.hpp"
  34 +#include "GModelSpatialRadialDisk.hpp"
  35 +#include "GModelSpatialRegistry.hpp"
  36 +
  37 +/* __ Constants __________________________________________________________ */
  38 +
  39 +/* __ Globals ____________________________________________________________ */
  40 +const GModelSpatialEllipticalDisk g_elliptical_disk_seed;
  41 +const GModelSpatialRegistry g_elliptical_disk_registry(&g_elliptical_disk_seed);
  42 +
  43 +/* __ Method name definitions ____________________________________________ */
  44 +#define G_READ "GModelSpatialEllipticalDisk::read(GXmlElement&)"
  45 +#define G_WRITE "GModelSpatialEllipticalDisk::write(GXmlElement&)"
  46 +
  47 +/* __ Macros _____________________________________________________________ */
  48 +
  49 +/* __ Coding definitions _________________________________________________ */
  50 +
  51 +/* __ Debug definitions __________________________________________________ */
  52 +
  53 +
  54 +/*==========================================================================
  55 + = =
  56 + = Constructors/destructors =
  57 + = =
  58 + ==========================================================================*/
  59 +
  60 +/***********************************************************************//**
  61 + * @brief Void constructor
  62 + ***************************************************************************/
  63 +GModelSpatialEllipticalDisk::GModelSpatialEllipticalDisk(void) : GModelSpatialElliptical()
  64 +{
  65 + // Initialise members
  66 + init_members();
  67 +
  68 + // Return
  69 + return;
  70 +}
  71 +
  72 +
  73 +/***********************************************************************//**
  74 + * @brief Disk constructor
  75 + *
  76 + * @param[in] dir Sky position of disk centre.
  77 + * @param[in] minor Semi minor axis (degrees).
  78 + * @param[in] major Semi major axis (degrees).
  79 + * @param[in] posangle Position angle of major axis (degrees).
  80 + ***************************************************************************/
  81 +GModelSpatialEllipticalDisk::GModelSpatialEllipticalDisk(const GSkyDir& dir,
  82 + const double& minor,const double& major, const double & posangle) :
  83 + GModelSpatialElliptical()
  84 +{
  85 + // Initialise members
  86 + init_members();
  87 +
  88 + // Assign parameters
  89 + this->dir(dir);
  90 + this->minor(minor);
  91 + this->major(major);
  92 + this->posangle(posangle);
  93 +
  94 + // Return
  95 + return;
  96 +}
  97 +
  98 +
  99 +/***********************************************************************//**
  100 + * @brief XML constructor
  101 + *
  102 + * @param[in] xml XML element.
  103 + *
  104 + * Creates instance of Elliptical disk model by extracting information from an
  105 + * XML element. See GModelSpatialEllipticalDisk::read() for more information about the
  106 + * expected structure of the XML element.
  107 + ***************************************************************************/
  108 +GModelSpatialEllipticalDisk::GModelSpatialEllipticalDisk(const GXmlElement& xml) :
  109 + GModelSpatialElliptical()
  110 +{
  111 + // Initialise members
  112 + init_members();
  113 +
  114 + // Read information from XML element
  115 + read(xml);
  116 +
  117 + // Return
  118 + return;
  119 +}
  120 +
  121 +
  122 +/***********************************************************************//**
  123 + * @brief Copy constructor
  124 + *
  125 + * @param[in] model Elliptical disk model.
  126 + ***************************************************************************/
  127 +GModelSpatialEllipticalDisk::GModelSpatialEllipticalDisk(const GModelSpatialEllipticalDisk& model) :
  128 + GModelSpatialElliptical(model)
  129 +{
  130 + // Initialise members
  131 + init_members();
  132 +
  133 + // Copy members
  134 + copy_members(model);
  135 +
  136 + // Return
  137 + return;
  138 +}
  139 +
  140 +
  141 +/***********************************************************************//**
  142 + * @brief Destructor
  143 + ***************************************************************************/
  144 +GModelSpatialEllipticalDisk::~GModelSpatialEllipticalDisk(void)
  145 +{
  146 + // Free members
  147 + free_members();
  148 +
  149 + // Return
  150 + return;
  151 +}
  152 +
  153 +
  154 +/*==========================================================================
  155 + = =
  156 + = Operators =
  157 + = =
  158 + ==========================================================================*/
  159 +
  160 +/***********************************************************************//**
  161 + * @brief Assignment operator
  162 + *
  163 + * @param[in] model Elliptical disk model.
  164 + * @return Elliptical disk model.
  165 + ***************************************************************************/
  166 +GModelSpatialEllipticalDisk& GModelSpatialEllipticalDisk::operator=(const GModelSpatialEllipticalDisk& model)
  167 +{
  168 + // Execute only if object is not identical
  169 + if (this != &model) {
  170 +
  171 + // Copy base class members
  172 + this->GModelSpatialElliptical::operator=(model);
  173 +
  174 + // Free members
  175 + free_members();
  176 +
  177 + // Initialise members
  178 + init_members();
  179 +
  180 + // Copy members
  181 + copy_members(model);
  182 +
  183 + } // endif: object was not identical
  184 +
  185 + // Return
  186 + return *this;
  187 +}
  188 +
  189 +
  190 +/*==========================================================================
  191 + = =
  192 + = Public methods =
  193 + = =
  194 + ==========================================================================*/
  195 +
  196 +/***********************************************************************//**
  197 + * @brief Clear instance
  198 + ***************************************************************************/
  199 +void GModelSpatialEllipticalDisk::clear(void)
  200 +{
  201 + // Free class members (base and derived classes, derived class first)
  202 + free_members();
  203 + this->GModelSpatialElliptical::free_members();
  204 + this->GModelSpatial::free_members();
  205 +
  206 + // Initialise members
  207 + this->GModelSpatial::init_members();
  208 + this->GModelSpatialElliptical::init_members();
  209 + init_members();
  210 +
  211 + // Return
  212 + return;
  213 +}
  214 +
  215 +
  216 +/***********************************************************************//**
  217 + * @brief Clone instance
  218 + ***************************************************************************/
  219 +GModelSpatialEllipticalDisk* GModelSpatialEllipticalDisk::clone(void) const
  220 +{
  221 + return new GModelSpatialEllipticalDisk(*this);
  222 +}
  223 +
  224 +
  225 +/***********************************************************************//**
  226 + * @brief Evaluate function (in units of sr^-1)
  227 + *
  228 + * @param[in] theta Angular distance from disk centre (radians).
  229 + * @param[in] posangle Position angle (clockwise from North) (radians).
  230 + *
  231 + * Evaluates the spatial part for a elliptical disk source model. The disk source
  232 + * model is a Elliptical function \f$f(\theta, \phi)\f$, where \f$\theta\f$ is the
  233 + * angular separation between elliptical disk centre and the actual location and \f$\phi\f$ the
  234 + * rotation angle (clockwise from North):@TODO edit formula
  235 + * \f[
  236 + * f(\theta,\phi) = \left \{
  237 + * \begin{array}{l l}
  238 + * \displaystyle
  239 + * {\tt m\_norm}
  240 + * & \mbox{if $\theta \le $ radius} \\
  241 + * \\
  242 + * \displaystyle
  243 + * 0 & \mbox{if $\theta > $ radius}
  244 + * \end{array}
  245 + * \right .
  246 + * \f]
  247 + *
  248 + * where
  249 + *
  250 + * \f[ {\tt m\_norm} = \frac{1}{2 \pi (1 - \cos r)} \f]
  251 + *
  252 + * is a normalization constant (see the update() method).
  253 + *
  254 + * @TODO check formula
  255 + ***************************************************************************/
  256 +double GModelSpatialEllipticalDisk::eval(const double& theta, const double& posangle) const
  257 +{
  258 + // Update precomputation cache
  259 + update();
  260 +
  261 +
  262 + double diff_angle = posangle - m_posangle.real_value()*deg2rad;
  263 + double cosinus = std::cos(diff_angle);
  264 + double sinus = std::sin(diff_angle);
  265 +
  266 + double r_ell = m_minor_rad*m_major_rad/std::sqrt(std::pow(m_minor_rad*cosinus,2)+std::pow(m_major_rad*sinus,2));
  267 + double value = (theta <= r_ell) ? m_norm : 0.0;
  268 +
  269 + // Compile option: Check for NaN/Inf
  270 + #if defined(G_NAN_CHECK)
  271 + if (isnotanumber(value) || isinfinite(value)) {
  272 + std::cout << "*** ERROR: GModelSpatialEllipticalDisk::eval";
  273 + std::cout << "(theta=" << theta << "): NaN/Inf encountered";
  274 + std::cout << "(posangle=" << posangle << "): NaN/Inf encountered";
  275 + std::cout << " (value=" << value;
  276 + std::cout << ", R_ellipse=" << r_ell;
  277 + std::cout << ", diff_angle=" << diff_angle;
  278 + std::cout << ", m_norm=" << m_norm;
  279 + std::cout << ")" << std::endl;
  280 + }
  281 + #endif
  282 +
  283 + // Return value
  284 + return value;
  285 +}
  286 +
  287 +
  288 +/***********************************************************************//**
  289 + * @brief Evaluate function and gradients (in units of sr^-1)
  290 + *
  291 + * @param[in] theta Angular distance from disk centre (radians).
  292 + * @param[in] posangle Position angle for sky position(radians).
  293 + *
  294 + * Evaluates the function value. No gradient computation is implemented as
  295 + * Elliptical models will be convolved with the instrument response and thus
  296 + * require the numerical computation of the derivatives.
  297 + *
  298 + * See the eval() method for more information.
  299 + ***************************************************************************/
  300 +double GModelSpatialEllipticalDisk::eval_gradients(const double& theta, const double& posangle) const
  301 +{
  302 + // Return value
  303 + return (eval(theta, posangle));
  304 +}
  305 +
  306 +
  307 +/***********************************************************************//**
  308 + * @brief Returns MC sky direction
  309 + *
  310 + * @param[in] ran Random number generator.
  311 + *
  312 + * Draws an arbitrary sky position from the 2D disk distribution.
  313 + * @TODO test function
  314 + ***************************************************************************/
  315 +GSkyDir GModelSpatialEllipticalDisk::mc(GRan& ran) const
  316 +{
  317 + update();
  318 +
  319 + // Initialise a radial disk model with radius of the major axis
  320 + GModelSpatialRadialDisk disk = GModelSpatialRadialDisk(dir(),major());
  321 +
  322 + // Initialise SkyDir
  323 + GSkyDir dir = GSkyDir();
  324 +
  325 + // Draw randomly from the radial disk
  326 + // and reject the value if its outside the ellipse
  327 + do {
  328 +
  329 + // Set SkyDir to random position
  330 + dir = disk.mc(ran);
  331 +
  332 + }while(GModelSpatialElliptical::eval(dir)<=0);
  333 +
  334 + // Return SkyDir
  335 + return dir;
  336 +
  337 +}
  338 +
  339 +
  340 +/***********************************************************************//**
  341 + * @brief Return maximum model radius (in radians)
  342 + ***************************************************************************/
  343 +double GModelSpatialEllipticalDisk::theta_max(void) const
  344 +{
  345 + // Return value
  346 + return (major()*deg2rad);
  347 +}
  348 +
  349 +
  350 +/***********************************************************************//**
  351 + * @brief Read model from XML element
  352 + *
  353 + * @param[in] xml XML element.
  354 + *
  355 + * @exception GException::model_invalid_parnum
  356 + * Invalid number of model parameters found in XML element.
  357 + * @exception GException::model_invalid_parnames
  358 + * Invalid model parameter names found in XML element.
  359 + *
  360 + * Read the elliptical disk information from an XML element. The XML element
  361 + * is required to have 5 parameters.
  362 + * The position is named either "RA" and "DEC" or "GLON" and "GLAT", the
  363 + * rotation angle "PA", the semi-minor axis "MinorRadius" and the semi-major
  364 + * axis "MajorRadius"
  365 + *
  366 + * @todo Implement a test of the ellipse boundary. The axes
  367 + * and axes minimum should be >0.
  368 + ***************************************************************************/
  369 +void GModelSpatialEllipticalDisk::read(const GXmlElement& xml)
  370 +{
  371 + // Determine number of parameter nodes in XML element
  372 + int npars = xml.elements("parameter");
  373 +
  374 + // Verify that XML element has exactly 3 parameters
  375 + if (xml.elements() != 5 || npars != 5) {
  376 + throw GException::model_invalid_parnum(G_READ, xml,
  377 + "Elliptical Disk model requires exactly 5 parameters.");
  378 + }
  379 +
  380 + // Read disk location
  381 + GModelSpatialElliptical::read(xml);
  382 +
  383 + // Extract model parameters
  384 + int npar[2] = {0};
  385 + for (int i = 0; i < npars; ++i) {
  386 +
  387 + // Get parameter element
  388 + GXmlElement* par = static_cast<GXmlElement*>(xml.element("parameter", i));
  389 +
  390 + // Handle Radius
  391 + if (par->attribute("name") == "MinorRadius") {
  392 +
  393 + // Read parameter
  394 + m_minor.read(*par);
  395 +
  396 + // Increment parameter counter
  397 + npar[0]++;
  398 + }
  399 +
  400 + else if (par->attribute("name") == "MajorRadius") {
  401 +
  402 + // Read parameter
  403 + m_major.read(*par);
  404 +
  405 + // Increment parameter counter
  406 + npar[1]++;
  407 + }
  408 +
  409 + } // endfor: looped over all parameters
  410 +
  411 + // Verify that all parameters were found
  412 + if (npar[0] != 1 || npar[1] != 1) {
  413 + throw GException::model_invalid_parnames(G_READ, xml,
  414 + "EllipticalDisk requires \"MinorRadius\" and \"MajorRadius\" parameter.");
  415 + }
  416 +
  417 + // Return
  418 + return;
  419 +}
  420 +
  421 +
  422 +/***********************************************************************//**
  423 + * @brief Write model into XML element
  424 + *
  425 + * @param[in] xml XML element into which model information is written.
  426 + *
  427 + * @exception GException::model_invalid_parnum
  428 + * Invalid number of model parameters found in XML element.
  429 + * @exception GException::model_invalid_parnames
  430 + * Invalid model parameter names found in XML element.
  431 + *
  432 + * Write the Disk source information into an XML element. The XML element
  433 + * will have 5 parameter leafs named "RA", "DEC", "PA", "MinorRadius" and
  434 + * "MajorRadius"
  435 + ***************************************************************************/
  436 +void GModelSpatialEllipticalDisk::write(GXmlElement& xml) const
  437 +{
  438 + // Write disk location
  439 + GModelSpatialElliptical::write(xml);
  440 +
  441 + // If XML element has 2 nodes (which should be the location nodes)
  442 + // then append 1 parameter node
  443 + if (xml.elements() == 2) {
  444 + xml.append(new GXmlElement("parameter name=\"MinorRadius\""));
  445 + xml.append(new GXmlElement("parameter name=\"MajorRadius\""));
  446 + }
  447 +
  448 + // Determine number of parameter nodes in XML element
  449 + int npars = xml.elements("parameter");
  450 +
  451 + // Verify that XML element has exactly 3 parameters
  452 + if (xml.elements() != 5 || npars != 5) {
  453 + throw GException::model_invalid_parnum(G_WRITE, xml,
  454 + "Elliptical Disk model requires exactly 2 parameters.");
  455 + }
  456 +
  457 + // Set or update model parameter attributes
  458 + int npar[1] = {0};
  459 + for (int i = 0; i < npars; ++i) {
  460 +
  461 + // Get parameter element
  462 + GXmlElement* par = static_cast<GXmlElement*>(xml.element("parameter", i));
  463 +
  464 + // Handle Radius
  465 + if (par->attribute("name") == "MinorRadius") {
  466 +
  467 + // Write parameter
  468 + m_minor.write(*par);
  469 +
  470 + // Increment parameter counter
  471 + npar[0]++;
  472 + }
  473 +
  474 + else if (par->attribute("name") == "MajorRadius") {
  475 +
  476 + // Write parameter
  477 + m_major.write(*par);
  478 +
  479 + // Increment parameter counter
  480 + npar[1]++;
  481 + }
  482 +
  483 + } // endfor: looped over all parameters
  484 +
  485 + // Check of all required parameters are present
  486 + if (npar[0] != 1 || npar[1] != 1) {
  487 + throw GException::model_invalid_parnames(G_WRITE, xml,
  488 + "Require \"MinorRadius\" and \"MajorRadius\" parameter.");
  489 + }
  490 +
  491 + // Return
  492 + return;
  493 +}
  494 +
  495 +
  496 +/***********************************************************************//**
  497 + * @brief Print information
  498 + *
  499 + * @return String containing model information.
  500 + ***************************************************************************/
  501 +std::string GModelSpatialEllipticalDisk::print(void) const
  502 +{
  503 + // Initialise result string
  504 + std::string result;
  505 +
  506 + // Append header
  507 + result.append("=== GModelSpatialEllipticalDisk ===\n");
  508 +
  509 + // Append parameters
  510 + result.append(parformat("Number of parameters")+str(size()));
  511 + for (int i = 0; i < size(); ++i) {
  512 + result.append("\n"+m_pars[i]->print());
  513 + }
  514 +
  515 + // Return result
  516 + return result;
  517 +}
  518 +
  519 +
  520 +/*==========================================================================
  521 + = =
  522 + = Private methods =
  523 + = =
  524 + ==========================================================================*/
  525 +
  526 +/***********************************************************************//**
  527 + * @brief Initialise class members
  528 + ***************************************************************************/
  529 +void GModelSpatialEllipticalDisk::init_members(void)
  530 +{
  531 + // Initialise minor axis
  532 + m_minor.clear();
  533 + m_minor.name("MinorRadius");
  534 + m_minor.unit("deg");
  535 + m_minor.value(2.778e-4); // 1 arcsec
  536 + m_minor.min(2.778e-4); // 1 arcsec
  537 + m_minor.free();
  538 + m_minor.scale(1.0);
  539 + m_minor.gradient(0.0);
  540 + m_minor.hasgrad(false); // Elliptical components never have gradients
  541 +
  542 + m_major.clear();
  543 + m_major.name("MajorRadius");
  544 + m_major.unit("deg");
  545 + m_major.value(2.778e-4); // 1 arcsec
  546 + m_major.min(2.778e-4); // 1 arcsec
  547 + m_major.free();
  548 + m_major.scale(1.0);
  549 + m_major.gradient(0.0);
  550 + m_major.hasgrad(false); // Elliptical components never have gradients
  551 +
  552 + // Set parameter pointer(s)
  553 + m_pars.push_back(&m_minor);
  554 + m_pars.push_back(&m_major);
  555 +
  556 + // Initialise precomputation cache. Note that zero values flag
  557 + // uninitialised as a zero radius is not meaningful
  558 + m_last_minor = 0.0;
  559 + m_last_major = 0.0;
  560 + m_minor_rad = 0.0;
  561 + m_major_rad = 0.0;
  562 + m_norm = 0.0;
  563 +
  564 + // Return
  565 + return;
  566 +}
  567 +
  568 +
  569 +/***********************************************************************//**
  570 + * @brief Copy class members
  571 + *
  572 + * @param[in] model Elliptical disk model.
  573 + *
  574 + * We do not have to push back the members on the parameter stack as this
  575 + * should have been done by init_members() that was called before. Otherwise
  576 + * we would have the radius twice on the stack.
  577 + ***************************************************************************/
  578 +void GModelSpatialEllipticalDisk::copy_members(const GModelSpatialEllipticalDisk& model)
  579 +{
  580 + // Copy members
  581 + m_minor = model.m_minor;
  582 + m_major = model.m_major;
  583 +
  584 + // Copy precomputation cache
  585 + m_last_minor = model.m_last_minor;
  586 + m_last_major = model.m_last_major;
  587 + m_minor_rad = model.m_minor_rad;
  588 + m_major_rad = model.m_major_rad;
  589 + m_norm = model.m_norm;
  590 +
  591 + // Return
  592 + return;
  593 +}
  594 +
  595 +
  596 +/***********************************************************************//**
  597 + * @brief Delete class members
  598 + ***************************************************************************/
  599 +void GModelSpatialEllipticalDisk::free_members(void)
  600 +{
  601 + // Return
  602 + return;
  603 +}
  604 +
  605 +
  606 +/***********************************************************************//**
  607 + * @brief Update precomputation cache
  608 + *
  609 + * Computes the normalization
  610 + * \f[{\tt m\_norm} = \frac{1}{2 \pi (1 - \cos a) (1 - \cos b)}\f]
  611 + * TODO check this formula
  612 + *
  613 + * Note that this is the correct normalization on the sphere for any
  614 + * ellipse with a semi-minor axis a and semi-major axis b. For small a and b it is very similar to the
  615 + * cartesian
  616 + * approximation you might have expected:
  617 + * \f[{\tt m\_norm} = \frac{1}{\pi ab}\f]
  618 + ***************************************************************************/
  619 +void GModelSpatialEllipticalDisk::update() const
  620 +{
  621 + // Update if one axis has changed has changed
  622 + if (m_last_minor != minor() || m_last_major != major()) {
  623 +
  624 + // Store last values
  625 + m_last_minor = minor();
  626 + m_last_major = major();
  627 +
  628 + // Compute axes in radians
  629 + m_minor_rad = minor() * deg2rad;
  630 + m_major_rad = major() * deg2rad;
  631 +
  632 + // Perform precomputations
  633 + double denom = twopi * std::sqrt(1 - std::cos(m_minor_rad))*std::sqrt(1 - std::cos(m_major_rad));
  634 + m_norm = (denom > 0.0) ? 1.0 / denom : 0.0;
  635 +
  636 + } // endif: update required
  637 +
  638 + // Return
  639 + return;
  640 +}
... ...
src/model/Makefile.am
  View file @244ac5b
... ... @@ -20,6 +20,7 @@ sources = GModelPar.cpp \
20 20 GModelSpatialRadialGauss.cpp \
21 21 GModelSpatialRadialShell.cpp \
22 22 GModelSpatialElliptical.cpp \
  23 + GModelSpatialEllipticalDisk.cpp \
23 24 GModelSpatialDiffuse.cpp \
24 25 GModelSpatialDiffuseConst.cpp \
25 26 GModelSpatialDiffuseCube.cpp \
... ...