first go at flux from circular region, does not work for small sources

Luigi Tibaldo
1 parent 4518c159
Showing 3 changed files with 241 additions and 0 deletions
include/GModelSpatial.hpp
pyext/GModelSpatial.i
src/model/GModelSpatial.cpp
@@ -34,10 +34,13 @@
 #include "GModelPar.hpp"
 #include "GPhoton.hpp"
 #include "GSkyDir.hpp"
+#include "GSkyRegionCircle.hpp"
 #include "GEnergy.hpp"
 #include "GTime.hpp"
 #include "GXmlElement.hpp"
 #include "GRan.hpp"
+#include "GFunction.hpp"
+#include "GIntegral.hpp"
  
 /* __ Forward declarations _______________________________________________ */
 class GSkyRegion;
@@ -102,6 +105,7 @@ public:
     virtual void           write(GXmlElement& xml) const = 0;
     virtual GSkyRegion*    region(void) const = 0;
     virtual std::string    print(const GChatter& chatter = NORMAL) const = 0;
+				      
  
     // Methods
     GModelPar&       at(const int& index);
@@ -110,6 +114,9 @@ public:
     bool             has_free_pars(void) const;
     int              size(void) const;
     void             autoscale(void);
+    double           flux(const GSkyRegionCircle& reg,
+  			  const GEnergy& srcEng = GEnergy(),
+                          const GTime&   srcTime = GTime()) const;
  
 protected:
     // Protected methods
@@ -117,6 +124,87 @@ protected:
     void copy_members(const GModelSpatial& model);
     void free_members(void);
  
+    /***********************************************************************//**
+ * @class circle_int_kern_rho
+ *
+ * @brief Kernel for circular sky region radial integration
+ *
+ * This class implements the integration kernel \f$K(\rho)\f$ for the
+ * integration
+ *
+ * \f[
+ *    \int_{\rho_{\rm min}}^{\rho_{\rm max}} K(\rho | E, t) d\rho
+ * \f]
+ *
+ * of a spatial model over a circular region. The eval() method computes
+ *
+ * \f[
+ *    K(\rho | E, t) = \sin \rho \times
+ *                     \int_{\omega} M(\rho, \omega | E, t) d\omega
+ * \f]
+ *
+ * where
+ * \f$M(\rho, \omega | E, t)\f$ is the spatial model,
+ * \f$\rho\f$ is the distance from the region centre, and
+ * \f$\omega\f$ is the position angle with respect to the connecting line
+ * between the region centre and the direction on the sky.
+ ***************************************************************************/
+  class circle_int_kern_rho : public GFunction {
+  public:
+    circle_int_kern_rho(const GModelSpatial*      model,
+			  const GSkyRegionCircle& reg,
+			  const GEnergy&          srcEng,
+			  const GTime&            srcTime) :
+      m_model(model),
+      m_reg(reg),
+      m_srcEng(srcEng),
+      m_srcTime(srcTime) { }
+    double eval(const double& rho);
+  public:
+    const GModelSpatial*       m_model;   //!< Model
+    GSkyRegionCircle    m_reg;     //!< Region
+    GEnergy                    m_srcEng;  //!< Photon energy
+    GTime                      m_srcTime; //!< Photon time
+  };
+
+  /***********************************************************************//**
+ * @class circle_int_kern_omega
+ *
+ * @brief Kernel for circular sky region azimuth angle integration
+ *
+ * This class implements the computation of
+ *
+ * \f[
+ *    K(\omega | \rho, E, t) = M(\omega | \rho)
+ * \f]
+ *
+ * where
+ * \fM(\omega | \rho)\f$ is the spatial model,
+ * \f$\rho\f$ is the distance from the region centre, and
+ * \f$\omega\f$ is the position angle with respect to the connecting line
+ * between the region centre and the direction on the sky.
+ ***************************************************************************/
+  class circle_int_kern_omega : public GFunction {
+  public:
+    circle_int_kern_omega(const GModelSpatial*    model,
+			  const GSkyRegionCircle& reg,
+			  const double&           rho,
+			  const GEnergy&          srcEng,
+			  const GTime&            srcTime) :
+      m_model(model),
+      m_reg(reg),
+      m_rho(rho),
+      m_srcEng(srcEng),
+      m_srcTime(srcTime) { }
+    double eval(const double& omega);
+  public:
+    const GModelSpatial*       m_model;   //!< Model
+    GSkyRegionCircle     m_reg;     //!< Region
+    double                     m_rho;     //!< Offset from center of the region
+    GEnergy                    m_srcEng;  //!< Photon energy
+    GTime                      m_srcTime; //!< Photon time
+  };
+
     // Proteced members
     std::vector<GModelPar*> m_pars;   //!< Parameter pointers
 };
@@ -64,6 +64,9 @@ public:
     bool       has_free_pars(void) const;
     int        size(void) const;
     void       autoscale(void);
+    double     flux(const GSkyRegionCircle& reg,
+    		    const GEnergy& srcEng = GEnergy(),
+    		    const GTime&   srcTime = GTime()) const;
 };
  
  
@@ -307,6 +307,57 @@ void GModelSpatial::autoscale(void)
     return;
 }
  
