com_helpers_response_vector.cpp

Go to the documentation of this file.

/***************************************************************************
 *                COMPTEL helper classes for vector response               *
 * ----------------------------------------------------------------------- *
 *  copyright (C) 2021 by Juergen Knoedlseder                              *
 * ----------------------------------------------------------------------- *
 *                                                                         *
 *  This program is free software: you can redistribute it and/or modify   *
 *  it under the terms of the GNU General Public License as published by   *
 *  the Free Software Foundation, either version 3 of the License, or      *
 *  (at your option) any later version.                                    *
 *                                                                         *
 *  This program is distributed in the hope that it will be useful,        *
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
 *  GNU General Public License for more details.                           *
 *                                                                         *
 *  You should have received a copy of the GNU General Public License      *
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.  *
 *                                                                         *
 ***************************************************************************/
/**
 * @file com_helpers_response_vector.cpp
 * @brief Implementation of COMPTEL helper classes for vector response
 * @author Juergen Knoedlseder
 */

/* __ Includes ___________________________________________________________ */
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <cmath>
#include "com_helpers_response_vector.hpp"
#include "GCOMResponse.hpp"
#include "GIntegral.hpp"
#include "GIntegrals.hpp"

/* __ Method name definitions ____________________________________________ */

/* __ Macros _____________________________________________________________ */

/* __ Coding definitions _________________________________________________ */

/* __ Debug definitions __________________________________________________ */

/* __ Constants __________________________________________________________ */


/*==========================================================================
 =                                                                         =
 =                             Helper methods                              =
 =                                                                         =
 ==========================================================================*/

/***********************************************************************//**
 * @brief Kernel for radial integration of radial models
 *
 * @param[in] rho Rho angle (radians).
 * @return Azimuthally integrated radial model.
 ***************************************************************************/
GVector com_radial_kerns_rho::eval(const double& rho)
{
    // Initialise kernel values
    m_irfs = 0.0;

    // Continue only if rho is positive
    if (rho > 0.0) {

        // Compute radial model value
        double model = m_model.eval(rho, m_bin->energy(), m_bin->time());

        // Continue only if model is positive
        if (model > 0.0) {

            // Precompute cosine and sine terms for azimuthal integration
            double sin_rho = std::sin(rho);
            double cos_rho = std::cos(rho);

            // Setup azimuthal integration kernel
            com_radial_kerns_omega integrands(m_iaq,
                                              m_irfs,
                                              m_bin,
                                              m_rot,
                                              m_drx,
                                              m_phigeo_bin_size,
                                              m_phigeo_bins,
                                              m_phibar_bins,
                                              sin_rho,
                                              cos_rho);

            // Setup integrator
            GIntegrals integral(&integrands);
            integral.fixed_iter(m_iter);

            // Integrate over Omega angle
            m_irfs = integral.romberg(0.0, gammalib::twopi, m_iter) *
                     model * sin_rho;

        } // endif: radial model was positive

    } // endif: phigeo was positive

    // Return kernel values
    return m_irfs;
}


/***********************************************************************//**
 * @brief Kernel for azimuthal integration of radial models
 *
 * @param[in] omega Omega angle (radians).
 * @return Kernel value for radial model.
 ***************************************************************************/
GVector com_radial_kerns_omega::eval(const double& omega)
{
    // Initialise kernel values
    m_irfs = 0.0;

    // Compute sine and cosine of azimuth angle
    double sin_omega = std::sin(omega);
    double cos_omega = std::cos(omega);

    // Get sky direction
    GVector native(-cos_omega*m_sin_rho, sin_omega*m_sin_rho, m_cos_rho);
    GVector dir = m_rot * native;
    GSkyDir skyDir;
    skyDir.celvector(dir);

    // Compute Phigeo of current position
    double phigeo = m_bin->dir().dir().dist(skyDir);

    // Precompute interpolated IAQ vector for Phigeo angle
    double phirat  = phigeo / m_phigeo_bin_size; // 0.5 at bin centre
    int    iphigeo = int(phirat);                // index into which Phigeo falls
    double eps     = phirat - iphigeo - 0.5;     // 0.0 at bin centre [-0.5, 0.5[

    // Continue only if Phigeo index is valid
    if (iphigeo < m_phigeo_bins) {

        // Get DRX value
        double intensity = m_drx(skyDir);

        // Multiply result with IAQ for all Phibar layers
        for (int iphibar = 0; iphibar < m_phibar_bins; ++iphibar) {

            // Initialise IAQ
            double iaq = 0.0;

            // 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 = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
                }
                else {
                    iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
                }
            }
            else {           // interpolate towards right
                if (iphigeo < m_phigeo_bins-1) {
                    iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
                }
                else {
                    iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
                }
            }

            // If interpolated IAQ value is positive then compute IRF value
            if (iaq > 0.0) {
                m_irfs[iphibar] = intensity * iaq;
            }

        } // endfor: looped over Phibar

    } // endif: Phigeo index was valid

    // Return kernel values
    return m_irfs;
}


