Justine Devin / gammalib

Commit 8e42fabd875d7ee468c3aece01be66af1f7df1dc

Authored by Jürgen Knödlseder
1 parent 96ba366e

Compute polygon inner angles in GWcs::solidangle method using spherical triangles.

Inline Side-by-side
Showing 1 changed file with 38 additions and 2 deletions
src/sky/GWcs.cpp
  View file @8e42fab
... ... @@ -435,6 +435,7 @@ double GWcs::solidangle(const GSkyPixel& pixel) const
435 435 GSkyDir dir3 = pix2dir(GSkyPixel(pixel.x()+0.5, pixel.y()+0.5));
436 436 GSkyDir dir4 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()+0.5));
437 437  
  438 + /*
438 439 GVector vec1 = dir1.celvector();
439 440 GVector vec2 = dir2.celvector();
440 441 GVector vec3 = dir3.celvector();
... ... @@ -444,9 +445,44 @@ double GWcs::solidangle(const GSkyPixel& pixel) const
444 445 double angle2 = std::acos(cross(vec3, (cross(vec2, vec3))) * cross(vec3, (cross(vec4, vec3))));
445 446 double angle3 = std::acos(cross(vec4, (cross(vec3, vec4))) * cross(vec4, (cross(vec1, vec4))));
446 447 double angle4 = std::acos(cross(vec1, (cross(vec4, vec1))) * cross(vec1, (cross(vec2, vec1))));
447   -
  448 + */
  449 +
  450 + //
  451 + double a;
  452 + double b;
  453 + double c;
  454 +
  455 + // Angle 1
  456 + a = dir1.dist(dir3);
  457 + b = dir4.dist(dir3);
  458 + c = dir4.dist(dir1);
  459 + double angle1 = std::acos((std::cos(a) - std::cos(b)*std::cos(c))/(std::sin(b)*std::sin(c)));
  460 +
  461 + // Angle 2
  462 + a = dir4.dist(dir2);
  463 + b = dir2.dist(dir3);
  464 + c = dir3.dist(dir4);
  465 + double angle2 = std::acos((std::cos(a) - std::cos(b)*std::cos(c))/(std::sin(b)*std::sin(c)));
  466 +
  467 + // Angle 3
  468 + a = dir1.dist(dir3);
  469 + b = dir1.dist(dir2);
  470 + c = dir2.dist(dir3);
  471 + double angle3 = std::acos((std::cos(a) - std::cos(b)*std::cos(c))/(std::sin(b)*std::sin(c)));
  472 +
  473 + // Angle 4
  474 + a = dir4.dist(dir2);
  475 + b = dir4.dist(dir1);
  476 + c = dir1.dist(dir2);
  477 + double angle4 = std::acos((std::cos(a) - std::cos(b)*std::cos(c))/(std::sin(b)*std::sin(c)));
  478 +
  479 + std::cout << angle1 << " ";
  480 + std::cout << angle2 << " ";
  481 + std::cout << angle3 << " ";
  482 + std::cout << angle4 << " " << std::endl;
  483 +
448 484 // http://mathworld.wolfram.com/SphericalPolygon.html
449   - double solidangle = (angle1 + angle2 + angle3 + angle4) - (2 * gammalib::pi);
  485 + double solidangle = (angle1 + angle2 + angle3 + angle4) - (2.0 * gammalib::pi);
450 486  
451 487 // Return solid angle
452 488 return solidangle;
... ...