Merge branch 'issue_1017' into issue_1017_jk

Conflicts: src/sky/GWcs.cpp

Merge branch 'issue_1017' into issue_1017_jk
Conflicts: src/sky/GWcs.cpp
Jürgen Knödlseder
2 parents d2abdb17 7d2423aa
Showing 3 changed files with 146 additions and 38 deletions
ChangeLog
NEWS
src/sky/GWcs.cpp
@@ -23,6 +23,7 @@
         Put likelihood computation in GObservation class
         Add broken power law spectral model
         Add Gaussian spectral model
+        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
  
@@ -560,7 +560,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
@@ -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,7 +426,27 @@ void GWcs::write(GFitsHDU&amp; hdu) const
  *
  * @param[in] pixel Pixel index (x,y)
  *
-TODO:: UPDATE THIS DOCSTRING (+link to wikipedia)
+ * 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
 {
@@ -435,47 +456,128 @@ double GWcs::solidangle(const GSkyPixel&amp; pixel) const
     GSkyDir dir3 = pix2dir(GSkyPixel(pixel.x()+0.5, pixel.y()+0.5));
     GSkyDir dir4 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()+0.5));
  
+    // 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 a13 = dir1.dist(dir3);
     double a14 = dir1.dist(dir4);
     double a23 = dir2.dist(dir3);
-    double a24 = dir2.dist(dir4);
     double a34 = dir3.dist(dir4);
  
-    // Angle 1
-    double a = a13;
-    double b = a34;
-    double c = a14;
-    double angle1 = std::acos((std::cos(a) - std::cos(b)*std::cos(c))/(std::sin(b)*std::sin(c)));
-
-    // Angle 2
-    a = a24;
-    b = a23;
-    c = a34;
-    double angle2 = std::acos((std::cos(a) - std::cos(b)*std::cos(c))/(std::sin(b)*std::sin(c)));
-
-    // Angle 3
-    a = a13;
-    b = a12;
-    c = a23;
-    double angle3 = std::acos((std::cos(a) - std::cos(b)*std::cos(c))/(std::sin(b)*std::sin(c)));
-
-    // Angle 4
-    a = a24;
-    b = a14;
-    c = a12;
-    double angle4 = std::acos((std::cos(a) - std::cos(b)*std::cos(c))/(std::sin(b)*std::sin(c)));
-
-    /*
-    std::cout << angle1 << " ";
-    std::cout << angle2 << " ";
-    std::cout << angle3 << " ";
-    std::cout << angle4 << " " << std::endl;
-    */
-    
-    // http://mathworld.wolfram.com/SphericalPolygon.html
-    double solidangle = (angle1 + angle2 + angle3 + angle4) - (2.0 * gammalib::pi);
+    // 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 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;