+/***********************************************************************//**
+ * @brief Returns model flux integrated in circular sky region
+ *
+ * @param[in] reg Sky region.
+ * @param[in] srcEng Energy.
+ * @param[in] srcTime Time.
+ * @return flux (adimensional or or ph/cm2/s).
+ *
+ * Computes
+ *
+ * \f[
+ *    \int_{\rho_{\rm min}}^{\rho_{\rm max}}  \sin \rho \times
+ *                     \int_{\omega} M(\rho, \omega | E, t) d\omega d\rho
+ * \f]
+ *
+ * where
+ * \f$M(\rho, \omega | E, t)\f$ is the spatial model,
+ * \f$\rho\f$ is the distance from the region centre, and
+ * \f$\omega\f$ is the position angle with respect to the connecting line
+ * between the region centre and the direction on the sky.
+ ***************************************************************************/
+double GModelSpatial::flux(const GSkyRegionCircle& reg,
+  			   const GEnergy& srcEng,
+  			   const GTime&   srcTime) const
+{
+    // Set number of iterations for Romberg integration.
+    // Need to test!!!!
+    static const int iter_rho = 5;
+    
+    // Initialise flux
+    double flux = 0.0;
+
+    // Define radial integration boundaries from 0 to ROI radius (in radians)
+    double rho_min = 0.;
+    double rho_max = reg.radius() * gammalib::deg2rad;
+
+    // Setup integration kernel
+    GModelSpatial::circle_int_kern_rho integrand(this,
+                                                 reg,
+    				                 srcEng,
+    				                 srcTime);
+    GIntegral integral(&integrand);
+
+    // Perform integration
+    flux = integral.romberg(rho_min, rho_max);
+
+    // Return
+    return flux;
+
+}
+
  
 /*==========================================================================
  =                                                                         =
@@ -350,3 +401,102 @@ void GModelSpatial::free_members(void)
     // Return
     return;
 }
+
+/***********************************************************************//**
+ * @brief Kernel for circular sky region radial integration
+ *
+ * @param[in] rho Radial distance from region centre (radians).
+ * @return Integration kernel defined as
+ *
+ * \f[
+ *    \int_{\rho_{\rm min}}^{\rho_{\rm max}} K(\rho | E, t) d\rho
+ * \f]
+ *
+ * of a spatial model over a circular region. The eval() method computes
+ *
+ * \f[
+ *    K(\rho | E, t) = \sin \rho \times
+ *                     \int_{\omega} M(\rho, \omega | E, t) d\omega
+ * \f]
+ *
+ * where
+ * \f$M(\rho, \omega | E, t)\f$ is the spatial model,
+ * \f$\rho\f$ is the distance from the region centre, and
+ * \f$\omega\f$ is the position angle with respect to the connecting line
+ * between the region centre and the direction on the sky.
+ ***************************************************************************/
+double GModelSpatial::circle_int_kern_rho::eval(const double& rho)
+{
+    // Initialise model value
+    double flux = 0.0;
+
+     // Continue only if rho is positive
+    if (rho > 0.0) {
+      
+      // Integrate over entire circular region
+      // omega between 0 and 2pi
+      double omega_min = 0.;
+      double omega_max = 2 * gammalib::pi;
+
+      // Setup integration kernel for azimuth integration
+      GModelSpatial::circle_int_kern_omega integrand(m_model,
+					             m_reg,
+						     rho,
+					             m_srcEng,
+					             m_srcTime);
+
+     // Setup integrator
+     GIntegral integral(&integrand);
+
+     // Compute sine term for radial integration
+     double sin_rho = std::sin(rho);
+
+     // Integrate over omega
+     flux = integral.romberg(omega_min, omega_max) * sin_rho;
+
+    } // endif: rho was positive
+
+      
+     // Return
+    return flux;
+}
+
+/***********************************************************************//**
+ * @brief Kernel for circular sky region azimuth angle integration
+ *
+ * @param[in] omega Azimuth angle (radians).
+* Computes
+ *
+ * \f[
+ *    K(\omega | \rho, E, t) = M(\omega | \rho)
+ * \f]
+ *
+ * where
+ * \fM(\omega | \rho)\f$ is the spatial model,
+ * \f$\rho\f$ is the distance from the region centre, and
+ * \f$\omega\f$ is the position angle with respect to the connecting line
+ * between the region centre and the direction on the sky.
+ ***************************************************************************/
+double GModelSpatial::circle_int_kern_omega::eval(const double& omega)
+{
+    // Initialise model value
+    double flux = 0.0;
+
+    // Precompute rho and omega in degrees
+    double omega_deg = omega * gammalib::rad2deg;
+    double rho_deg = m_rho * gammalib::rad2deg;
+
+    // Compute sky direction corresponding to given rho, omega
+    // Deep copy of region centre
+    GSkyDir dir = GSkyDir(m_reg.centre());
+    // Rotate to otain requested sky direction
+    dir.rotate_deg(omega_deg,rho_deg);
+
+    // Evaluate model for this sky direction
+    GPhoton photon = GPhoton(dir,m_srcEng,m_srcTime);
+    flux = m_model->eval(photon);
+
+     // Return
+    return flux;
+  
+}