Sylvie Brau-Nogué / gammalib

Commit 96ba366ec5ba0b76b58f1d4fb2183848e36c963b

Authored by ellisowen
1 parent 6c5b8ca3

first attempt to improve the solid angle method in GWcs::solidangle

Inline Side-by-side
Showing 1 changed file with 13 additions and 19 deletions
src/sky/GWcs.cpp
  View file @96ba366
... ... @@ -425,34 +425,28 @@ void GWcs::write(GFitsHDU& hdu) const
425 425 *
426 426 * @param[in] pixel Pixel index (x,y)
427 427 *
428   - * Estimate solid angles of pixel by computing the coordinates in the 4 pixel
429   - * corners. The surface is computed using a cartesian approximation:
430   - * a
431   - * 1-----5-----2 a+b
432   - * / | h \ where A = h ---
433   - * 4-------6-------3 2
434   - * b
435   - * This is a brute force technique that works sufficiently well for non-
436   - * rotated sky maps. Something more intelligent should be implemented in
437   - * the future.
  428 +TODO:: UPDATE THIS DOCSTRING (+link to wikipedia)
438 429 ***************************************************************************/
439 430 double GWcs::solidangle(const GSkyPixel& pixel) const
440 431 {
441   - // Get the sky directions of the 6 points
  432 + // Get the sky directions of the 4 points
442 433 GSkyDir dir1 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()-0.5));
443 434 GSkyDir dir2 = pix2dir(GSkyPixel(pixel.x()+0.5, pixel.y()-0.5));
444 435 GSkyDir dir3 = pix2dir(GSkyPixel(pixel.x()+0.5, pixel.y()+0.5));
445 436 GSkyDir dir4 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()+0.5));
446   - GSkyDir dir5 = pix2dir(GSkyPixel(pixel.x(), pixel.y()-0.5));
447   - GSkyDir dir6 = pix2dir(GSkyPixel(pixel.x(), pixel.y()+0.5));
448 437  
449   - // Compute distances between sky directions
450   - double a = dir1.dist(dir2);
451   - double b = dir3.dist(dir4);
452   - double h = dir5.dist(dir6);
  438 + GVector vec1 = dir1.celvector();
  439 + GVector vec2 = dir2.celvector();
  440 + GVector vec3 = dir3.celvector();
  441 + GVector vec4 = dir4.celvector();
453 442  
454   - // Compute solid angle
455   - double solidangle = 0.5*(h*(a+b));
  443 + double angle1 = std::acos(cross(vec2, (cross(vec1, vec2))) * cross(vec2, (cross(vec3, vec2))));
  444 + double angle2 = std::acos(cross(vec3, (cross(vec2, vec3))) * cross(vec3, (cross(vec4, vec3))));
  445 + double angle3 = std::acos(cross(vec4, (cross(vec3, vec4))) * cross(vec4, (cross(vec1, vec4))));
  446 + double angle4 = std::acos(cross(vec1, (cross(vec4, vec1))) * cross(vec1, (cross(vec2, vec1))));
  447 +
  448 + // http://mathworld.wolfram.com/SphericalPolygon.html
  449 + double solidangle = (angle1 + angle2 + angle3 + angle4) - (2 * gammalib::pi);
456 450  
457 451 // Return solid angle
458 452 return solidangle;
... ...