/***********************************************************************//**
 * @brief Kernel for radial integration of elliptical models
 *
 * @param[in] rho Rho angle (radians).
 * @return Azimuthally integrated elliptical model.
 ***************************************************************************/
GVector com_elliptical_kerns_rho::eval(const double& rho)
{
    // Initialise kernel values
    m_irfs = 0.0;

    // Continue only if rho is positive
    if (rho > 0.0) {

        // Precompute cosine and sine terms for azimuthal integration
        double sin_rho = std::sin(rho);
        double cos_rho = std::cos(rho);

        // Setup azimuthal integration kernel
        com_elliptical_kerns_omega integrands(m_iaq,
                                              m_model,
                                              m_irfs,
                                              m_bin,
                                              m_rot,
                                              m_drx,
                                              m_phigeo_bin_size,
                                              m_phigeo_bins,
                                              m_phibar_bins,
                                              sin_rho,
                                              cos_rho);

        // Setup integrator
        GIntegrals integral(&integrands);
        integral.fixed_iter(m_iter);

        // Integrate over Omega angle
        m_irfs = integral.romberg(0.0, gammalib::twopi, m_iter) * sin_rho;

    } // endif: phigeo was positive

    // Return kernel values
    return m_irfs;
}


/***********************************************************************//**
 * @brief Kernel for azimuthal integration of elliptical models
 *
 * @param[in] omega Omega angle (radians).
 * @return Kernel value for elliptical model.
 ***************************************************************************/
GVector com_elliptical_kerns_omega::eval(const double& omega)
{
    // Initialise kernel values
    m_irfs = 0.0;

    // Compute sine and cosine of azimuth angle
    double sin_omega = std::sin(omega);
    double cos_omega = std::cos(omega);

    // Get sky direction
    GVector native(-cos_omega*m_sin_rho, sin_omega*m_sin_rho, m_cos_rho);
    GVector dir = m_rot * native;
    GSkyDir skyDir;
    skyDir.celvector(dir);

    // Compute Phigeo of current position
    double phigeo = m_bin->dir().dir().dist(skyDir);

    // Precompute interpolated IAQ vector for Phigeo angle
    double phirat  = phigeo / m_phigeo_bin_size; // 0.5 at bin centre
    int    iphigeo = int(phirat);                // index into which Phigeo falls
    double eps     = phirat - iphigeo - 0.5;     // 0.0 at bin centre [-0.5, 0.5[

    // Continue only if Phigeo index is valid
    if (iphigeo < m_phigeo_bins) {

        // Set photon
        GPhoton photon(skyDir, m_bin->energy(), m_bin->time());

        // Get model sky intensity for photon (unit: sr^-1)
        double intensity = m_model.spatial()->eval(photon);

        // Continue only if intensity is positive
        if (intensity > 0.0) {

            // Multiply-in DRX value
            intensity *= m_drx(skyDir);

            // Multiply result with IAQ for all Phibar layers
            for (int iphibar = 0; iphibar < m_phibar_bins; ++iphibar) {

                // Initialise IAQ
                double iaq = 0.0;

                // 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 = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
                    }
                    else {
                        iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
                    }
                }
                else {           // interpolate towards right
                    if (iphigeo < m_phigeo_bins-1) {
                        iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
                    }
                    else {
                        iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
                    }
                }

                // If interpolated IAQ value is positive then compute IRF value
                if (iaq > 0.0) {
                    m_irfs[iphibar] = intensity * iaq;
                }

            } // endfor: looped over Phibar

        } // endif: intensity was valid

    } // endif: Phigeo index was valid

    // Return kernel values
    return m_irfs;
}