Justine Devin / gammalib

Commit 4e9565bd2393541dd0249f0cdfa13d49fd6886f4

Authored by Jürgen Knödlseder
1 parent 9aae907d

Correct solidangle computation for triangular pixels.

Inline Side-by-side
Showing 1 changed file with 9 additions and 23 deletions
src/sky/GWcs.cpp
  View file @4e9565b
... ... @@ -564,46 +564,32 @@ double GWcs::solidangle(const GSkyPixel& pixel) const
564 564 double a12 = dir1.dist(dir2);
565 565 double a14 = dir1.dist(dir4);
566 566 double a23 = dir2.dist(dir3);
  567 + double a24 = dir2.dist(dir4);
567 568 double a34 = dir3.dist(dir4);
568 569  
569 570 // Special case: a12 or a14 is zero, then pixel is a triangle composed
570 571 // of [2,3,4]
571 572 if (a12 <= 0.0 || a14 <= 0.0) {
572   -
573   - // Compute diagonals
574   - double a24 = dir2.dist(dir4);
575   -
576   - // Triangle [2,3,4]
577 573 double s = 0.5 * (a23 + a34 + a24);
578   - solidangle = std::atan(std::sqrt(std::tan(0.5*s) *
579   - std::tan(0.5*(s-a23)) *
580   - std::tan(0.5*(s-a34)) *
581   - std::tan(0.5*(s-a24))));
582   -
  574 + solidangle = 4.0 * std::atan(std::sqrt(std::tan(0.5*s) *
  575 + std::tan(0.5*(s-a23)) *
  576 + std::tan(0.5*(s-a34)) *
  577 + std::tan(0.5*(s-a24))));
583 578 }
584 579  
585 580 // Special case: a23 or a 34 is zero, then pixel is a triangle composed
586 581 // of [1,2,4]
587 582 else if (a23 <= 0.0 || a34 <= 0.0) {
588   -
589   - // Compute diagonals
590   - double a24 = dir2.dist(dir4);
591   -
592   - // Triangle [1,2,4]
593 583 double s = 0.5 * (a12 + a24 + a14);
594   - solidangle = std::atan(std::sqrt(std::tan(0.5*s) *
595   - std::tan(0.5*(s-a12)) *
596   - std::tan(0.5*(s-a24)) *
597   - std::tan(0.5*(s-a14))));
598   -
  584 + solidangle = 4.0 * std::atan(std::sqrt(std::tan(0.5*s) *
  585 + std::tan(0.5*(s-a12)) *
  586 + std::tan(0.5*(s-a24)) *
  587 + std::tan(0.5*(s-a14))));
599 588 }
600 589  
601 590 // Otherwise we have a polygon
602 591 else {
603 592  
604   - // Compute diagonal
605   - double a24 = dir2.dist(dir4);
606   -
607 593 // Triangle 1 [1,2,4]
608 594 double s1 = 0.5 * (a12 + a24 + a14);
609 595 double excess1 = std::atan(std::sqrt(std::tan(0.5*s1) *
... ...