Merge branch 'issue_1017_jk' into integration

Conflicts: ChangeLog

Merge branch 'issue_1017_jk' into integration
Conflicts: ChangeLog
Jürgen Knödlseder
2 parents 3c035092 9c75e11c
Showing 5 changed files with 161 additions and 28 deletions
ChangeLog
NEWS
inst/com/test/test_COM.cpp
inst/cta/test/test_CTA.cpp
src/sky/GWcs.cpp
@@ -24,6 +24,7 @@
         Add broken power law spectral model
         Add Gaussian spectral model
         Update parameter files if necessary
+        Implement exact sky map pixel solid angle computation
         Rename the following classes:
         - GWcs => GSkyProjection
         - GWcslib => GWcs
-New Features and Important Changes in GammaLib 00-08-00
+vNew Features and Important Changes in GammaLib 00-08-00
 Author(s) : Juergen Knoedlseder
@@ -565,7 +565,12 @@ The GSkyPixel class now supports 1D and 2D pixelisations, allowing
 a more generic usage of the class. Integer and double precision
 values can be assigned from and to GSkyPixel for 1D pixelisations.
-The GWcsAIT class has been added for Aitoff-Hammer projections.
+The GWcsAIT class has been added for Aitoff projections.
+
+The solid angle of sky map pixels are now computed using the spherical
+excess. This provides accurate solid angles even in case that the
+pixels are highly distorted, for example at the poles of the Aitoff
+projections.
 13. Support module
 /***************************************************************************
  *                       test_COM.cpp - test COM classes                   *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2012-2013 by Juergen Knoedlseder                         *
+ *  copyright (C) 2012-2014 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -580,7 +580,7 @@ void TestGCOMOptimize::test_binned_optimizer(void)
                             47, 0,
                             29.9617, 0.377506,
                             49, 0,
-                            27.9993, 0.372603};
+                            27.9992, 0.372603};
     test_try("Perform LM optimization");
     try {
         GOptimizerLM opt;
 /***************************************************************************
  *                       test_CTA.cpp - Test CTA classes                   *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010-2013 by Juergen Knoedlseder                         *
+ *  copyright (C) 2010-2014 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -335,7 +335,8 @@ void TestGCTAResponse::test_response_irf_diffuse(void)
 {
     // Set reference value
     //double ref = 13803.800313356;
-    const double ref = 13803.5186374;
+    //const double ref = 13803.5186374;
+    const double ref = 13796.8200974; // After GWcs::solidangle improvement
     // Set parameters
     double src_ra  = 201.3651;
@@ -412,7 +413,8 @@ void TestGCTAResponse::test_response_npred_diffuse(void)
 {
     // Set reference value
     //double ref = 11212.26274; // npred_spec precision of 1e-6
-    const double ref = 11212.437464; // npred_spec precision of 1e-5
+    //const double ref = 11212.437464; // npred_spec precision of 1e-5
+    const double ref = 11206.8536502; // After GWcs::solidangle improvement
     // Set parameters
     double src_ra  = 201.3651;
@@ -766,10 +768,10 @@ void TestGCTAOptimize::test_binned_optimizer(void)
     double fit_results[] = {83.6331, 0,
                             22.0145, 0,
                             5.616410411e-16, 1.904730785e-17,
-                            -2.481781246, -0.02580905077,
+                            -2.47925, -0.02580905077,
                             300000, 0,
                             1, 0,
-                            2.933677595, 0.06639644824,
+                            2.93353, 0.06639644824,
                             6.550723074e-05, 1.945714239e-06,
                             -1.833781187, -0.0160819,
                             1000000, 0,
@@ -56,6 +56,7 @@
 /* __ Coding definitions _________________________________________________ */
 //#define G_LIN_MATINV_FORCE_PC                             // Force PC usage
+#define G_USE_TRIGO_FOR_SOLIDANGLE           // Use trigonometric functions
 /* __ Debug definitions __________________________________________________ */
 //#define G_DIR2XY_DEBUG                                      // Debug dir2xy
@@ -425,34 +426,158 @@ void GWcs::write(GFitsHDU&amp; hdu) const
  *
  * @param[in] pixel Pixel index (x,y)
  *
- * Estimate solid angles of pixel by computing the coordinates in the 4 pixel
- * corners. The surface is computed using a cartesian approximation:
- *           a
- *     1-----5-----2                 a+b
- *    /      | h    \    where A = h ---
- *   4-------6-------3                2
- *           b
- * This is a brute force technique that works sufficiently well for non-
- * rotated sky maps. Something more intelligent should be implemented in
- * the future.
+ * Estimates solid angles of pixels using the Girard equation for excess
+ * area - see: http://mathworld.wolfram.com/SphericalPolygon.html
+ *
+ *
+ * Below, the definiton of the pixel cornes and sides are shown as used
+ * within the code.
+ *
+ *             a12
+ *         1---------2
+ *         |\       /|
+ *         | \a13  / |
+ *         |  \   /  |
+ *         |   \ /   |
+ *      a14|    X    |a23
+ *         |   / \   |
+ *         |  /   \  |
+ *         | /a24  \ |
+ *         |/       \|
+ *         4---------3
+ *             a34
+ *
  ***************************************************************************/
 double GWcs::solidangle(const GSkyPixel& pixel) const
 {
-    // Get the sky directions of the 6 points
+    // Get the sky directions of the 4 points
     GSkyDir dir1 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()-0.5));
     GSkyDir dir2 = pix2dir(GSkyPixel(pixel.x()+0.5, pixel.y()-0.5));
     GSkyDir dir3 = pix2dir(GSkyPixel(pixel.x()+0.5, pixel.y()+0.5));
     GSkyDir dir4 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()+0.5));
