diff --git a/ChangeLog b/ChangeLog index d99f54c..c872871 100644 --- a/ChangeLog +++ b/ChangeLog @@ -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 diff --git a/NEWS b/NEWS index 5eaa763..28fdcf0 100644 --- a/NEWS +++ b/NEWS @@ -1,4 +1,4 @@ -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 diff --git a/inst/com/test/test_COM.cpp b/inst/com/test/test_COM.cpp index d053a99..29d55e6 100644 --- a/inst/com/test/test_COM.cpp +++ b/inst/com/test/test_COM.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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; diff --git a/inst/cta/test/test_CTA.cpp b/inst/cta/test/test_CTA.cpp index 02cc534..f417a8a 100644 --- a/inst/cta/test/test_CTA.cpp +++ b/inst/cta/test/test_CTA.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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, diff --git a/src/sky/GWcs.cpp b/src/sky/GWcs.cpp index 29482cd..231983b 100644 --- a/src/sky/GWcs.cpp +++ b/src/sky/GWcs.cpp @@ -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& 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;