Olga Sergijenko / gammalib

Commit db5e2b950466700ea02e64985908119d4cc59350

Authored by Andreas Specovius
Committed by Jürgen Knödlseder
1 parent 714a002b

Improved algorithm to transform between global<->local coordinate systems (now u… 

…se dist/posang). In terms of that, sorted out previous issues towards the poles. No longer need precomputed sin/cos, hence removed cache. Added 1e-10 tolerance for skydir containment check. Changed corners indexing, now starting with +RA+DEC and counting CCW. Improved documentation.
Inline Side-by-side
Showing 2 changed files with 40 additions and 67 deletions
include/GSkyRegionRect.hpp
  View file @db5e2b9
... ... @@ -49,7 +49,13 @@
49 49 * region and computed through internal method compute_solid().
50 50 *
51 51 * The position angle counts counter-clockwise from celestial North and is
52   - * aligned with the height axis of the rectangle.
  52 + * aligned with the height axis of the rectangle. Note, that height and width
  53 + * describe the diameter of the region.
  54 + *
  55 + * Sky directions can be transformed to the local/global coordinate system via
  56 + * the methods transform_to_local() / transform_to_global().
  57 + * The position (in global coordinates) of the 4 corners of the rectangle can
  58 + * be received via the method get_corner(index).
53 59 *
54 60 ***************************************************************************/
55 61 class GSkyRegionRect : public GSkyRegion {
... ... @@ -107,16 +113,12 @@ protected:
107 113 void copy_members(const GSkyRegionRect& region);
108 114 void free_members(void);
109 115 void compute_solid_angle(void);
110   - void update_cache(void);
111 116  
112 117 // Protected members
113 118 GSkyDir m_centre; //!< Centre or reference point of the region
114 119 double m_halfwidth; //!< Half width of the region [deg]
115 120 double m_halfheight; //!< Half height of the region [deg]
116 121 double m_posang; //!< Position angle, counterclockwise from North [radians]
117   -
118   - double m_posang_cos; //!< Cache for transform_to_local
119   - double m_posang_sin; //!< Cache for transform_to_local
120 122 };
121 123  
122 124  
... ...
src/sky/GSkyRegionRect.cpp
  View file @db5e2b9