-    GSkyDir dir5 = pix2dir(GSkyPixel(pixel.x(),     pixel.y()-0.5));
-    GSkyDir dir6 = pix2dir(GSkyPixel(pixel.x(),     pixel.y()+0.5));
-    // Compute distances between sky directions
-    double a = dir1.dist(dir2);
-    double b = dir3.dist(dir4);
-    double h = dir5.dist(dir6);
+    // Compile option 1: Use triginometric function
+    #if defined(G_USE_TRIGO_FOR_SOLIDANGLE)
+
+    // Initialise solid angle
+    double solidangle = 0.0;
+
+    // Compute angular distances between pixel corners
+    double a12 = dir1.dist(dir2);
+    double a14 = dir1.dist(dir4);
+    double a23 = dir2.dist(dir3);
+    double a34 = dir3.dist(dir4);
+
+    // Special case: a12 or a14 is zero, then pixel is a triangle composed
+    // of [2,3,4]
+    if (a12 <= 0.0 || a14 <= 0.0) {
-    // Compute solid angle
-    double solidangle = 0.5*(h*(a+b));
+        // Compute diagonals
+        double a24 = dir2.dist(dir4);
+
+        // Compute sines
+        double sin_a23 = std::sin(a23);
+        double sin_a24 = std::sin(a24);
+        double sin_a34 = std::sin(a34);
+
+        // Compute cosines
+        double cos_a23 = std::cos(a23);
+        double cos_a24 = std::cos(a24);
+        double cos_a34 = std::cos(a34);
+
+        // Compute angles
+        double angle2 = std::acos((cos_a34-cos_a24*cos_a23)/(sin_a24*sin_a23));
+        double angle3 = std::acos((cos_a24-cos_a34*cos_a23)/(sin_a34*sin_a23));
+        double angle4 = std::acos((cos_a23-cos_a24*cos_a34)/(sin_a24*sin_a34));
+
+        // Compute excess area to determine solid angle
+        solidangle = (angle2 + angle3 + angle4) - gammalib::pi;
+
+    }
+
+    // Special case: a23 or a 34 is zero, then pixel is a triangle composed
+    // of [1,2,4]
+    else if (a23 <= 0.0 || a34 <= 0.0) {
+
+        // Compute diagonals
+        double a24 = dir2.dist(dir4);
+
+        // Compute sines
+        double sin_a12 = std::sin(a12);
+        double sin_a14 = std::sin(a14);
+        double sin_a24 = std::sin(a24);
+
+        // Compute cosines
+        double cos_a12 = std::cos(a12);
+        double cos_a14 = std::cos(a14);
+        double cos_a24 = std::cos(a24);
+
+        // Compute angles
+        double angle1 = std::acos((cos_a24-cos_a12*cos_a14)/(sin_a12*sin_a14));
+        double angle2 = std::acos((cos_a14-cos_a12*cos_a24)/(sin_a12*sin_a24));
+        double angle4 = std::acos((cos_a12-cos_a14*cos_a24)/(sin_a14*sin_a24));
+
+        // Compute excess area to determine solid angle
+        solidangle = (angle1 + angle2 + angle4) - gammalib::pi;
+
+    }
+
+    // Otherwise we have a polygon
+    else {
+
+        // Compute diagonals
+        double a13 = dir1.dist(dir3);
+        double a24 = dir2.dist(dir4);
+
+        // Compute sines
+        double sin_a12 = std::sin(a12);
+        double sin_a14 = std::sin(a14);
+        double sin_a23 = std::sin(a23);
+        double sin_a34 = std::sin(a34);
+
+        // Compute cosines
+        double cos_a12 = std::cos(a12);
+        double cos_a13 = std::cos(a13);
+        double cos_a14 = std::cos(a14);
+        double cos_a23 = std::cos(a23);
+        double cos_a24 = std::cos(a24);
+        double cos_a34 = std::cos(a34);
+
+        // Compute angles
+        double angle1 = std::acos((cos_a13-cos_a34*cos_a14)/(sin_a34*sin_a14));
+        double angle2 = std::acos((cos_a24-cos_a23*cos_a34)/(sin_a23*sin_a34));
+        double angle3 = std::acos((cos_a13-cos_a12*cos_a23)/(sin_a12*sin_a23));
+        double angle4 = std::acos((cos_a24-cos_a14*cos_a12)/(sin_a14*sin_a12));
+
+        // Use Girard equation for excess area to determine solid angle
+        solidangle = (angle1 + angle2 + angle3 + angle4) - gammalib::twopi;
+
+    } // endif: we had a polynom
+
+    // Compile option 2: Use vectors
+    #else
+
+    // Get vectors to pixel corners
+    GVector vec1 = dir1.celvector();
+    GVector vec2 = dir2.celvector();
+    GVector vec3 = dir3.celvector();
+    GVector vec4 = dir4.celvector();
+
+    // Compute inner angles of pixel corners
+    double angle1 = gammalib::acos(cross(vec2, (cross(vec1, vec2))) * 
+                                   cross(vec2, (cross(vec3, vec2))));
+    double angle2 = gammalib::acos(cross(vec3, (cross(vec2, vec3))) *
+                                   cross(vec3, (cross(vec4, vec3))));
+    double angle3 = gammalib::acos(cross(vec4, (cross(vec3, vec4))) *
+                                   cross(vec4, (cross(vec1, vec4))));
+    double angle4 = gammalib::acos(cross(vec1, (cross(vec4, vec1))) *
+                                   cross(vec1, (cross(vec2, vec1))));
+
+    // Use Girard equation for excess area to determine solid angle
+    double solidangle = gammalib::deg2rad * gammalib::deg2rad +
+                        (angle1 + angle2 + angle3 + angle4) -
+                        gammalib::twopi;
+    #endif
     // Return solid angle
     return solidangle;