From 5ce529679a40ebdd5590abd4a9dc6423d31c869a Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Sat, 16 Oct 2021 22:08:48 +0200 Subject: [PATCH] Optimise computation of COMPTEL response for extended sources --- inst/com/src/com_helpers_response_vector.cpp | 87 +++++++++++++++++++++++++++++++-------------------------------------------------------- inst/com/src/com_helpers_response_vector.hpp | 12 +++--------- 2 files changed, 34 insertions(+), 65 deletions(-) diff --git a/inst/com/src/com_helpers_response_vector.cpp b/inst/com/src/com_helpers_response_vector.cpp index 628e3d8..d27e388 100644 --- a/inst/com/src/com_helpers_response_vector.cpp +++ b/inst/com/src/com_helpers_response_vector.cpp @@ -31,6 +31,7 @@ #include #include "com_helpers_response_vector.hpp" #include "GCOMResponse.hpp" +#include "GIntegral.hpp" #include "GIntegrals.hpp" /* __ Method name definitions ____________________________________________ */ @@ -58,6 +59,9 @@ ***************************************************************************/ GVector com_extended_kerns_phigeo::eval(const double& phigeo) { + // Initialise kernel values + m_irfs = 0.0; + // Continue only if phigeo is positive if (phigeo > 0.0) { @@ -90,56 +94,54 @@ GVector com_extended_kerns_phigeo::eval(const double& phigeo) // Continue only if Phigeo index is valid if (iphigeo < m_phigeo_bins) { - // Initialise IAQ vector to hold precomputation results - GVector iaq(m_phibar_bins); + // Setup kernel for azimuthal integration + com_extended_kerns_phi integrand(m_model, + m_srcEng, + m_srcTime, + m_rot, + m_drx, + sin_phigeo, + cos_phigeo); + + // Setup integrator + GIntegral integral(&integrand); + integral.fixed_iter(m_iter); + + // Integrate over azimuth angle + double irf = integral.romberg(phi_min, phi_max, m_iter) * sin_phigeo; - // Loop over Phibar + // Multiply result with IAQ for all Phibar layers for (int iphibar = 0; iphibar < m_phibar_bins; ++iphibar) { + // Initialise IAQ + double iaq; + // Get IAQ index int i = iphibar * m_phigeo_bins + iphigeo; // Get interpolated IAQ value if (eps < 0.0) { // interpolate towards left if (iphigeo > 0) { - iaq[iphibar] = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1]; + iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1]; } else { - iaq[iphibar] = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1]; + iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1]; } } else { // interpolate towards right if (iphigeo < m_phigeo_bins-1) { - iaq[iphibar] = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1]; + iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1]; } else { - iaq[iphibar] = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1]; + iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1]; } } - // Normalise IAQ value - iaq[iphibar] *= m_iaq_norm; + // Compute IRF + m_irfs[iphibar] = irf * iaq * m_iaq_norm; } // endfor: looped over Phibar - // Setup kernel for azimuthal integration - com_extended_kerns_phi integrands(m_model, - m_irfs, - m_srcEng, - m_srcTime, - m_rot, - m_drx, - iaq, - sin_phigeo, - cos_phigeo); - - // Setup integrator - GIntegrals integral(&integrands); - integral.fixed_iter(m_iter); - - // Integrate over azimuth - m_irfs = integral.romberg(phi_min, phi_max, m_iter) * sin_phigeo; - } // endif: Phigeo bin was valid } // endif: arc length was positive @@ -157,11 +159,8 @@ GVector com_extended_kerns_phigeo::eval(const double& phigeo) * @param[in] phi Azimuth angle (radians). * @return Vector of azimuthally integrated model. ***************************************************************************/ -GVector com_extended_kerns_phi::eval(const double& phi) +double com_extended_kerns_phi::eval(const double& phi) { - // Initialise kernel values - m_irfs = 0.0; - // Compute sine and cosine of azimuth angle double sin_phi = std::sin(phi); double cos_phi = std::cos(phi); @@ -186,32 +185,8 @@ GVector com_extended_kerns_phi::eval(const double& phi) intensity *= (*m_drx)(skyDir); } - // Continue only if intensity is still positive - if (intensity > 0.0) { - - // Get number of Phibar layers - int nphibar = m_iaq.size(); - - // Loop over Phibar - for (int iphibar = 0; iphibar < nphibar; ++iphibar) { - - // Continue only if IAQ value is positive - if (m_iaq[iphibar] > 0.0) { - - // Compute IRF value (unit: cm^2) - double irf = m_iaq[iphibar] * intensity; - - // Set IRF value - m_irfs[iphibar] = irf; - - } // endif: IAQ value was valid - - } // endfor: looped over Phibar - - } // endif: intensity was valid - } // endif: intensity was valid // Return kernel values - return m_irfs; + return intensity; } diff --git a/inst/com/src/com_helpers_response_vector.hpp b/inst/com/src/com_helpers_response_vector.hpp index ed7187c..35142ba 100644 --- a/inst/com/src/com_helpers_response_vector.hpp +++ b/inst/com/src/com_helpers_response_vector.hpp @@ -32,6 +32,7 @@ #include "GTime.hpp" #include "GVector.hpp" #include "GMatrix.hpp" +#include "GFunction.hpp" #include "GFunctions.hpp" #include "GModelSky.hpp" #include "GSkyMap.hpp" @@ -121,36 +122,29 @@ protected: * This class provides the kernel for the azimuth angle integration of * extended sky models. ***************************************************************************/ -class com_extended_kerns_phi : public GFunctions { +class com_extended_kerns_phi : public GFunction { public: com_extended_kerns_phi(const GModelSky& model, - GVector& irfs, const GEnergy& srcEng, const GTime& srcTime, const GMatrix& rot, const GSkyMap* drx, - const GVector& iaq, const double& sin_phigeo, const double& cos_phigeo) : m_model(model), - m_irfs(irfs), m_srcEng(srcEng), m_srcTime(srcTime), m_rot(rot), m_drx(drx), - m_iaq(iaq), m_sin_phigeo(sin_phigeo), m_cos_phigeo(cos_phigeo) { } - int size(void) const { return m_irfs.size(); } - GVector eval(const double& phi); + double eval(const double& phi); protected: const GModelSky& m_model; //!< Sky model - GVector& m_irfs; //!< IRF vector to update const GEnergy& m_srcEng; //!< Source energy const GTime& m_srcTime; //!< Source time const GMatrix& m_rot; //!< Rotation matrix const GSkyMap* m_drx; //!< DRX - const GVector m_iaq; //!< Precomputed IAQ vector const double& m_sin_phigeo; //!< Sine of Phigeo const double& m_cos_phigeo; //!< Cosine of Phigeo }; -- libgit2 0.21.4