... ... @@ -266,9 +266,6 @@ void GSkyRegionRect::posang(const double&amp; posang)
266 266 // Set the position angle
267 267 m_posang = posang;
268 268  
269   - // Update the cache
270   - update_cache();
271   -
272 269 // Return
273 270 return;
274 271 }
... ... @@ -586,10 +583,10 @@ bool GSkyRegionRect::contains_local(const GSkyDir&amp; locdir) const
586 583 bool dir_is_in = false;
587 584  
588 585 // Check containment: declination axis
589   - if (std::abs(locdir.dec_deg()) <= m_halfheight) {
  586 + if ((std::abs(locdir.dec_deg()) - 1.0e-10) <= m_halfheight) {
590 587  
591 588 // Check containment: right ascension axis
592   - if (std::abs(locdir.ra_deg()) <= m_halfwidth) {
  589 + if ((std::abs(locdir.ra_deg()) - 1.0e-10) <= m_halfwidth) {
593 590  
594 591 // Sky direction is inside the rectangle
595 592 dir_is_in = true;
... ... @@ -752,40 +749,22 @@ bool GSkyRegionRect::overlaps(const GSkyRegion&amp; reg) const
752 749 * Transform the sky direction :skydir: to the local cartesian coordinate system.
753 750 * The origin of the local coordinate system is fixed to the center of the
754 751 * rectangle and aligned in +ra (width) and +dec (height).
  752 + * Please note, that a GSkyDir object is returned, even if the local coordsys
  753 + * is cartesian.
755 754 ***************************************************************************/
756 755 GSkyDir GSkyRegionRect::transform_to_local(const GSkyDir& skydir) const
757 756 {
758   - // Get rid of negative right ascension coordinates
759   - // to obtain angles in [0,2pi]
760   - // (Surely there is a more clever way to to this)
761   - double skydir_ra = skydir.ra();
762   - if (skydir_ra < 0) {
763   - skydir_ra += gammalib::twopi;
764   - }
765   - double centre_ra = m_centre.ra();
766   - if (centre_ra < 0) {
767   - centre_ra += gammalib::twopi;
768   - }
769   -
770   - // Compute separation (in radians)
771   - double dx = skydir_ra - centre_ra;
772   - double dy = skydir.dec() - m_centre.dec();
773   -
774   - // Make sure the closest angular separation is chosen
775   - if (std::abs(dx) > gammalib::pi) {
776   - dx = gammalib::twopi - std::abs(dx);
777   - }
  757 + // Get distance and polar angle relative to rectangle centre
  758 + double dist = m_centre.dist(skydir);
  759 + double pa = m_centre.posang(skydir) - m_posang;
778 760  
779   - // Correction for spherical coordinate system
780   - dx *= std::cos(skydir.dec());
781   -
782   - // Rotate with neg. PA
783   - double new_x = m_posang_cos*dx - m_posang_sin*dy;
784   - double new_y = m_posang_sin*dx + m_posang_cos*dy;
  761 + // Compute corresponding x and y coordinate from polar coords
  762 + double x = dist * std::sin(pa);
  763 + double y = dist * std::cos(pa);
785 764  
786 765 // Create output sky direction
787 766 GSkyDir transformed;
788   - transformed.radec(new_x, new_y);
  767 + transformed.radec(x, y);
789 768  
790 769 // Return
791 770 return transformed;
... ... @@ -804,17 +783,25 @@ GSkyDir GSkyRegionRect::transform_to_local(const GSkyDir&amp; skydir) const
804 783 ***************************************************************************/
805 784 GSkyDir GSkyRegionRect::transform_to_global(const GSkyDir& locdir) const
806 785 {
807   - // Rotate local coordinates with positive PA
808   - double dx = (m_posang_cos *locdir.ra()) + (m_posang_sin *locdir.dec());
809   - double dy = (m_posang_cos *locdir.dec()) - (m_posang_sin *locdir.ra());
  786 + // // Initialise centre in local coordinates (maybe make class member?)
  787 + GSkyDir centre_local;
  788 + centre_local.radec(0,0);
810 789  
811   - // Shift by rectangle center, apply spherical coordinate scaling
812   - double new_y = m_centre.dec() + dy;
813   - double new_x = m_centre.ra() + dx/std::cos(new_y);
  790 + // Access x and y in local coordsys
  791 + double x = locdir.ra();
  792 + double y = locdir.dec();
814 793  
815   - // Create output sky direction
816   - GSkyDir transformed;
817   - transformed.radec(new_x, new_y);
  794 + // Compute dist and position angle in polar coords
  795 + double dist = std::sqrt(x*x + y*y);
  796 + double pa = gammalib::pihalf - std::atan2(y, x) + m_posang;
  797 +
  798 + // // Compute distance and position angle for locdir in global coordinates
  799 + // double dist = centre_local.dist(locdir);
  800 + // double pa = centre_local.posang(locdir) + m_posang;
  801 +
  802 + // Create output global sky direction
  803 + GSkyDir transformed(m_centre);
  804 + transformed.rotate(pa, dist);
818 805  
819 806 // Return
820 807 return transformed;
... ... @@ -829,6 +816,8 @@ GSkyDir GSkyRegionRect::transform_to_global(const GSkyDir&amp; locdir) const
829 816 *
830 817 * @exception GException::out_of_range
831 818 * Corner index is not in [0,3].
  819 + *
  820 + * The index counts CCW starting from +RA,+DEC: +R+D, +R-D, -R-D, -R+D.
832 821 ***************************************************************************/
833 822 GSkyDir GSkyRegionRect::get_corner(const int& index) const
834 823 {
... ... @@ -838,9 +827,9 @@ GSkyDir GSkyRegionRect::get_corner(const int&amp; index) const
838 827 }
839 828  
840 829 // Compute the offset to the corners
841   - // Will be tr, br, tl, bl for indices 0-3 respectively
842   - double dx = m_halfwidth * (1 - 2*int(index < 2));
843   - double dy = m_halfheight * (1 - 2*int(index % 2));
  830 + // Index counts CCW starting from +RA,+DEC: +R+D, +R-D, -R-D, -R+D
  831 + double dx = m_halfwidth * (1 - 2*int(index >= 2));
  832 + double dy = m_halfheight * (1 - 2*int((index==1) || (index==2)));
844 833  
845 834 // Define local corner sky direction
846 835 GSkyDir corner;
... ... @@ -891,10 +880,6 @@ void GSkyRegionRect::copy_members(const GSkyRegionRect&amp; region)
891 880 m_halfheight = region.m_halfheight;
892 881 m_posang = region.m_posang;
893 882  
894   - // Copy cache
895   - m_posang_cos = region.m_posang_cos;
896   - m_posang_sin = region.m_posang_sin;
897   -
898 883 // Return
899 884 return;
900 885 }
... ... @@ -922,17 +907,3 @@ void GSkyRegionRect::compute_solid_angle(void)
922 907 // Return
923 908 return;
924 909 }
925   -
926   -
927   -/***********************************************************************//**
928   - * @brief Update the cache.
929   - ***************************************************************************/
930   -void GSkyRegionRect::update_cache(void)
931   -{
932   - // Precompute cos and sin of position angle
933   - m_posang_cos = std::cos(m_posang);
934   - m_posang_sin = std::sin(m_posang);
935   -
936   - // Return
937   - return;
938   -}
... ...