C++ GMatrix::eulerz方法代码示例

/***********************************************************************//**
 * @brief Rotate sky direction by zenith and azimuth angle
 *
 * @param[in] phi Azimuth angle (deg).
 * @param[in] theta Zenith angle (deg).
 *
 * Rotate sky direction by a zenith and azimuth angle given in the system
 * of the sky direction and aligned in celestial coordinates. 
 * The azimuth angle is counted counter clockwise from celestial north
 * (this is identical to the astronomical definition of a position angle).
 ***************************************************************************/
void GSkyDir::rotate_deg(const double& phi, const double& theta)
{
    // If we have no equatorial coordinates then get them now
    if (!m_has_radec && m_has_lb)
        gal2equ();

    // Allocate Euler and rotation matrices
    GMatrix ry;
    GMatrix rz;
    GMatrix rot;

    // Set up rotation matrix to rotate from native coordinates to
    // celestial coordinates
    ry.eulery(m_dec*rad2deg - 90.0);
    rz.eulerz(-m_ra*rad2deg);
    rot = transpose(ry * rz);

    // Set up native coordinate vector
    double phi_rad   = phi   * deg2rad;
    double theta_rad = theta * deg2rad;
    double cos_phi   = std::cos(phi_rad);
    double sin_phi   = std::sin(phi_rad);
    double cos_theta = std::cos(theta_rad);
    double sin_theta = std::sin(theta_rad);
    GVector native(-cos_phi*sin_theta, sin_phi*sin_theta, cos_theta);

    // Rotate vector into celestial coordinates
    GVector dir = rot * native;

    // Convert vector into sky position
    celvector(dir);

    // Return
    return;
}

/***********************************************************************//**
 * @brief Return spatially integrated data model
 *
 * @param[in] obsEng Measured event energy.
 * @param[in] obsTime Measured event time.
 * @param[in] obs Observation.
 * @return Spatially integrated model.
 *
 * @exception GException::invalid_argument
 *            No CTA event list found in observation.
 *            No CTA pointing found in observation.
 *
 * Spatially integrates the data model for a given measured event energy and
 * event time. This method also applies a deadtime correction factor, so that
 * the normalization of the model is a real rate (counts/exposure time).
 ***************************************************************************/
double GCTAModelBackground::npred(const GEnergy&      obsEng,
                                  const GTime&        obsTime,
                                  const GObservation& obs) const
{
    // Initialise result
    double npred     = 0.0;
    bool   has_npred = false;

    // Build unique identifier
    std::string id = obs.instrument() + "::" + obs.id();

    // Check if Npred value is already in cache
    #if defined(G_USE_NPRED_CACHE)
    if (!m_npred_names.empty()) {

        // Search for unique identifier, and if found, recover Npred value
		// and break
		for (int i = 0; i < m_npred_names.size(); ++i) {
			if (m_npred_names[i] == id && m_npred_energies[i] == obsEng) {
				npred     = m_npred_values[i];
				has_npred = true;
				#if defined(G_DEBUG_NPRED)
				std::cout << "GCTAModelBackground::npred:";
				std::cout << " cache=" << i;
				std::cout << " npred=" << npred << std::endl;
				#endif
				break;
			}
		}

    } // endif: there were values in the Npred cache
    #endif

    // Continue only if no Npred cache value was found
    if (!has_npred) {

        // Evaluate only if model is valid
        if (valid_model()) {

            // Get CTA event list
			const GCTAEventList* events = dynamic_cast<const GCTAEventList*>(obs.events());
            if (events == NULL) {
                std::string msg = "No CTA event list found in observation.\n" +
                                  obs.print();
                throw GException::invalid_argument(G_NPRED, msg);
            }

            #if !defined(G_NPRED_AROUND_ROI)
			// Get CTA pointing direction
			GCTAPointing* pnt = dynamic_cast<GCTAPointing*>(obs.pointing());
            if (pnt == NULL) {
                std::string msg = "No CTA pointing found in observation.\n" +
                                  obs.print();
                throw GException::invalid_argument(G_NPRED, msg);
            }
            #endif

            // Get reference to ROI centre
            const GSkyDir& roi_centre = events->roi().centre().dir();

			// Get ROI radius in radians
			double roi_radius = events->roi().radius() * gammalib::deg2rad;

			// Get distance from ROI centre in radians
            #if defined(G_NPRED_AROUND_ROI)
			double roi_distance = 0.0;
            #else
			double roi_distance = roi_centre.dist(pnt->dir());
            #endif

			// Initialise rotation matrix to transform from ROI system to
            // celestial coordinate system
			GMatrix ry;
			GMatrix rz;
			ry.eulery(roi_centre.dec_deg() - 90.0);
			rz.eulerz(-roi_centre.ra_deg());
			GMatrix rot = (ry * rz).transpose();

			// Compute position angle of ROI centre with respect to model
			// centre (radians)
            #if defined(G_NPRED_AROUND_ROI)
            double omega0 = 0.0;
            #else
			double omega0 = pnt->dir().posang(events->roi().centre().dir());
//.........这里部分代码省略.........