GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GWcs.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GWcs.cpp - Abstract world coordinate system base class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2011-2019 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GWcs.cpp
23  * @brief Abstract world coordinate system base class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <cstdlib>
32 #include <cmath>
33 #include "GException.hpp"
34 #include "GTools.hpp"
35 #include "GMath.hpp"
36 #include "GWcs.hpp"
37 
38 /* __ Method name definitions ____________________________________________ */
39 #define G_READ "GWcs::read(GFitsHDU&)"
40 #define G_PIX2DIR "GWcs::pix2dir(int&)"
41 #define G_DIR2PIX "GWcs::dir2pix(GSkyDir&)"
42 #define G_CRVAL "GWcs::crval(int&)"
43 #define G_CRPIX "GWcs::crpix(int&)"
44 #define G_CDELT "GWcs::cdelt(int&)"
45 #define G_WCS_SET_CTYPE "GWcs::wcs_set_ctype()"
46 #define G_WCS_P2S "GWcs::wcs_s2p(int, int, double*, double*, double*,"\
47  " double*, double*, int*)"
48 #define G_WCS_S2P "GWcs::wcs_s2p(int, int, double*, double*, double*,"\
49  " double*, double*, int*)"
50 #define G_CEL_SET "GWcs::cel_set()"
51 #define G_LIN_MATINV "GWcs::lin_matinv(std::vector<double>&,"\
52  " std::vector<double>&)"
53 
54 
55 /* __ Macros _____________________________________________________________ */
56 
57 /* __ Coding definitions _________________________________________________ */
58 //#define G_LIN_MATINV_FORCE_PC // Force PC usage
59 #define G_SOLIDANGLE_OPTION 2 // 1=sin/cos, 2=tan, 3=vectors
60 
61 /* __ Debug definitions __________________________________________________ */
62 //#define G_DIR2XY_DEBUG // Debug dir2xy
63 //#define G_XY2DIR_DEBUG // Debug xy2dir
64 
65 /* __ Local prototypes ___________________________________________________ */
66 
67 /* __ Constants __________________________________________________________ */
68 const double GWcs::UNDEFINED = 987654321.0e99;
69 
70 
71 /*==========================================================================
72  = =
73  = Constructors/destructors =
74  = =
75  ==========================================================================*/
76 
77 /***********************************************************************//**
78  * @brief Void constructor
79  ***************************************************************************/
81 {
82  // Initialise class members
83  init_members();
84 
85  // Return
86  return;
87 }
88 
89 
90 /***********************************************************************//**
91  * @brief Standard WCS sky map constructor
92  *
93  * @param[in] coords Coordinate system.
94  * @param[in] crval1 X value of reference pixel.
95  * @param[in] crval2 Y value of reference pixel.
96  * @param[in] crpix1 X index of reference pixel (first pixel is 1).
97  * @param[in] crpix2 Y index of reference pixel (first pixel is 1).
98  * @param[in] cdelt1 Increment in x direction at reference pixel (deg).
99  * @param[in] cdelt2 Increment in y direction at reference pixel (deg).
100  *
101  * Construct standard WCS sky map from standard definition parameters.
102  ***************************************************************************/
103 GWcs::GWcs(const std::string& coords,
104  const double& crval1, const double& crval2,
105  const double& crpix1, const double& crpix2,
106  const double& cdelt1, const double& cdelt2) : GSkyProjection()
107 {
108  // Initialise class members
109  init_members();
110 
111  // Set standard parameters
112  set_members(coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2);
113 
114  // Return
115  return;
116 }
117 
118 
119 /***********************************************************************//**
120  * @brief Construct from FITS HDU table
121  *
122  * @param[in] hdu FITS HDU.
123  ***************************************************************************/
125 {
126  // Initialise class members
127  init_members();
128 
129  // Read WCS definition from FITS HDU
130  read(hdu);
131 
132  // Return
133  return;
134 }
135 
136 
137 /***********************************************************************//**
138  * @brief Copy constructor
139  *
140  * @param wcs World Coordinate System.
141  ***************************************************************************/
142 GWcs::GWcs(const GWcs& wcs) : GSkyProjection(wcs)
143 {
144  // Initialise class members for clean destruction
145  init_members();
146 
147  // Copy members
148  copy_members(wcs);
149 
150  // Return
151  return;
152 }
153 
154 
155 /***********************************************************************//**
156  * @brief Destructor
157  ***************************************************************************/
159 {
160  // Free members
161  free_members();
162 
163  // Return
164  return;
165 }
166 
167 
168 /*==========================================================================
169  = =
170  = Operators =
171  = =
172  ==========================================================================*/
173 
174 /***********************************************************************//**
175  * @brief Assignment operator
176  *
177  * @param[in] wcs World Coordinate System.
178  * @return World Coordinate System.
179  ***************************************************************************/
181 {
182  // Execute only if object is not identical
183  if (this != &wcs) {
184 
185  // Copy base class members
186  this->GSkyProjection::operator=(wcs);
187 
188  // Free members
189  free_members();
190 
191  // Initialise private members for clean destruction
192  init_members();
193 
194  // Copy members
195  copy_members(wcs);
196 
197  } // endif: object was not identical
198 
199  // Return this object
200  return *this;
201 }
202 
203 
204 /*==========================================================================
205  = =
206  = Public methods =
207  = =
208  ==========================================================================*/
209 
210 /***********************************************************************//**
211  * @brief Read WCS definition from FITS header.
212  *
213  * @param[in] hdu FITS HDU.
214  *
215  * @exception GException::wcs_bad_coords
216  * Coordinate system is not of valid type.
217  *
218  * This method reads the WCS definition from the FITS header.
219  ***************************************************************************/
220 void GWcs::read(const GFitsHDU& hdu)
221 {
222  // Clear object
223  clear();
224 
225  // Continue only if there are standard keywords
226  if (hdu.has_card("CTYPE1") && hdu.has_card("CTYPE2") &&
227  hdu.has_card("CRVAL1") && hdu.has_card("CRVAL2") &&
228  hdu.has_card("CRPIX1") && hdu.has_card("CRPIX2")) {
229 
230  // Get standard keywords. Get bin size either from CDELT1/CDELT2 or
231  // CD1_1/CD2_2
232  std::string ctype1 = hdu.string("CTYPE1");
233  std::string ctype2 = hdu.string("CTYPE2");
234  double crval1 = hdu.real("CRVAL1");
235  double crval2 = hdu.real("CRVAL2");
236  double crpix1 = hdu.real("CRPIX1");
237  double crpix2 = hdu.real("CRPIX2");
238  double cdelt1 = (hdu.has_card("CDELT1")) ? hdu.real("CDELT1") : hdu.real("CD1_1");
239  double cdelt2 = (hdu.has_card("CDELT2")) ? hdu.real("CDELT2") : hdu.real("CD2_2");
240 
241  // Determine coordinate system
242  std::string coords;
243  std::string xcoord = ctype1.substr(0,4);
244  std::string ycoord = ctype2.substr(0,4);
245  if (xcoord == "RA--" && ycoord == "DEC-") {
246  m_lng = 0;
247  m_lat = 1;
248  coords = "EQU";
249  }
250  else if (xcoord == "DEC-" && ycoord == "RA--") {
251  m_lng = 1;
252  m_lat = 0;
253  coords = "EQU";
254  }
255  else if (xcoord == "GLON" && ycoord == "GLAT") {
256  m_lng = 0;
257  m_lat = 1;
258  coords = "GAL";
259  }
260  else if (xcoord == "GLAT" && ycoord == "GLON") {
261  m_lng = 1;
262  m_lat = 0;
263  coords = "GAL";
264  }
265  else if (xcoord == "ELON" && ycoord == "ELAT") {
266  m_lng = 0;
267  m_lat = 1;
268  coords = "ECL";
269  }
270  else if (xcoord == "ELAT" && ycoord == "ELON") {
271  m_lng = 1;
272  m_lat = 0;
273  coords = "ECL";
274  }
275  else if (xcoord == "HLON" && ycoord == "HLAT") {
276  m_lng = 0;
277  m_lat = 1;
278  coords = "HEL";
279  }
280  else if (xcoord == "HLAT" && ycoord == "HLON") {
281  m_lng = 1;
282  m_lat = 0;
283  coords = "HEL";
284  }
285  else if (xcoord == "SLON" && ycoord == "SLAT") {
286  m_lng = 0;
287  m_lat = 1;
288  coords = "SGL";
289  }
290  else if (xcoord == "SLAT" && ycoord == "SLON") {
291  m_lng = 1;
292  m_lat = 0;
293  coords = "SGL";
294  }
295  else {
297  }
298 
299  // Set standard parameters
300  set_members(coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2);
301 
302  // Optionally read "RADESYS" and "EQUINOX" keywords (needs to come
303  // before wcs_set() since the method will eventually complement
304  // missing keywords)
305  if (hdu.has_card("RADESYS")) {
306  m_radesys = hdu.string("RADESYS");
307  }
308  if (hdu.has_card("EQUINOX")) {
309  m_equinox = hdu.real("EQUINOX");
310  }
311 
312  // Setup WCS derived parameters
313  wcs_set();
314 
315  } // endif: there were standard keywords
316 
317  // Return
318  return;
319 }
320 
321 
322 /***********************************************************************//**
323  * @brief Write WCS definition into FITS HDU header
324  *
325  * @param[in] hdu FITS HDU.
326  *
327  * This method writes the World Coordinate System definition into the FITS
328  * HDU header.
329  *
330  * This method has been adapted from wcshdr.c::wcshdo().
331  ***************************************************************************/
332 void GWcs::write(GFitsHDU& hdu) const
333 {
334  // Continue only if there are axes
335  if (m_naxis > 0) {
336 
337  // Write reference pixel coordinates
338  for (int i = 0; i < m_naxis; ++i) {
339  std::string keyname = "CRPIX"+gammalib::str(i+1);
340  std::string comment = "Pixel coordinate of reference point (starting from 1)";
341  hdu.card(keyname, m_crpix.at(i), comment);
342  }
343 
344  //TODO: Write linear transformation matrix
345 
346  // Write coordinate increment at reference point
347  for (int i = 0; i < m_naxis; ++i) {
348  std::string keyname = "CDELT"+gammalib::str(i+1);
349  std::string comment;
350  if (m_cunit.at(i).length() > 0) {
351  comment += "["+gammalib::strip_whitespace(m_cunit.at(i))+"] ";
352  }
353  comment += "Coordinate increment at reference point";
354  hdu.card(keyname, m_cdelt.at(i), comment);
355  }
356 
357  // Write units of coordinate increment and reference value
358  for (int i = 0; i < m_naxis; ++i) {
359  if (m_cunit.at(i).length() > 0) {
360  std::string keyname = "CUNIT"+gammalib::str(i+1);
361  std::string comment = "Units of coordinate increment and value";
362  hdu.card(keyname, m_cunit.at(0), comment);
363  }
364  }
365 
366  // Write coordinate type
367  wcs_set_ctype();
368  for (int i = 0; i < m_naxis; ++i) {
369  if (i == m_lng) {
370  std::string keyname = "CTYPE"+gammalib::str(i+1);
371  hdu.card(keyname, m_ctype.at(i), m_ctype_c.at(i));
372  }
373  if (i == m_lat) {
374  std::string keyname = "CTYPE"+gammalib::str(i+1);
375  hdu.card(keyname, m_ctype.at(i), m_ctype_c.at(i));
376  }
377  if (i == m_spec) {
378  std::string keyname = "CTYPE"+gammalib::str(i+1);
379  hdu.card(keyname, m_ctype.at(i), m_ctype_c.at(i));
380  }
381  }
382 
383  // Write coordinate value at reference point
384  for (int i = 0; i < m_naxis; ++i) {
385  std::string keyname = "CRVAL"+gammalib::str(i+1);
386  std::string comment;
387  if (m_cunit.at(i).length() > 0) {
388  comment += "["+gammalib::strip_whitespace(m_cunit.at(i))+"] ";
389  }
390  comment += "Coordinate value at reference point";
391  hdu.card(keyname, m_crval.at(i), comment);
392  }
393 
394  //TODO: Parameter values
395  hdu.card("CROTA2", 0.0, "[deg] Rotation Angle"); // Old style, use PV instead
396  //hdu.card("PV2_1", 0.0, "Projection parameter 1");
397  //hdu.card("PV2_2", 0.0, "Projection parameter 2");
398 
399  // Celestial and spectral transformation parameters
400  if (!undefined(m_lonpole)) {
401  hdu.card("LONPOLE", m_lonpole, "[deg] Native longitude of celestial pole");
402  }
403  if (!undefined(m_latpole)) {
404  hdu.card("LATPOLE", m_latpole, "[deg] Native latitude of celestial pole");
405  }
406  if (!undefined(m_restfrq)) {
407  hdu.card("RESTFRQ", m_restfrq, "[Hz] Line rest frequency");
408  }
409  if (!undefined(m_restwav)) {
410  hdu.card("RESTWAV", m_restwav, "[Hz] Line rest wavelength");
411  }
412 
413  // Equatorial coordinate system type
414  if (m_radesys.length() > 0) {
415  hdu.card("RADESYS", m_radesys, "Equatorial coordinate system");
416  }
417 
418  // Equinox of equatorial coordinate system
419  if (!undefined(m_equinox)) {
420  hdu.card("EQUINOX", m_equinox, "[yr] Equinox of equatorial coordinates");
421  }
422 
423  //TODO: Reference frame of spectral coordinates
424 
425  //TODO: Reference frame of spectral observation
426 
427  //TODO: Observer's velocity towards source
428 
429  //TODO: Reference frame of source redshift
430 
431  //TODO: Redshift of the source
432 
433  //TODO: Observatory coordinates
434 
435  //TODO: MJD of observation
436 
437  //TODO: MJD mid-observation time
438 
439  //TODO: ISO-8601 date corresponding to MJD-OBS
440 
441  //TODO: ISO-8601 date corresponding to MJD-AVG
442 
443  } // endif: there were axes
444 
445  // Return
446  return;
447 }
448 
449 
450 /***********************************************************************//**
451  * @brief Returns solid angle of pixel in units of steradians
452  *
453  * @param[in] pixel Pixel index (x,y)
454  *
455  * Estimates solid angles of pixels using the Girard equation for excess
456  * area - see: http://mathworld.wolfram.com/SphericalPolygon.html
457  *
458  *
459  * Below, the definiton of the pixel cornes and sides are shown as used
460  * within the code.
461  *
462  * a12
463  * 1---------2
464  * |\ /|
465  * | \a13 / |
466  * | \ / |
467  * | \ / |
468  * a14| X |a23
469  * | / \ |
470  * | / \ |
471  * | /a24 \ |
472  * |/ \|
473  * 4---------3
474  * a34
475  *
476  ***************************************************************************/
477 double GWcs::solidangle(const GSkyPixel& pixel) const
478 {
479  // Get the sky directions of the 4 points
480  GSkyDir dir1 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()-0.5));
481  GSkyDir dir2 = pix2dir(GSkyPixel(pixel.x()+0.5, pixel.y()-0.5));
482  GSkyDir dir3 = pix2dir(GSkyPixel(pixel.x()+0.5, pixel.y()+0.5));
483  GSkyDir dir4 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()+0.5));
484 
485  // Compile option 1: Use triginometric function
486  #if G_SOLIDANGLE_OPTION == 1
487 
488  // Initialise solid angle
489  double solidangle = 0.0;
490 
491  // Compute angular distances between pixel corners
492  double a12 = dir1.dist(dir2);
493  double a14 = dir1.dist(dir4);
494  double a23 = dir2.dist(dir3);
495  double a34 = dir3.dist(dir4);
496 
497  // Special case: a12 or a14 is zero, then pixel is a triangle composed
498  // of [2,3,4]
499  if (a12 <= 0.0 || a14 <= 0.0) {
500 
501  // Compute diagonals
502  double a24 = dir2.dist(dir4);
503 
504  // Compute sines
505  double sin_a23 = std::sin(a23);
506  double sin_a24 = std::sin(a24);
507  double sin_a34 = std::sin(a34);
508 
509  // Compute cosines
510  double cos_a23 = std::cos(a23);
511  double cos_a24 = std::cos(a24);
512  double cos_a34 = std::cos(a34);
513 
514  // Compute angles
515  double angle2 = gammalib::acos((cos_a34-cos_a24*cos_a23)/(sin_a24*sin_a23));
516  double angle3 = gammalib::acos((cos_a24-cos_a34*cos_a23)/(sin_a34*sin_a23));
517  double angle4 = gammalib::acos((cos_a23-cos_a24*cos_a34)/(sin_a24*sin_a34));
518 
519  // Compute excess area to determine solid angle
520  solidangle = (angle2 + angle3 + angle4) - gammalib::pi;
521 
522  }
523 
524  // Special case: a23 or a 34 is zero, then pixel is a triangle composed
525  // of [1,2,4]
526  else if (a23 <= 0.0 || a34 <= 0.0) {
527 
528  // Compute diagonals
529  double a24 = dir2.dist(dir4);
530 
531  // Compute sines
532  double sin_a12 = std::sin(a12);
533  double sin_a14 = std::sin(a14);
534  double sin_a24 = std::sin(a24);
535 
536  // Compute cosines
537  double cos_a12 = std::cos(a12);
538  double cos_a14 = std::cos(a14);
539  double cos_a24 = std::cos(a24);
540 
541  // Compute angles
542  double angle1 = gammalib::acos((cos_a24-cos_a12*cos_a14)/(sin_a12*sin_a14));
543  double angle2 = gammalib::acos((cos_a14-cos_a12*cos_a24)/(sin_a12*sin_a24));
544  double angle4 = gammalib::acos((cos_a12-cos_a14*cos_a24)/(sin_a14*sin_a24));
545 
546  // Compute excess area to determine solid angle
547  solidangle = (angle1 + angle2 + angle4) - gammalib::pi;
548 
549  }
550 
551  // Otherwise we have a polygon
552  else {
553 
554  // Compute diagonals
555  double a13 = dir1.dist(dir3);
556  double a24 = dir2.dist(dir4);
557 
558  // Compute sines
559  double sin_a12 = std::sin(a12);
560  double sin_a14 = std::sin(a14);
561  double sin_a23 = std::sin(a23);
562  double sin_a34 = std::sin(a34);
563 
564  // Compute cosines
565  double cos_a12 = std::cos(a12);
566  double cos_a13 = std::cos(a13);
567  double cos_a14 = std::cos(a14);
568  double cos_a23 = std::cos(a23);
569  double cos_a24 = std::cos(a24);
570  double cos_a34 = std::cos(a34);
571 
572  // Compute angles
573  double angle1 = std::acos((cos_a13-cos_a34*cos_a14)/(sin_a34*sin_a14));
574  double angle2 = std::acos((cos_a24-cos_a23*cos_a34)/(sin_a23*sin_a34));
575  double angle3 = std::acos((cos_a13-cos_a12*cos_a23)/(sin_a12*sin_a23));
576  double angle4 = std::acos((cos_a24-cos_a14*cos_a12)/(sin_a14*sin_a12));
577 
578  // Use Girard equation for excess area to determine solid angle
579  solidangle = (angle1 + angle2 + angle3 + angle4) - gammalib::twopi;
580 
581  } // endif: we had a polynom
582 
583  // Compile option 2: Use Huilier's theorem
584  #elif G_SOLIDANGLE_OPTION == 2
585 
586  // Initialise solid angle
587  double solidangle = 0.0;
588 
589  // Compute angular distances between pixel corners
590  double a12 = dir1.dist(dir2);
591  double a14 = dir1.dist(dir4);
592  double a23 = dir2.dist(dir3);
593  double a24 = dir2.dist(dir4);
594  double a34 = dir3.dist(dir4);
595 
596  // Special case: a12 or a14 is zero, then pixel is a triangle composed
597  // of [2,3,4]
598  if (a12 <= 0.0 || a14 <= 0.0) {
599  double s = 0.5 * (a23 + a34 + a24);
600  solidangle = 4.0 * std::atan(std::sqrt(std::tan(0.5*s) *
601  std::tan(0.5*(s-a23)) *
602  std::tan(0.5*(s-a34)) *
603  std::tan(0.5*(s-a24))));
604  }
605 
606  // Special case: a23 or a 34 is zero, then pixel is a triangle composed
607  // of [1,2,4]
608  else if (a23 <= 0.0 || a34 <= 0.0) {
609  double s = 0.5 * (a12 + a24 + a14);
610  solidangle = 4.0 * std::atan(std::sqrt(std::tan(0.5*s) *
611  std::tan(0.5*(s-a12)) *
612  std::tan(0.5*(s-a24)) *
613  std::tan(0.5*(s-a14))));
614  }
615 
616  // Otherwise we have a polygon
617  else {
618 
619  // Triangle 1 [1,2,4]
620  double s1 = 0.5 * (a12 + a24 + a14);
621  double excess1 = std::atan(std::sqrt(std::tan(0.5*s1) *
622  std::tan(0.5*(s1-a12)) *
623  std::tan(0.5*(s1-a24)) *
624  std::tan(0.5*(s1-a14))));
625 
626  // Triangle 2 [2,3,4]
627  double s2 = 0.5 * (a23 + a34 + a24);
628  double excess2 = std::atan(std::sqrt(std::tan(0.5*s2) *
629  std::tan(0.5*(s2-a23)) *
630  std::tan(0.5*(s2-a34)) *
631  std::tan(0.5*(s2-a24))));
632 
633  // Determine solid angle
634  solidangle = 4.0 * (excess1 + excess2);
635 
636  } // endif: we had a polynom
637 
638  // Compile option 3: Use vectors
639  #else
640 
641  // Get vectors to pixel corners
642  GVector vec1 = dir1.celvector();
643  GVector vec2 = dir2.celvector();
644  GVector vec3 = dir3.celvector();
645  GVector vec4 = dir4.celvector();
646 
647  // Compute inner angles of pixel corners
648  double angle1 = gammalib::acos(cross(vec2, (cross(vec1, vec2))) *
649  cross(vec2, (cross(vec3, vec2))));
650  double angle2 = gammalib::acos(cross(vec3, (cross(vec2, vec3))) *
651  cross(vec3, (cross(vec4, vec3))));
652  double angle3 = gammalib::acos(cross(vec4, (cross(vec3, vec4))) *
653  cross(vec4, (cross(vec1, vec4))));
654  double angle4 = gammalib::acos(cross(vec1, (cross(vec4, vec1))) *
655  cross(vec1, (cross(vec2, vec1))));
656 
657  // Use Girard equation for excess area to determine solid angle
658  double solidangle = gammalib::deg2rad * gammalib::deg2rad +
659  (angle1 + angle2 + angle3 + angle4) -
661  #endif
662 
663  // Return solid angle
664  return solidangle;
665 }
666 
667 
668 /***********************************************************************//**
669  * @brief Returns sky direction of sky map pixel
670  *
671  * @param[in] pixel Sky map pixel.
672  * @return Sky direction.
673  *
674  * Returns the sky direction of a sky map @p pixel. Note that the sky
675  * map pixel values start from 0 while the WCS pixel reference starts from
676  * 1.
677  *
678  * A pre-computation cache is implemented so that successive calls with the
679  * same @p pixel value will returned the cached sky direction. The
680  * pre-computation cache is OMP thread safe.
681  ***************************************************************************/
682 GSkyDir GWcs::pix2dir(const GSkyPixel& pixel) const
683 {
684  // Initialise sky direction
685  GSkyDir dir;
686 
687  // If there is a pre-computation cache and the sky pixel is identical to
688  // the one used for the cache, then get the cached sky direction. We need
689  // to put this in a OMP critical zone so that no other thread will change
690  // the cache in the meantime
691  bool update_cache = true;
692  set_lock(2);
693  if (m_has_pix2dir_cache && (m_last_pix2dir_pix == pixel)) {
694  dir = m_last_pix2dir_dir;
695  update_cache = false;
696  }
697  unset_lock(2);
698 
699  // ... otherwise the computation of the sky direction is needed
700  if (update_cache) {
701 
702  // Allocate memory for transformation
703  double pixcrd[2];
704  double imgcrd[2];
705  double phi;
706  double theta;
707  double world[2];
708  int stat;
709 
710  // Set sky pixel. We have to add 1.0 here as the WCS pixel reference
711  // (CRPIX) starts from one while GSkyPixel starts from 0.
712  pixcrd[0] = pixel.x() + 1.0;
713  pixcrd[1] = pixel.y() + 1.0;
714 
715  // Trasform pixel-to-world coordinate
716  wcs_p2s(1, 2, pixcrd, imgcrd, &phi, &theta, world, &stat);
717 
718  // Set sky direction
719  if (m_coordsys == 0) {
720  dir.radec_deg(world[0], world[1]);
721  }
722  else {
723  dir.lb_deg(world[0], world[1]);
724  }
725 
726  // Store result in cache. We need to put this in a OMP critical zone
727  // to make sure that no other thread is fiddeling with the cache
728  set_lock(2);
729  {
730  m_has_pix2dir_cache = true;
731  m_last_pix2dir_pix = pixel;
732  m_last_pix2dir_dir = dir;
733  }
734  unset_lock(2);
735 
736  // Debug: Dump transformation steps
737  #if defined(G_XY2DIR_DEBUG)
738  std::cout << "xy2dir: pixel=" << pixel
739  << " (x,y)=(" << pixcrd[0] << "," << pixcrd[1] << ")"
740  << " (phi,theta)=(" << phi << "," << theta << ")"
741  << " (lng,lat)=(" << world[0] << "," << world[1] << ")"
742  << std::endl;
743  #endif
744 
745  } // endif: sky direction computation requested
746 
747  // Return
748  return dir;
749 }
750 
751 
752 /***********************************************************************//**
753  * @brief Returns sky map pixel of sky direction
754  *
755  * @param[in] dir Sky direction.
756  * @return Sky map pixel.
757  *
758  * Returns the sky map pixel for a given sky direction. Note that the sky
759  * map pixel values start from 0 while the WCS pixel reference starts from
760  * 1.
761  *
762  * A pre-computation cache is implemented so that successive calls with the
763  * same @p dir value will returned the cached sky pixel. The pre-computation
764  * cache is OMP thread safe.
765  ***************************************************************************/
766 GSkyPixel GWcs::dir2pix(const GSkyDir& dir) const
767 {
768  // Initialise sky pixel
769  GSkyPixel pixel;
770 
771  // If there is a pre-computation cache and the sky direction is identical
772  // to the one used for the cache, then get the cached sky pixel. We need
773  // to put this in a OMP critical zone so that no other thread will change
774  // the cache in the meantime
775  bool update_cache = true;
776  set_lock(1);
777  if (m_has_dir2pix_cache && (m_last_dir2pix_dir == dir)) {
778  pixel = m_last_dir2pix_pix;
779  update_cache = false;
780  }
781  unset_lock(1);
782 
783  // ... otherwise the computation of the sky pixel is needed
784  if (update_cache) {
785 
786  // Allocate memory for transformation
787  double pixcrd[2];
788  double imgcrd[2];
789  double phi;
790  double theta;
791  double world[2];
792  int stat;
793 
794  // Set world coordinate
795  if (m_coordsys == 0) {
796  world[0] = dir.ra_deg();
797  world[1] = dir.dec_deg();
798  }
799  else {
800  world[0] = dir.l_deg();
801  world[1] = dir.b_deg();
802  }
803 
804  // Transform world-to-pixel coordinate
805  wcs_s2p(1, 2, world, &phi, &theta, imgcrd, pixcrd, &stat);
806 
807  // Set sky pixel. We have to subtract 1 here as GSkyPixel starts from
808  // zero while the WCS reference (CRPIX) starts from one.
809  pixel.xy(pixcrd[0]-1.0, pixcrd[1]-1.0);
810 
811  // Store result in cache. We need to put this in a OMP critical zone
812  // to make sure that no other thread is fiddeling with the cache
813  set_lock(1);
814  {
815  m_has_dir2pix_cache = true;
816  m_last_dir2pix_dir = dir;
817  m_last_dir2pix_pix = pixel;
818  }
819  unset_lock(1);
820 
821  // Debug: Dump transformation steps
822  #if defined(G_DIR2XY_DEBUG)
823  std::cout << "dir2xy: dir=" << dir
824  << " (lng,lat)=(" << world[0] << "," << world[1] << ")"
825  << " (phi,theta)=(" << phi << "," << theta << ")"
826  << " (x,y)=" << pixel << std::endl;
827  #endif
828 
829  } // endif: sky pixel computation requested
830 
831  // Return sky pixel
832  return pixel;
833 }
834 
835 
836 /***********************************************************************//**
837  * @brief Set World Coordinate System parameters
838  *
839  * @param[in] coords Coordinate system.
840  * @param[in] crval1 X value of reference pixel.
841  * @param[in] crval2 Y value of reference pixel.
842  * @param[in] crpix1 X index of reference pixel (first pixel is 1).
843  * @param[in] crpix2 Y index of reference pixel (first pixel is 1).
844  * @param[in] cdelt1 Increment in x direction at reference pixel (deg).
845  * @param[in] cdelt2 Increment in y direction at reference pixel (deg).
846  *
847  * This method sets the WCS parameters.
848  ***************************************************************************/
849 void GWcs::set(const std::string& coords,
850  const double& crval1, const double& crval2,
851  const double& crpix1, const double& crpix2,
852  const double& cdelt1, const double& cdelt2)
853 {
854  // Clear any existing information
855  clear();
856 
857  // Set standard parameters
858  set_members(coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2);
859 
860  // Setup WCS derived parameters
861  wcs_set();
862 
863  // Return
864  return;
865 }
866 
867 
868 /***********************************************************************//**
869  * @brief Return value of reference pixel
870  *
871  * @param[in] inx Coordinate index.
872  * @return Value of reference pixel.
873  *
874  * @exception GException::out_of_range
875  * Index is out of valid range.
876  ***************************************************************************/
877 double GWcs::crval(const int& inx) const
878 {
879  // Compile option: raise an exception if index is out of range
880  #if defined(G_RANGE_CHECK)
881  if (inx < 0 || inx >= m_naxis) {
883  }
884  #endif
885 
886  // Return crval
887  return (m_crval[inx]);
888 }
889 
890 
891 /***********************************************************************//**
892  * @brief Return reference pixel
893  *
894  * @param[in] inx Coordinate index.
895  * @return Reference pixel.
896  *
897  * @exception GException::out_of_range
898  * Index is out of valid range.
899  ***************************************************************************/
900 double GWcs::crpix(const int& inx) const
901 {
902  // Compile option: raise an exception if index is out of range
903  #if defined(G_RANGE_CHECK)
904  if (inx < 0 || inx >= m_naxis) {
906  }
907  #endif
908 
909  // Return crpix
910  return (m_crpix[inx]);
911 }
912 
913 
914 /***********************************************************************//**
915  * @brief Return pixel size
916  *
917  * @param[in] inx Coordinate index.
918  * @return Pixel size.
919  *
920  * @exception GException::out_of_range
921  * Index is out of valid range.
922  ***************************************************************************/
923 double GWcs::cdelt(const int& inx) const
924 {
925  // Compile option: raise an exception if index is out of range
926  #if defined(G_RANGE_CHECK)
927  if (inx < 0 || inx >= m_naxis) {
929  }
930  #endif
931 
932  // Return cdelt
933  return (m_cdelt[inx]);
934 }
935 
936 
937 /*==========================================================================
938  = =
939  = Protected methods =
940  = =
941  ==========================================================================*/
942 
943 /***********************************************************************//**
944  * @brief Initialise class members
945  *
946  * Code adapted from wcs.c::wcsini(). In addition, the method sets up the
947  * World Coordinate System by calling wcs_set().
948  ***************************************************************************/
950 {
951  // Initialise World Coordinate System with 0 axes
952  wcs_ini(0);
953 
954  // Setup World Coordinate System
955  wcs_set();
956 
957  // Initialise pre-computation cache
958  m_has_pix2dir_cache = false;
959  m_has_dir2pix_cache = false;
964 
965  // Initialize the locks
966  init_lock(1);
967  init_lock(2);
968 
969  // Return
970  return;
971 }
972 
973 
974 /***********************************************************************//**
975  * @brief Copy class members
976  *
977  * @param[in] wcs GWcs instance from which members should be copied
978  ***************************************************************************/
979 void GWcs::copy_members(const GWcs& wcs)
980 {
981  // Copy WCS attributes
982  m_wcsset = wcs.m_wcsset;
983  m_naxis = wcs.m_naxis;
984  m_crval = wcs.m_crval;
985  m_cunit = wcs.m_cunit;
986  m_ctype = wcs.m_ctype;
987  m_ctype_c = wcs.m_ctype_c;
988  m_lonpole = wcs.m_lonpole;
989  m_latpole = wcs.m_latpole;
990  m_restfrq = wcs.m_restfrq;
991  m_restwav = wcs.m_restwav;
992  m_radesys = wcs.m_radesys;
993  m_equinox = wcs.m_equinox;
994  m_cd = wcs.m_cd;
995  m_crota = wcs.m_crota;
996  m_lng = wcs.m_lng;
997  m_lat = wcs.m_lat;
998  m_spec = wcs.m_spec;
999 
1000  // Copy linear transformation parameters
1001  m_linset = wcs.m_linset;
1002  m_unity = wcs.m_unity;
1003  m_crpix = wcs.m_crpix;
1004  m_pc = wcs.m_pc;
1005  m_cdelt = wcs.m_cdelt;
1006  m_piximg = wcs.m_piximg;
1007  m_imgpix = wcs.m_imgpix;
1008 
1009  // Copy celestial transformation parameters
1010  m_celset = wcs.m_celset;
1011  m_offset = wcs.m_offset;
1012  m_phi0 = wcs.m_phi0;
1013  m_theta0 = wcs.m_theta0;
1014  m_ref[0] = wcs.m_ref[0];
1015  m_ref[1] = wcs.m_ref[1];
1016  m_ref[2] = wcs.m_ref[2];
1017  m_ref[3] = wcs.m_ref[3];
1018  m_euler[0] = wcs.m_euler[0];
1019  m_euler[1] = wcs.m_euler[1];
1020  m_euler[2] = wcs.m_euler[2];
1021  m_euler[3] = wcs.m_euler[3];
1022  m_euler[4] = wcs.m_euler[4];
1023  m_latpreq = wcs.m_latpreq;
1024  m_isolat = wcs.m_isolat;
1025 
1026  // Copy projection parameters
1027  m_prjset = wcs.m_prjset;
1028  m_r0 = wcs.m_r0;
1029  m_bounds = wcs.m_bounds;
1030  m_x0 = wcs.m_x0;
1031  m_y0 = wcs.m_y0;
1032  m_w = wcs.m_w;
1033  for (int i = 0; i < PVN; ++i) {
1034  m_pv[i] = wcs.m_pv[i];
1035  }
1036 
1037  // Copy pre-computation cache
1044 
1045  // Return
1046  return;
1047 }
1048 
1049 
1050 /***********************************************************************//**
1051  * @brief Delete class members
1052  ***************************************************************************/
1054 {
1055  // Return
1056  return;
1057 }
1058 
1059 
1060 /***********************************************************************//**
1061  * @brief Set World Coordinate System parameters
1062  *
1063  * @param[in] coords Coordinate system.
1064  * @param[in] crval1 X value of reference pixel.
1065  * @param[in] crval2 Y value of reference pixel.
1066  * @param[in] crpix1 X index of reference pixel (first pixel is 1).
1067  * @param[in] crpix2 Y index of reference pixel (first pixel is 1).
1068  * @param[in] cdelt1 Increment in x direction at reference pixel (deg).
1069  * @param[in] cdelt2 Increment in y direction at reference pixel (deg).
1070  *
1071  * This method sets the WCS parameters. It does not call wcs_set(), however,
1072  * as wcs_set() may depend on the availability of derived class methods which
1073  * can not be used in a base class constructor.
1074  *
1075  * @todo Implement parameter validity check
1076  ***************************************************************************/
1077 void GWcs::set_members(const std::string& coords,
1078  const double& crval1, const double& crval2,
1079  const double& crpix1, const double& crpix2,
1080  const double& cdelt1, const double& cdelt2)
1081 
1082 {
1083  //TODO: Check parameters
1084 
1085  // Initialise WCS
1086  wcs_ini(2);
1087 
1088  // Set coordinate system
1089  coordsys(coords);
1090 
1091  // Set World Coordinate parameters
1092  m_lng = 0;
1093  m_lat = 1;
1094  m_crpix[0] = crpix1;
1095  m_crpix[1] = crpix2;
1096  m_cdelt[0] = cdelt1;
1097  m_cdelt[1] = cdelt2;
1098  m_crval[0] = crval1;
1099  m_crval[1] = crval2;
1100 
1101  // Return
1102  return;
1103 }
1104 
1105 
1106 /***********************************************************************//**
1107  * @brief Compares sky projection
1108  *
1109  * @param[in] proj Sky projection.
1110  * @return True if @p proj is identical to the actual object.
1111  *
1112  * This method is a helper for the World Coordinate Comparison friends. It
1113  * does not compare derived quantities as those may not have been initialised
1114  * in both objects.
1115  ***************************************************************************/
1116 bool GWcs::compare(const GSkyProjection& proj) const
1117 {
1118  // Initialise result
1119  bool result = false;
1120 
1121  // Continue only we compare to a GWcs object
1122  const GWcs* ptr = dynamic_cast<const GWcs*>(&proj);
1123  if (ptr != NULL) {
1124 
1125  // Perform comparion of standard (non derived) parameters
1126  result = ((code() == ptr->code()) &&
1127  (m_coordsys == ptr->m_coordsys) &&
1128  (m_lng == ptr->m_lng) &&
1129  (m_lat == ptr->m_lat) &&
1130  (m_crval == ptr->m_crval) &&
1131  (m_crpix == ptr->m_crpix) &&
1132  (m_cdelt == ptr->m_cdelt));
1133 
1134  }
1135 
1136  // Return result
1137  return result;
1138 }
1139 
1140 
1141 /*==========================================================================
1142  = =
1143  = World Coordinate Methods adapted from wcslib =
1144  = =
1145  ==========================================================================*/
1146 
1147 /***********************************************************************//**
1148  * @brief Initialise World Coordinate System
1149  *
1150  * @param[in] naxis Number of axes
1151  *
1152  * This method initialises the World Coordinate System information. The
1153  * method has been adapted from the wcslib function wcs.c::wcsini. In
1154  * contrast to the wcslib function, however, this method accepts naxis=0.
1155  * In this case, all vectors are cleared.
1156  ***************************************************************************/
1157 void GWcs::wcs_ini(int naxis)
1158 {
1159  // Signal that WCS information has not been set so far. This will be done
1160  // in wcs_set() upon request
1161  m_wcsset = false;
1162 
1163  // Clear vectors
1164  m_crpix.clear();
1165  m_pc.clear();
1166  m_cdelt.clear();
1167  m_crval.clear();
1168  m_cunit.clear();
1169  m_ctype.clear();
1170  m_ctype_c.clear();
1171  m_cd.clear();
1172  m_crota.clear();
1173 
1174  // Set number of axes
1175  m_naxis = naxis;
1176 
1177  // Set defaults for the linear transformation (sets m_Crpix, m_pc, m_Cdelt)
1178  lin_ini(naxis);
1179 
1180  // Set default axes parameters
1181  if (naxis > 0) {
1182  for (int i = 0; i < naxis; ++i) {
1183  m_crval.push_back(0.0);
1184  m_cunit.push_back("");
1185  m_ctype.push_back("");
1186  m_ctype_c.push_back("");
1187  }
1188  }
1189 
1190  // Set defaults for the celestial transformation parameters
1191  m_lonpole = UNDEFINED;
1192  m_latpole = +90.0;
1193 
1194  // Set defaults for the spectral transformation parameters
1195  m_restfrq = UNDEFINED; // 0.0 in wcslib
1196  m_restwav = UNDEFINED; // 0.0 in wcslib
1197 
1198  //TODO: Default parameter values
1199 
1200  // Defaults for alternate linear transformations
1201  if (naxis > 0) {
1202  for (int i = 0; i < naxis; ++i) {
1203  for (int j = 0; j < naxis; ++j) {
1204  m_cd.push_back(0.0);
1205  }
1206  m_crota.push_back(0.0);
1207  }
1208  }
1209 
1210  //TODO: Defaults for auxiliary coordinate system information
1211  m_radesys.clear();
1212  m_equinox = UNDEFINED;
1213 
1214  // Reset derived values
1215  m_lng = -1;
1216  m_lat = -1;
1217  m_spec = -1;
1218 
1219  // Initialise celestial transformation parameters
1220  cel_ini();
1221 
1222  // Initialise spectral transformation parameters
1223  spc_ini();
1224 
1225  // Return
1226  return;
1227 }
1228 
1229 
1230 /***********************************************************************//**
1231  * @brief Setup of World Coordinate System
1232  *
1233  * This method sets up the World Coordinate System information. In particular
1234  * it:
1235  * - initialises the celestial parameters (call to cel_ini)
1236  * - sets the m_ref information (reference value and native pole)
1237  * - sets the celestial parameters (call to cel_set)
1238  * - sets the linear transformation (call to lin_set)
1239  * The method has been adapted from the wcslib function wcs.c::wcsset.
1240  *
1241  * @todo Determine axis types from CTYPEia
1242  * @todo Convert to canonical units
1243  * @todo Do we have PVi_ma keyvalues?
1244  * @todo Do simple alias translations
1245  * @todo Update PVi_ma keyvalues
1246  * @todo Non-linear spectral axis present?
1247  * @todo Tabular axes present?
1248  ***************************************************************************/
1249 void GWcs::wcs_set(void) const
1250 {
1251  //TODO: Determine axis types from CTYPEia
1252 
1253  //TODO: Convert to canonical units
1254 
1255  // Non-linear celestial axes present?
1256  if (m_lng >= 0) {
1257 
1258  // Initialise celestial parameters
1259  cel_ini();
1260 
1261  // Set CRVALia, LONPOLEa, and LATPOLEa keyvalues
1262  m_ref[0] = m_crval[m_lng]; // Longitude reference value
1263  m_ref[1] = m_crval[m_lat]; // Latitude reference value
1264  m_ref[2] = m_lonpole; // LONPOLE
1265  m_ref[3] = m_latpole; // LATPOLE
1266 
1267  //TODO: Do we have PVi_ma keyvalues?
1268 
1269  // Do simple alias translations
1270  if (code() == "GLS") {
1271  m_offset = true;
1272  m_phi0 = 0.0;
1273  m_theta0 = m_crval[m_lat];
1274  }
1275 
1276  // Initialize the celestial transformation routines
1277  m_r0 = 0.0; // Forces initialisation
1278  cel_set();
1279 
1280  // Update LONPOLE and LATPOLE
1281  m_lonpole = m_ref[2];
1282  m_latpole = m_ref[3];
1283 
1284  //TODO: Update PVi_ma keyvalues
1285 
1286  } // endif: non-linear celestial axes were present
1287 
1288  //TODO: Non-linear spectral axis present?
1289 
1290  //TODO: Tabular axes present?
1291 
1292  // Initialize the linear transformation
1293  lin_set();
1294 
1295  // Set defaults for radesys and equinox
1296  wcs_set_radesys();
1297 
1298  // Signal that WCS is set
1299  m_wcsset = true;
1300 
1301  // Return
1302  return;
1303 }
1304 
1305 
1306 /***********************************************************************//**
1307  * @brief Set radesys and equinox members
1308  *
1309  * Set default values for radesys and equinox for equatorial, ecliptic or
1310  * helioecliptic coordinates. For other coordinates the radesys and equinox
1311  * members are set to empty/undefined values.
1312  *
1313  * This method has been inspired by code from wcs.c::wcsset.
1314  ***************************************************************************/
1315 void GWcs::wcs_set_radesys(void) const
1316 {
1317  // Set defaults for radesys and equinox for equatorial, ecliptic or
1318  // helioecliptic coordinates
1319  if ((m_coordsys == 0) || (m_coordsys == 2) || (m_coordsys == 3)) {
1320 
1321  // If radesys is undefined then set it dependend on the equinox.
1322  // If equinox is undefined then use "ICRS", otherwise use "FK4" if
1323  // equinox is smaller than 1984 and "FK5" for later equinoxes
1324  if (m_radesys.empty()) {
1325  if (undefined(m_equinox)) {
1326  m_radesys = "ICRS";
1327  }
1328  else if (m_equinox < 1984.0) {
1329  m_radesys = "FK4";
1330  }
1331  else {
1332  m_radesys = "FK5";
1333  }
1334  }
1335 
1336  // ... otherwise, if radesys is "ICRS" or "GAPPT" then set the equinox
1337  // to undefined since the equinox is not applicable for these coordinate
1338  // systems
1339  else if ((m_radesys == "ICRS") || (m_radesys == "GAPPT")) {
1340  m_equinox = UNDEFINED;
1341  }
1342 
1343  // ... otherwise, if equinox is undefined and radesys is "FK5" then
1344  // set equinox to 2000.0 and if RADESYS is "FK4" or "FK4-NO-E" then
1345  // set equinox to 1950.0
1346  else if (undefined(m_equinox)) {
1347  if (m_radesys == "FK5") {
1348  m_equinox = 2000.0;
1349  }
1350  else if ((m_radesys == "FK4") || (m_radesys == "FK4-NO-E")) {
1351  m_equinox = 1950.0;
1352  }
1353  }
1354 
1355  } // endif: coordinate system requires RADESYS
1356 
1357  // ... otherwise reset radesys and set equinox to undefined
1358  else {
1359  m_radesys.clear();
1360  m_equinox = UNDEFINED;
1361  }
1362 
1363  // Return
1364  return;
1365 }
1366 
1367 
1368 /***********************************************************************//**
1369  * @brief Set CTYPEa keywords
1370  *
1371  * @exception GException::wcs_invalid
1372  * WCS projection not valid.
1373  * @exception GException::wcs_bad_coords
1374  * Coordinate system is not of valid type.
1375  *
1376  * This method has been inspired by code from wcshdr.c::wcshdo.
1377  ***************************************************************************/
1378 void GWcs::wcs_set_ctype(void) const
1379 {
1380  // Check on correct type length
1381  if (code().length() != 3) {
1383  "3-character type required.");
1384  }
1385 
1386  // Set longitude keyword
1387  if (m_lng >= 0) {
1388 
1389  // Set coordinate system
1390  if (m_coordsys == 0) {
1391  m_ctype.at(m_lng) = "RA---" + code();
1392  m_ctype_c.at(m_lng) = "Right ascension, ";
1393  }
1394  else if (m_coordsys == 1) {
1395  m_ctype.at(m_lng) = "GLON-" + code();
1396  m_ctype_c.at(m_lng) = "Galactic longitude, ";
1397  }
1398  else if (m_coordsys == 2) {
1399  m_ctype.at(m_lng) = "ELON-" + code();
1400  m_ctype_c.at(m_lng) = "Ecliptic longitude, ";
1401  }
1402  else if (m_coordsys == 3) {
1403  m_ctype.at(m_lng) = "HLON-" + code();
1404  m_ctype_c.at(m_lng) = "Helioecliptic longitude, ";
1405  }
1406  else if (m_coordsys == 4) {
1407  m_ctype.at(m_lng) = "SLON-" + code();
1408  m_ctype_c.at(m_lng) = "Supergalactic longitude, ";
1409  }
1410  else {
1412  }
1413 
1414  // Add projection name to comment
1415  m_ctype_c.at(m_lng).append(name());
1416  m_ctype_c.at(m_lng).append(" projection");
1417  }
1418 
1419  // Set latitude keyword
1420  if (m_lat >= 0) {
1421 
1422  // Set coordinate system
1423  if (m_coordsys == 0) {
1424  m_ctype.at(m_lat) = "DEC--" + code();
1425  m_ctype_c.at(m_lat) = "Declination, ";
1426  }
1427  else if (m_coordsys == 1) {
1428  m_ctype.at(m_lat) = "GLAT-" + code();
1429  m_ctype_c.at(m_lat) = "Galactic latitude, ";
1430  }
1431  else if (m_coordsys == 2) {
1432  m_ctype.at(m_lat) = "ELAT-" + code();
1433  m_ctype_c.at(m_lat) = "Ecliptic latitude, ";
1434  }
1435  else if (m_coordsys == 3) {
1436  m_ctype.at(m_lat) = "HLAT-" + code();
1437  m_ctype_c.at(m_lat) = "Helioecliptic latitude, ";
1438  }
1439  else if (m_coordsys == 4) {
1440  m_ctype.at(m_lat) = "SLAT-" + code();
1441  m_ctype_c.at(m_lat) = "Supergalactic latitude, ";
1442  }
1443  else {
1445  }
1446 
1447  // Add projection name to comment
1448  m_ctype_c.at(m_lat).append(name());
1449  m_ctype_c.at(m_lat).append(" projection");
1450  }
1451 
1452  //TODO: Set spectral keyword
1453  if (m_spec >= 0) {
1454  m_ctype.at(m_spec) = "NONE";
1455  m_ctype_c.at(m_spec) = "Not yet implemented";
1456  }
1457 
1458  // Return
1459  return;
1460 }
1461 
1462 
1463 /***********************************************************************//**
1464  * @brief Pixel-to-world transformation
1465  *
1466  * @param[in] ncoord Number of coordinates.
1467  * @param[in] nelem Vector length of each coordinate (>=m_naxis).
1468  * @param[in] pixcrd Array [ncoord][nelem] of pixel coordinates.
1469  * @param[out] imgcrd Array [ncoord][nelem] of intermediate world coordinates.
1470  * @param[out] phi Array [ncoord] of longitudes in native coordinate system.
1471  * @param[out] theta Array [ncoord] of latitudes in native coordinate system.
1472  * @param[out] world Array [ncoord][nelem] of world coordinates.
1473  * @param[out] stat Array [ncoord] of pixel coordinates validity.
1474  *
1475  * @exception GException::wcs_invalid_parameter
1476  * Invalid input parameters provided
1477  *
1478  * This method transforms pixel coordinates to world coordinates. The method
1479  * has been adapted from wcs.c::wcsp2s(). Note that this method is extremely
1480  * simplified with respect to wcslib, but it does the job for now.
1481  *
1482  * For celestial axes, imgcrd[][wcs.lng] and imgcrd[][wcs.lat] are the
1483  * projected x-, and y-coordinates in pseudo "degrees" and world[][wcs.lng]
1484  * and world[][wcs.lat] are the celestial longitude and latitude [deg]
1485  * For spectral axes, imgcrd[][wcs.spec] is the intermediate spectral
1486  * coordinate, in SI units and world[][wcs.spec] is ...
1487  *
1488  * @todo Check for constant x and/or y to speed-up computations
1489  * @todo Zero the unused world coordinate elements
1490  ***************************************************************************/
1491 void GWcs::wcs_p2s(int ncoord, int nelem, const double* pixcrd, double* imgcrd,
1492  double* phi, double* theta, double* world, int* stat) const
1493 {
1494  // Initialize if required
1495  if (!m_wcsset) {
1496  wcs_set();
1497  }
1498 
1499  // Sanity check
1500  if (ncoord < 1 || (ncoord > 1 && nelem < m_naxis)) {
1501  std::string message;
1502  if (ncoord < 1) {
1503  message = "ncoord="+gammalib::str(ncoord)+", >0 required.";
1504  }
1505  else {
1506  message = "nelem="+gammalib::str(nelem)+", >="+gammalib::str(m_naxis)+" required.";
1507  }
1509  }
1510 
1511  // Apply pixel-to-world linear transformation
1512  lin_p2x(ncoord, nelem, pixcrd, imgcrd);
1513 
1514  //TODO: Check for constant x and/or y
1515  int nx = ncoord;
1516  int ny = 0;
1517 
1518  // Transform projection plane coordinates to celestial coordinates
1519  cel_x2s(nx, ny, nelem, nelem, imgcrd+m_lng, imgcrd+m_lat, phi, theta,
1520  world+m_lng, world+m_lat, stat);
1521 
1522  //TODO: Zero the unused world coordinate elements
1523 
1524  // Return
1525  return;
1526 }
1527 
1528 
1529 /***********************************************************************//**
1530  * @brief World-to-pixel transformation
1531  *
1532  * @param[in] ncoord Number of coordinates.
1533  * @param[in] nelem Vector length of each coordinate (>=m_naxis).
1534  * @param[in] world Array [ncoord][nelem] of world coordinates.
1535  * @param[out] phi Array [ncoord] of longitudes in native coordinate system.
1536  * @param[out] theta Array [ncoord] of latitudes in native coordinate system.
1537  * @param[out] imgcrd Array [ncoord][nelem] of intermediate world coordinates.
1538  * @param[out] pixcrd Array [ncoord][nelem] of pixel coordinates.
1539  * @param[out] stat Array [ncoord] of pixel coordinates validity.
1540  *
1541  * @exception GException::wcs_invalid_parameter
1542  * Invalid input parameters provided
1543  *
1544  * This method transforms world coordinates to pixel coordinates. The method
1545  * has been adapted from wcs.c::wcss2p(). Note that this method is extremely
1546  * simplified with respect to wcslib, but it does the job for now.
1547  *
1548  * For celestial axes, imgcrd[][wcs.lng] and imgcrd[][wcs.lat] are the
1549  * projected x-, and y-coordinates in pseudo "degrees" and world[][wcs.lng]
1550  * and world[][wcs.lat] are the celestial longitude and latitude [deg]
1551  * For spectral axes, imgcrd[][wcs.spec] is the intermediate spectral
1552  * coordinate, in SI units and world[][wcs.spec] is ...
1553  *
1554  * @todo Check for constant x and/or y to speed-up computations
1555  * @todo Zero the unused world coordinate elements
1556  ***************************************************************************/
1557 void GWcs::wcs_s2p(int ncoord, int nelem, const double* world,
1558  double* phi, double* theta, double* imgcrd,
1559  double* pixcrd, int* stat) const
1560 {
1561  // Initialize if required
1562  if (!m_wcsset) {
1563  wcs_set();
1564  }
1565 
1566  // Sanity check
1567  if (ncoord < 1 || (ncoord > 1 && nelem < m_naxis)) {
1568  std::string message;
1569  if (ncoord < 1) {
1570  message = "ncoord="+gammalib::str(ncoord)+", >0 required.";
1571  }
1572  else {
1573  message = "nelem="+gammalib::str(nelem)+", >="+gammalib::str(m_naxis)+" required.";
1574  }
1576  }
1577 
1578  //TODO: Check for constant x and/or y
1579  int nlng = ncoord;
1580  int nlat = 0;
1581 
1582  // Transform celestial coordinates to projection plane coordinates
1583  cel_s2x(nlng, nlat, nelem, nelem, world+m_lng, world+m_lat,
1584  phi, theta, imgcrd+m_lng, imgcrd+m_lat, stat);
1585 
1586  //TODO: Zero the unused world coordinate elements
1587 
1588  // Apply world-to-pixel linear transformation
1589  lin_x2p(ncoord, nelem, imgcrd, pixcrd);
1590 
1591  // Return
1592  return;
1593 }
1594 
1595 
1596 /***********************************************************************//**
1597  * @brief Print WCS information
1598  *
1599  * @param[in] chatter Chattiness.
1600  * @return String containing WCS information.
1601  ***************************************************************************/
1602 std::string GWcs::wcs_print(const GChatter& chatter) const
1603 {
1604  // Initialise result string
1605  std::string result;
1606 
1607  // Continue only if chatter is not silent
1608  if (chatter != SILENT) {
1609 
1610  // Append World Coordinate parameters
1611  result.append("\n"+gammalib::parformat("Number of axes"));
1612  result.append(gammalib::str(m_naxis));
1613  if (chatter >= EXPLICIT) {
1614  result.append("\n"+gammalib::parformat("Longitude axis"));
1615  result.append(gammalib::str(m_lng));
1616  result.append("\n"+gammalib::parformat("Latitude axis"));
1617  result.append(gammalib::str(m_lat));
1618  result.append("\n"+gammalib::parformat("Spectral axis"));
1619  result.append(gammalib::str(m_spec));
1620  }
1621 
1622  // Append coordinate system
1623  result.append("\n"+gammalib::parformat("Coodinate system")+coordsys());
1624 
1625  // Append projection parameters
1626  result.append("\n"+gammalib::parformat("Projection code")+code());
1627  result.append("\n"+gammalib::parformat("Projection name")+name());
1628  if (chatter >= EXPLICIT) {
1629  result.append("\n"+gammalib::parformat("Radius of the gen. sphere"));
1630  result.append(gammalib::str(m_r0)+" deg");
1631  }
1632 
1633  // Append coordinates
1634  result.append("\n"+gammalib::parformat("Reference coordinate")+"(");
1635  for (int i = 0; i < m_crval.size(); ++i) {
1636  if (i > 0) {
1637  result.append(", ");
1638  }
1639  result.append(gammalib::str(m_crval[i]));
1640  if (m_cunit[i].length() > 0) {
1641  result.append(" "+m_cunit[i]);
1642  }
1643  }
1644  result.append(")");
1645  result.append("\n"+gammalib::parformat("Reference pixel")+"(");
1646  for (int i = 0; i < m_crpix.size(); ++i) {
1647  if (i > 0) {
1648  result.append(", ");
1649  }
1650  result.append(gammalib::str(m_crpix[i]));
1651  }
1652  result.append(")");
1653  result.append("\n"+gammalib::parformat("Increment at reference")+"(");
1654  for (int i = 0; i < m_cdelt.size(); ++i) {
1655  if (i > 0) {
1656  result.append(", ");
1657  }
1658  result.append(gammalib::str(m_cdelt[i]));
1659  if (m_cunit[i].length() > 0) {
1660  result.append(" "+m_cunit[i]);
1661  }
1662  }
1663  result.append(")");
1664 
1665  // Append origin
1666  result.append("\n"+gammalib::parformat("(Phi_0, Theta_0)")+"(");
1667  result.append(wcs_print_value(m_phi0)+", ");
1668  result.append(wcs_print_value(m_theta0)+") deg");
1669 
1670  // Append native pole
1671  result.append("\n"+gammalib::parformat("(Phi_p, Theta_p)")+"(");
1672  result.append(wcs_print_value(m_lonpole)+", ");
1673  result.append(wcs_print_value(m_latpole)+") deg");
1674 
1675  // Append details
1676  if (chatter >= EXPLICIT) {
1677 
1678  // Append LATPOLEa keyword usage
1679  result.append("\n"+gammalib::parformat("LATPOLE keyword usage"));
1680  switch (m_latpreq) {
1681  case 0:
1682  result.append("Not used. Theta_p determined uniquely by"
1683  " CRVALia and LONPOLEa keywords.");
1684  break;
1685  case 1:
1686  result.append("Required to select between two valid solutions"
1687  " of Theta_p.");
1688  break;
1689  case 2:
1690  result.append("Theta_p was specified solely by LATPOLE.");
1691  break;
1692  default:
1693  result.append("UNDEFINED");
1694  break;
1695  }
1696 
1697  // Append celestial transformation parameters
1698  result.append("\n"+gammalib::parformat("Reference vector (m_ref)")+"(");
1699  for (int k = 0; k < 4; ++k) {
1700  if (k > 0) {
1701  result.append(", ");
1702  }
1703  result.append(wcs_print_value(m_ref[k]));
1704  }
1705  result.append(") deg");
1706 
1707  // Append Euler angles
1708  result.append("\n"+gammalib::parformat("Euler angles")+"(");
1709  for (int k = 0; k < 5; ++k) {
1710  if (k > 0) {
1711  result.append(", ");
1712  }
1713  result.append(gammalib::str(m_euler[k]));
1714  if (k < 3) {
1715  result.append(" deg");
1716  }
1717  }
1718  result.append(")");
1719 
1720  // Append latitude preservement flag
1721  result.append("\n"+gammalib::parformat("Latitude preserved"));
1722  if (m_isolat) {
1723  result.append("True");
1724  }
1725  else {
1726  result.append("False");
1727  }
1728 
1729  // Append linear transformation parameters
1730  result.append("\n"+gammalib::parformat("Unity PC matrix"));
1731  if (m_unity) {
1732  result.append("True");
1733  }
1734  else {
1735  result.append("False");
1736  }
1737  result.append("\n"+gammalib::parformat("Pixel-to-image trafo")+"(");
1738  for (int k = 0; k < m_piximg.size(); ++k) {
1739  if (k > 0) {
1740  result.append(", ");
1741  }
1742  result.append(gammalib::str(m_piximg[k]));
1743  }
1744  result.append(")");
1745  result.append("\n"+gammalib::parformat("Image-to-pixel trafo")+"(");
1746  for (int k = 0; k < m_imgpix.size(); ++k) {
1747  if (k > 0) {
1748  result.append(", ");
1749  }
1750  result.append(gammalib::str(m_imgpix[k]));
1751  }
1752  result.append(")");
1753 
1754  // Append boundary checking information
1755  result.append("\n"+gammalib::parformat("Strict bounds checking"));
1756  if (m_bounds) {
1757  result.append("True");
1758  }
1759  else {
1760  result.append("False");
1761  }
1762 
1763  // Append fiducial offset information
1764  result.append("\n"+gammalib::parformat("Use fiducial offset"));
1765  if (m_offset) {
1766  result.append("True");
1767  }
1768  else {
1769  result.append("False");
1770  }
1771  result.append("\n"+gammalib::parformat("Fiducial offset"));
1772  result.append("("+gammalib::str(m_x0)+", "+gammalib::str(m_y0)+")");
1773 
1774  // Append spectral transformation parameters
1775  result.append("\n"+gammalib::parformat("Rest frequency"));
1776  result.append(wcs_print_value(m_restfrq));
1777  result.append("\n"+gammalib::parformat("Rest wavelength"));
1778  result.append(wcs_print_value(m_restwav));
1779 
1780  } // endif: print details
1781 
1782  } // endif: chatter was not silent
1783 
1784  // Return
1785  return result;
1786 }
1787 
1788 
1789 /***********************************************************************//**
1790  * @brief Helper function for value printing
1791  *
1792  * @param[in] value Double precision value
1793  *
1794  * This helper function either prints a double precision value or the
1795  * word 'UNDEFINED', depending on whether the value is defined or not.
1796  ***************************************************************************/
1797 std::string GWcs::wcs_print_value(const double& value) const
1798 {
1799  // Initialise result
1800  std::string result;
1801 
1802  // Depend value dependent string
1803  if (undefined(value)) {
1804  result.append("UNDEFINED");
1805  }
1806  else {
1807  result.append(gammalib::str(value));
1808  }
1809 
1810  // Return result
1811  return result;
1812 }
1813 
1814 
1815 /*==========================================================================
1816  = =
1817  = Celestial transformation methods adapted from wcslib =
1818  = =
1819  ==========================================================================*/
1820 
1821 /***********************************************************************//**
1822  * @brief Initialise celestial transformation parameters
1823  *
1824  * Code adapted from cel.c::celini().
1825  ***************************************************************************/
1826 void GWcs::cel_ini(void) const
1827 {
1828  // Initialise parameters
1829  m_celset = false;
1830  m_offset = false;
1831  m_phi0 = UNDEFINED;
1832  m_theta0 = UNDEFINED;
1833  m_ref[0] = 0.0;
1834  m_ref[1] = 0.0;
1835  m_ref[2] = UNDEFINED;
1836  m_ref[3] = +90.0;
1837  m_latpreq = -1;
1838  m_isolat = false;
1839 
1840  // Clear Euler angles
1841  for (int k = 0; k < 5; ++k) {
1842  m_euler[k] = 0.0;
1843  }
1844 
1845  // Initialise projection parameters
1846  prj_ini();
1847 
1848  // Return
1849  return;
1850 }
1851 
1852 
1853 /***********************************************************************//**
1854  * @brief Setup of celestial transformation
1855  *
1856  * @exception GException::wcs_invalid_parameter
1857  * Invalid World Coordinate System parameter encountered.
1858  *
1859  * This method sets up the celestial transformation information. The method
1860  * has been adapted from the wcslib function cel.c::celset. It:
1861  * - sets the projection (using prj_set)
1862  * - computes the native longitude and latitude and stores the results
1863  * in m_ref[2] and m_ref[3]
1864  * - computes the Euler angles for the celestial transformation and stores
1865  * the result in m_euler
1866  * - it sets the m_latpreq (0=, 1=two solutions)
1867  * - it sets the m_isolat flag
1868  *
1869  * The m_latpreq member is for informational purposes and indicates how the
1870  * LATPOLEa keyword was used
1871  * - 0: Not required, theta_p (== delta_p) was determined uniquely by the
1872  * CRVALia and LONPOLEa keywords.
1873  * - 1: Required to select between two valid solutions of theta_p.
1874  * - 2: theta_p was specified solely by LATPOLEa.
1875  *
1876  * The m_isolat flag is true if the spherical rotation preserves the
1877  * magnitude of the latitude, which occurs if the axes of the native and
1878  * celestial coordinates are coincident. It signals an opportunity to cache
1879  * intermediate calculations common to all elements in a vector computation.
1880  ***************************************************************************/
1881 void GWcs::cel_set(void) const
1882 {
1883  // Set tolerance
1884  const double tol = 1.0e-10;
1885 
1886  // If no fiducial offsets are requested then make sure that (phi0,theta0)
1887  // are undefined
1888  if (!m_offset) {
1889  m_phi0 = UNDEFINED;
1890  m_theta0 = UNDEFINED;
1891  }
1892 
1893  // Setup projection. This method will set (phi0,theta0) if they were
1894  // undefined.
1895  prj_set();
1896 
1897  // Initialise celestial parameters
1898  double lng0 = m_ref[0];
1899  double lat0 = m_ref[1];
1900  double phip = m_ref[2];
1901  double latp = m_ref[3];
1902  double lngp = 0.0;
1903 
1904  // Set default for native longitude of the celestial pole?
1905  if (undefined(phip) || phip == 999.0) {
1906  phip = (lat0 < m_theta0) ? 180.0 : 0.0;
1907  phip += m_phi0;
1908  if (phip < -180.0) {
1909  phip += 360.0;
1910  }
1911  else if (phip > 180.0) {
1912  phip -= 360.0;
1913  }
1914  m_ref[2] = phip;
1915  }
1916 
1917  // Initialise
1918  m_latpreq = 0;
1919 
1920  // Compute celestial coordinates of the native pole
1921  // Fiducial point at the native pole?
1922  if (m_theta0 == 90.0) {
1923  lngp = lng0;
1924  latp = lat0;
1925  }
1926 
1927  // ... otherwise fiducial point is away from the native pole
1928  else {
1929 
1930  // Compute sine and cosine of lat0 and theta0
1931  double slat0;
1932  double clat0;
1933  double sthe0;
1934  double cthe0;
1935  gammalib::sincosd(lat0, &slat0, &clat0);
1936  gammalib::sincosd(m_theta0, &sthe0, &cthe0);
1937 
1938  // ...
1939  double sphip;
1940  double cphip;
1941  double u;
1942  double v;
1943  if (phip == m_phi0) {
1944  sphip = 0.0;
1945  cphip = 1.0;
1946  u = m_theta0;
1947  v = 90.0 - lat0;
1948  }
1949 
1950  // ...
1951  else {
1952 
1953  // Compute sine and cosine of Phi_p - Phi_0
1954  gammalib::sincosd(phip - m_phi0, &sphip, &cphip);
1955 
1956  // ...
1957  double x = cthe0 * cphip;
1958  double y = sthe0;
1959  double z = std::sqrt(x*x + y*y);
1960  if (z == 0.0) {
1961 
1962  // Check of an invalid coordinate transformation parameter
1963  // has been encountered. Since z=0 we exlect sin(lat0)=0.
1964  if (slat0 != 0.0) {
1965  std::string message = "sin(lat0)="+gammalib::str(slat0)+", expected 0.";
1967  }
1968 
1969  // latp determined solely by LATPOLEa in this case
1970  m_latpreq = 2;
1971  if (latp > 90.0) {
1972  latp = 90.0;
1973  }
1974  else if (latp < -90.0) {
1975  latp = -90.0;
1976  }
1977 
1978  }
1979 
1980  // ... otherwise ...
1981  else {
1982  double slz = slat0/z;
1983  if (std::abs(slz) > 1.0) {
1984  if ((std::abs(slz) - 1.0) < tol) {
1985  if (slz > 0.0) {
1986  slz = 1.0;
1987  }
1988  else {
1989  slz = -1.0;
1990  }
1991  }
1992  else {
1993  std::string message;
1994  message = "abs(slz)-1 >= "+gammalib::str(std::abs(slz) - 1.0);
1995  message += +", expected <"+gammalib::str(tol)+".";
1997  }
1998  }
1999  u = gammalib::atan2d(y,x);
2000  v = gammalib::acosd(slz);
2001  } // endelse: z != 0
2002  } // endelse: ...
2003 
2004  // ...
2005  if (m_latpreq == 0) {
2006  double latp1 = u + v;
2007  if (latp1 > 180.0) {
2008  latp1 -= 360.0;
2009  }
2010  else if (latp1 < -180.0) {
2011  latp1 += 360.0;
2012  }
2013 
2014  // ...
2015  double latp2 = u - v;
2016  if (latp2 > 180.0) {
2017  latp2 -= 360.0;
2018  }
2019  else if (latp2 < -180.0) {
2020  latp2 += 360.0;
2021  }
2022 
2023  // Are there two valid solutions for latp?
2024  if (std::abs(latp1) < 90.0+tol && std::abs(latp2) < 90.0+tol) {
2025  m_latpreq = 1;
2026  }
2027 
2028  // Select solution
2029  if (std::abs(latp-latp1) < std::abs(latp-latp2)) {
2030  latp = (std::abs(latp1) < 90.0+tol) ? latp1 : latp2;
2031  }
2032  else {
2033  latp = (std::abs(latp2) < 90.0+tol) ? latp2 : latp1;
2034  }
2035 
2036  // Account for rounding error
2037  if (std::abs(latp) < 90.0+tol) {
2038  if (latp > 90.0) {
2039  latp = 90.0;
2040  }
2041  else if (latp < -90.0) {
2042  latp = -90.0;
2043  }
2044  }
2045  } // endif: ...
2046 
2047  // ...
2048  double z = gammalib::cosd(latp) * clat0;
2049  if (std::abs(z) < tol) {
2050 
2051  // Celestial pole at the fiducial point
2052  if (std::abs(clat0) < tol) {
2053  lngp = lng0;
2054  }
2055 
2056  // Celestial north pole at the native pole
2057  else if (latp > 0.0) {
2058  lngp = lng0 + phip - m_phi0 - 180.0;
2059  }
2060 
2061  // Celestial south pole at the native pole
2062  else {
2063  lngp = lng0 - phip + m_phi0;
2064  }
2065 
2066  }
2067 
2068  // ...
2069  else {
2070  double x = (sthe0 - gammalib::sind(latp)*slat0)/z;
2071  double y = sphip*cthe0/clat0;
2072  if (x == 0.0 && y == 0.0) {
2073  std::string message;
2074  message = "x=0 and y=0";
2075  message += +", expected at least that one differs from 0.";
2077  }
2078  lngp = lng0 - gammalib::atan2d(y,x);
2079  }
2080 
2081  // Make celestial longitude at the pole the same sign as at the
2082  // fiducial point
2083  if (lng0 >= 0.0) {
2084  if (lngp < 0.0) {
2085  lngp += 360.0;
2086  }
2087  else if (lngp > 360.0) {
2088  lngp -= 360.0;
2089  }
2090  }
2091  else {
2092  if (lngp > 0.0) {
2093  lngp -= 360.0;
2094  }
2095  else if (lngp < -360.0) {
2096  lngp += 360.0;
2097  }
2098  }
2099 
2100  } // endelse: fiducial point was away from native pole
2101 
2102  // Store LATPOLEa
2103  m_ref[3] = latp;
2104 
2105  // Set the Euler angles
2106  m_euler[0] = lngp;
2107  m_euler[1] = 90.0 - latp;
2108  m_euler[2] = phip;
2109  gammalib::sincosd(m_euler[1], &m_euler[4], &m_euler[3]);
2110 
2111  // Signal if |latitude| is preserved
2112  m_isolat = (m_euler[4] == 0.0);
2113 
2114  // Check for ill-conditioned parameters
2115  if (std::abs(latp) > 90.0+tol) {
2116  std::string message;
2117  message = "abs(latp) > 90, coordinate transformation is ill-conditioned.";
2119  }
2120 
2121  // Celestial parameter have been set
2122  m_celset = true;
2123 
2124  // Return
2125  return;
2126 }
2127 
2128 
2129 /***********************************************************************//**
2130  * @brief Pixel-to-world celestial transformation
2131  *
2132  * @param[in] nx X pixel vector length.
2133  * @param[in] ny Y pixel vector length (0=no replication, ny=nx).
2134  * @param[in] sxy Input vector step.
2135  * @param[in] sll Output vector step.
2136  * @param[in] x Vector of projected x coordinates.
2137  * @param[in] y Vector of projected y coordinates.
2138  * @param[out] phi Longitude of the projected point in native spherical
2139  * coordinates [deg].
2140  * @param[out] theta Latitude of the projected point in native spherical
2141  * coordinates [deg].
2142  * @param[out] lng Celestial longitude of the projected point [deg].
2143  * @param[out] lat Celestial latitude of the projected point [deg].
2144  * @param[out] stat Status return value for each vector element
2145  * (0=success, 1=invalid).
2146  *
2147  * This method transforms (x,y) coordinates in the plane of projection to
2148  * celestial coordinates (lng,lat). The method has been adapted from the
2149  * wcslib function cel.c::celx2s.
2150  ***************************************************************************/
2151 void GWcs::cel_x2s(int nx, int ny, int sxy, int sll,
2152  const double* x, const double *y,
2153  double* phi, double* theta,
2154  double* lng, double* lat, int* stat) const
2155 {
2156  // Initialize celestial transformations if required
2157  if (!m_celset) {
2158  cel_set();
2159  }
2160 
2161  // Apply spherical deprojection
2162  prj_x2s(nx, ny, sxy, 1, x, y, phi, theta, stat);
2163 
2164  // Compute number of Phi values
2165  int nphi = (ny > 0) ? (nx*ny) : nx;
2166 
2167  // Compute celestial coordinates
2168  sph_x2s(nphi, 0, 1, sll, phi, theta, lng, lat);
2169 
2170  // Return
2171  return;
2172 }
2173 
2174 
2175 /***********************************************************************//**
2176  * @brief World-to-pixel celestial transformation
2177  *
2178  * @param[in] nlng Longitude vector length.
2179  * @param[in] nlat Latitude vector length (0=no replication, nlat=nlng).
2180  * @param[in] sll Input vector step.
2181  * @param[in] sxy Output vector step.
2182  * @param[in] lng Celestial longitude of the projected point [deg].
2183  * @param[in] lat Celestial latitude of the projected point [deg].
2184  * @param[out] phi Longitude of the projected point in native spherical
2185  * coordinates [deg].
2186  * @param[out] theta Latitude of the projected point in native spherical
2187  * coordinates [deg].
2188  * @param[out] x Vector of projected x coordinates.
2189  * @param[out] y Vector of projected y coordinates.
2190  * @param[out] stat Status return value for each vector element
2191  * (0=success, 1=invalid).
2192  *
2193  * This method transforms (x,y) coordinates in the plane of projection to
2194  * celestial coordinates (lng,lat). The method has been adapted from the
2195  * wcslib function cel.c::celx2s.
2196  ***************************************************************************/
2197 void GWcs::cel_s2x(int nlng, int nlat, int sll, int sxy,
2198  const double* lng, const double* lat,
2199  double* phi, double* theta,
2200  double* x, double* y, int* stat) const
2201 {
2202  // Initialize celestial transformations if required
2203  if (!m_celset) {
2204  cel_set();
2205  }
2206 
2207  // Compute native coordinates
2208  sph_s2x(nlng, nlat, sll, 1, lng, lat, phi, theta);
2209 
2210  // Constant celestial latitude -> constant native latitude
2211  int nphi;
2212  int ntheta;
2213  if (m_isolat) {
2214  nphi = nlng;
2215  ntheta = nlat;
2216  }
2217  else {
2218  nphi = (nlat > 0) ? (nlng*nlat) : nlng;
2219  ntheta = 0;
2220  }
2221 
2222  // Apply spherical projection
2223  prj_s2x(nphi, ntheta, 1, sxy, phi, theta, x, y, stat);
2224 
2225  // Return
2226  return;
2227 }
2228 
2229 
2230 /*==========================================================================
2231  = =
2232  = Spherical transformation methods adapted from wcslib =
2233  = =
2234  ==========================================================================*/
2235 
2236 /***********************************************************************//**
2237  * @brief Rotation in the pixel-to-world direction
2238  *
2239  * @param[in] nphi Phi vector length.
2240  * @param[in] ntheta Theta vector length (0=no replication).
2241  * @param[in] spt Input vector step.
2242  * @param[in] sll Output vector step.
2243  * @param[in] phi Longitude in the native coordinate system of the
2244  * projection [deg].
2245  * @param[in] theta Latitude in the native coordinate system of the
2246  * projection [deg].
2247  * @param[out] lng Celestial longitude [deg].
2248  * @param[out] lat Celestial latitude [deg].
2249  *
2250  * This method has been adapted from the wcslib function sph.c::sphx2s().
2251  * The interface follows very closely that of wcslib.
2252  ***************************************************************************/
2253 void GWcs::sph_x2s(int nphi, int ntheta, int spt, int sll,
2254  const double* phi, const double* theta,
2255  double* lng, double* lat) const
2256 {
2257  // Set tolerance
2258  const double tol = 1.0e-5;
2259 
2260  // Set value replication length mphi,mtheta
2261  int mphi;
2262  int mtheta;
2263  if (ntheta > 0) {
2264  mphi = nphi;
2265  mtheta = ntheta;
2266  }
2267  else {
2268  mphi = 1;
2269  mtheta = 1;
2270  ntheta = nphi;
2271  }
2272 
2273  // Check for a simple change in origin of longitude
2274  if (m_euler[4] == 0.0) {
2275 
2276  // Initialise pointers
2277  double* lngp = lng;
2278  double* latp = lat;
2279  const double* phip = phi;
2280  const double* thetap = theta;
2281 
2282  // Case A: ...
2283  if (m_euler[1] == 0.0) {
2284 
2285  // ...
2286  double dlng = fmod(m_euler[0] + 180.0 - m_euler[2], 360.0);
2287 
2288  // ...
2289  for (int itheta = 0; itheta < ntheta; ++itheta, phip += spt, thetap += spt) {
2290  for (int iphi = 0; iphi < mphi; ++iphi, lngp += sll, latp += sll) {
2291 
2292  // Shift longitude
2293  *lngp = *phip + dlng;
2294  *latp = *thetap;
2295 
2296  // Normalize the celestial longitude
2297  if (m_euler[0] >= 0.0) {
2298  if (*lngp < 0.0) {
2299  *lngp += 360.0;
2300  }
2301  }
2302  else {
2303  if (*lngp > 0.0) {
2304  *lngp -= 360.0;
2305  }
2306  }
2307  if (*lngp > 360.0) {
2308  *lngp -= 360.0;
2309  }
2310  else if (*lngp < -360.0) {
2311  *lngp += 360.0;
2312  }
2313 
2314  } // endfor: looped over phi
2315  } // endfor: looped over theta
2316 
2317  } // endif: ...
2318 
2319  // Case B: ...
2320  else {
2321 
2322  // ...
2323  double dlng = fmod(m_euler[0] + m_euler[2], 360.0);
2324 
2325  // ...
2326  for (int itheta = 0; itheta < ntheta; ++itheta, phip += spt, thetap += spt) {
2327  for (int iphi = 0; iphi < mphi; ++iphi, lngp += sll, latp += sll) {
2328 
2329  // Shift longitude and inverte latitude
2330  *lngp = dlng - *phip;
2331  *latp = -(*thetap);
2332 
2333  // Normalize the celestial longitude
2334  if (m_euler[0] >= 0.0) {
2335  if (*lngp < 0.0) {
2336  *lngp += 360.0;
2337  }
2338  }
2339  else {
2340  if (*lngp > 0.0) {
2341  *lngp -= 360.0;
2342  }
2343  }
2344  if (*lngp > 360.0) {
2345  *lngp -= 360.0;
2346  }
2347  else if (*lngp < -360.0) {
2348  *lngp += 360.0;
2349  }
2350 
2351  } // endfor: looped over phi
2352  } // endfor: looped over theta
2353  } // endelse: ...
2354 
2355  } // endif: there was a simple change of longitude
2356 
2357  // ... otherwise to general transformation
2358  else {
2359 
2360  // Do phi dependency
2361  const double* phip = phi;
2362  int rowoff = 0;
2363  int rowlen = nphi*sll;
2364  for (int iphi = 0; iphi < nphi; ++iphi, rowoff += sll, phip += spt) {
2365  double dphi = *phip - m_euler[2];
2366  double* lngp = lng + rowoff;
2367  for (int itheta = 0; itheta < mtheta; ++itheta, lngp += rowlen) {
2368  *lngp = dphi;
2369  }
2370  }
2371 
2372  // Do theta dependency
2373  const double* thetap = theta;
2374  double* lngp = lng;
2375  double* latp = lat;
2376  for (int itheta = 0; itheta < ntheta; ++itheta, thetap += spt) {
2377 
2378  // ...
2379  double sinthe;
2380  double costhe;
2381  gammalib::sincosd(*thetap, &sinthe, &costhe);
2382 
2383  // ...
2384  double costhe3 = costhe * m_euler[3];
2385  double costhe4 = costhe * m_euler[4];
2386  double sinthe3 = sinthe * m_euler[3];
2387  double sinthe4 = sinthe * m_euler[4];
2388 
2389  // Loop over Phi
2390  for (int iphi = 0; iphi < mphi; ++iphi, lngp += sll, latp += sll) {
2391 
2392  //
2393  double dphi = *lngp;
2394  double sinphi;
2395  double cosphi;
2396  gammalib::sincosd(dphi, &sinphi, &cosphi);
2397 
2398  // Compute the celestial longitude
2399  double x = sinthe4 - costhe3*cosphi;
2400  double y = -costhe*sinphi;
2401 
2402  // Rearrange longitude formula to reduce roundoff errors
2403  if (std::abs(x) < tol) {
2404  x = -gammalib::cosd(*thetap + m_euler[1]) +
2405  costhe3*(1.0 - cosphi);
2406  }
2407 
2408  // Compute longitude shift
2409  double dlng;
2410  if (x != 0.0 || y != 0.0) {
2411  dlng = gammalib::atan2d(y, x);
2412  }
2413  else {
2414  dlng = (m_euler[1] < 90.0) ? dphi + 180.0 : -dphi;
2415  }
2416 
2417  // Set celestial longitude
2418  *lngp = m_euler[0] + dlng;
2419 
2420  // Normalize the celestial longitude
2421  if (m_euler[0] >= 0.0) {
2422  if (*lngp < 0.0) {
2423  *lngp += 360.0;
2424  }
2425  }
2426  else {
2427  if (*lngp > 0.0) {
2428  *lngp -= 360.0;
2429  }
2430  }
2431  if (*lngp > 360.0) {
2432  *lngp -= 360.0;
2433  }
2434  else if (*lngp < -360.0) {
2435  *lngp += 360.0;
2436  }
2437 
2438  // Compute the celestial latitude. First handle the case
2439  // of longitude shifts by 180 deg
2440  if (fmod(dphi,180.0) == 0.0) {
2441  *latp = *thetap + cosphi*m_euler[1];
2442  if (*latp > 90.0) *latp = 180.0 - *latp;
2443  if (*latp < -90.0) *latp = -180.0 - *latp;
2444  }
2445 
2446  // ... then handle the case of general longitude shifts
2447  else {
2448  // Use alternative formulae for greater accuracy
2449  double z = sinthe3 + costhe4*cosphi;
2450  if (z > 0.99) {
2451  *latp = gammalib::acosd(sqrt(x*x+y*y));
2452  }
2453  else if (z < -0.99) {
2454  *latp = -gammalib::acosd(sqrt(x*x+y*y));
2455  }
2456  else {
2457  *latp = gammalib::asind(z);
2458  }
2459  } // endelse: general longitude shift
2460 
2461  } // endfor: looped over phi
2462 
2463  } // endfor: looped over theta
2464 
2465  } // endelse: did general transformation
2466 
2467  // Return
2468  return;
2469 }
2470 
2471 
2472 /***********************************************************************//**
2473  * @brief Rotation in the pixel-to-world direction
2474  *
2475  * @param[in] nlng Longitude vector length.
2476  * @param[in] nlat Latitude vector length (0=no replication).
2477  * @param[in] sll Input vector step.
2478  * @param[in] spt Output vector step.
2479  * @param[in] lng Celestial longitude [deg].
2480  * @param[in] lat Celestial latitude [deg].
2481  * @param[out] phi Longitude in the native coordinate system of the
2482  * projection [deg].
2483  * @param[out] theta Latitude in the native coordinate system of the
2484  * projection [deg].
2485  *
2486  * This method has been adapted from the wcslib function sph.c::sphs2x().
2487  * The interface follows very closely that of wcslib.
2488  ***************************************************************************/
2489 void GWcs::sph_s2x(int nlng, int nlat, int sll, int spt,
2490  const double* lng, const double* lat,
2491  double* phi, double* theta) const
2492 {
2493  // Set tolerance
2494  const double tol = 1.0e-5;
2495 
2496  // Set value replication length mlng,mlat
2497  int mlng;
2498  int mlat;
2499  if (nlat > 0) {
2500  mlng = nlng;
2501  mlat = nlat;
2502  }
2503  else {
2504  mlng = 1;
2505  mlat = 1;
2506  nlat = nlng;
2507  }
2508 
2509  // Check for a simple change in origin of longitude
2510  if (m_euler[4] == 0.0) {
2511 
2512  // Initialise pointers
2513  const double* lngp = lng;
2514  const double* latp = lat;
2515  double* phip = phi;
2516  double* thetap = theta;
2517 
2518  // Case A: ...
2519  if (m_euler[1] == 0.0) {
2520 
2521  // Compute longitude shift
2522  double dphi = fmod(m_euler[2] - 180.0 - m_euler[0], 360.0);
2523 
2524  // Apply longitude shift
2525  for (int ilat = 0; ilat < nlat; ++ilat, lngp += sll, latp += sll) {
2526  for (int ilng = 0; ilng < mlng; ++ilng, phip += spt, thetap += spt) {
2527 
2528  // Shift longitude, keep Theta
2529  *phip = fmod(*lngp + dphi, 360.0);
2530  *thetap = *latp;
2531 
2532  // Normalize the native longitude
2533  if (*phip > 180.0) {
2534  *phip -= 360.0;
2535  }
2536  else if (*phip < -180.0) {
2537  *phip += 360.0;
2538  }
2539 
2540  } // endfor: looped over longitude
2541  } // endfor: looped over latitude
2542 
2543  } // endif: Case A
2544 
2545  // Case B: ...
2546  else {
2547 
2548  // Compute longitude shift
2549  double dphi = fmod(m_euler[2] + m_euler[0], 360.0);
2550 
2551  // Apply longitude shift
2552  for (int ilat = 0; ilat < nlat; ++ilat, lngp += sll, latp += sll) {
2553  for (int ilng = 0; ilng < mlng; ++ilng, phip += spt, thetap += spt) {
2554 
2555  // Shift longitude, flip Theta
2556  *phip = fmod(dphi - *lngp, 360.0);
2557  *thetap = -(*latp);
2558 
2559  // Normalize the native longitude
2560  if (*phip > 180.0) {
2561  *phip -= 360.0;
2562  }
2563  else if (*phip < -180.0) {
2564  *phip += 360.0;
2565  }
2566 
2567  } // endfor: looped over longitude
2568  } // endfor: looped over latitude
2569 
2570  } // endelse: Case B
2571 
2572  } // endif: simple change in origin
2573 
2574  // ... otherwise compute general transformation
2575  else {
2576 
2577  // Do lng dependency
2578  const double* lngp = lng;
2579  int rowoff = 0;
2580  int rowlen = nlng * spt;
2581  for (int ilng = 0; ilng < nlng; ++ilng, rowoff += spt, lngp += sll) {
2582  double dlng = *lngp - m_euler[0];
2583  double* phip = phi + rowoff;
2584  for (int ilat = 0; ilat < mlat; ++ilat, phip += rowlen) {
2585  *phip = dlng;
2586  }
2587  }
2588 
2589  // Do lat dependency
2590  const double* latp = lat;
2591  double* phip = phi;
2592  double* thetap = theta;
2593  for (int ilat = 0; ilat < nlat; ++ilat, latp += sll) {
2594 
2595  // ...
2596  double sinlat;
2597  double coslat;
2598  gammalib::sincosd(*latp, &sinlat, &coslat);
2599  double coslat3 = coslat*m_euler[3];
2600  double coslat4 = coslat*m_euler[4];
2601  double sinlat3 = sinlat*m_euler[3];
2602  double sinlat4 = sinlat*m_euler[4];
2603 
2604  // Loop over longitudes
2605  for (int ilng = 0; ilng < mlng; ++ilng, phip += spt, thetap += spt) {
2606 
2607  // ...
2608  double dlng = *phip;
2609  double sinlng;
2610  double coslng;
2611  gammalib::sincosd(dlng, &sinlng, &coslng);
2612 
2613  // Compute the native longitude
2614  double x = sinlat4 - coslat3*coslng;
2615  double y = -coslat*sinlng;
2616 
2617  // Rearrange formula to reduce roundoff errors
2618  if (std::abs(x) < tol) {
2619  x = -gammalib::cosd(*latp+m_euler[1]) +
2620  coslat3*(1.0 - coslng);
2621  }
2622 
2623  // Compute Phi shift
2624  double dphi;
2625  if (x != 0.0 || y != 0.0) {
2626  dphi = gammalib::atan2d(y, x);
2627  }
2628  else { // Change of origin of longitude
2629  if (m_euler[1] < 90.0) {
2630  dphi = dlng - 180.0;
2631  }
2632  else {
2633  dphi = -dlng;
2634  }
2635  }
2636 
2637  // Set Phi
2638  *phip = fmod(m_euler[2] + dphi, 360.0);
2639 
2640  // Normalize the native longitude
2641  if (*phip > 180.0) {
2642  *phip -= 360.0;
2643  }
2644  else if (*phip < -180.0) {
2645  *phip += 360.0;
2646  }
2647 
2648  // Compute the native latitude. First handle the case
2649  // of longitude shifts by 180 deg
2650  if (fmod(dlng, 180.0) == 0.0) {
2651  *thetap = *latp + coslng*m_euler[1];
2652  if (*thetap > 90.0) *thetap = 180.0 - *thetap;
2653  if (*thetap < -90.0) *thetap = -180.0 - *thetap;
2654  }
2655 
2656  // ... then handle the case of general longitude shifts
2657  else {
2658  // Use alternative formulae for greater accuracy
2659  double z = sinlat3 + coslat4*coslng;
2660  if (z > 0.99) {
2661  *thetap = gammalib::acosd(sqrt(x*x+y*y));
2662  }
2663  else if (z < -0.99) {
2664  *thetap = -gammalib::acosd(sqrt(x*x+y*y));
2665  }
2666  else {
2667  *thetap = gammalib::asind(z);
2668  }
2669  }
2670 
2671  } // endfor: looped over longitudes
2672 
2673  } // endfor: looped over latitude
2674 
2675  } // endelse: handled general case
2676 
2677  // Return
2678  return;
2679 }
2680 
2681 
2682 /*==========================================================================
2683  = =
2684  = Spectral methods adapted from wcslib =
2685  = =
2686  ==========================================================================*/
2687 
2688 /***********************************************************************//**
2689  * @brief Initialise spectral transformation parameters
2690  *
2691  * Code adapted from spc.c::spcini(). The actual version of the code does
2692  * nothing.
2693  ***************************************************************************/
2694 void GWcs::spc_ini(void)
2695 {
2696  // Return
2697  return;
2698 }
2699 
2700 
2701 /*==========================================================================
2702  = =
2703  = Linear transformation methods adapted from wcslib =
2704  = =
2705  ==========================================================================*/
2706 
2707 /***********************************************************************//**
2708  * @brief Initialise linear transformation parameters
2709  *
2710  * Code adapted from lin.c::linini().
2711  ***************************************************************************/
2712 void GWcs::lin_ini(int naxis)
2713 {
2714  // Initialise parameters
2715  m_linset = false;
2716  m_unity = false;
2717 
2718  // Clear vectors
2719  m_crpix.clear();
2720  m_pc.clear();
2721  m_cdelt.clear();
2722  m_piximg.clear();
2723  m_imgpix.clear();
2724 
2725  // Continue only if there are axes
2726  if (naxis > 0) {
2727 
2728  // Loop over axes
2729  for (int i = 0; i < naxis; ++i) {
2730 
2731  // CRPIXja defaults to 0.0
2732  m_crpix.push_back(0.0);
2733 
2734  // PCi_ja defaults to the unit matrix
2735  for (int j = 0; j < naxis; ++j) {
2736  if (j == i) {
2737  m_pc.push_back(1.0);
2738  }
2739  else {
2740  m_pc.push_back(0.0);
2741  }
2742  }
2743 
2744  // CDELTia defaults to 1.0
2745  m_cdelt.push_back(1.0);
2746 
2747  } // endfor: looped over axes
2748 
2749  } // endif: there were axes
2750 
2751  // Return
2752  return;
2753 }
2754 
2755 
2756 /***********************************************************************//**
2757  * @brief Initialise linear transformation parameters
2758  *
2759  * Code adapted from lin.c::linset(). It requires m_pc to be set and to be
2760  * of size m_naxis*m_naxis.
2761  ***************************************************************************/
2762 void GWcs::lin_set(void) const
2763 {
2764  // Clear transformation matrices
2765  m_piximg.clear();
2766  m_imgpix.clear();
2767 
2768  // Check for unity PC matrix
2769  m_unity = true;
2770  for (int i = 0, index = 0; i < m_naxis; ++i) {
2771  for (int j = 0; j < m_naxis; ++j, ++index) {
2772  if (j == i) {
2773  if (m_pc[index] != 1.0) {
2774  m_unity = false;
2775  break;
2776  }
2777  }
2778  else {
2779  if (m_pc[index] != 0.0) {
2780  m_unity = false;
2781  break;
2782  }
2783  }
2784  }
2785  if (!m_unity) {
2786  break;
2787  }
2788  }
2789 
2790  // Debug option: force PC usage for testing purposes
2791  #if defined(G_LIN_MATINV_FORCE_PC)
2792  m_unity = false;
2793  std::cout << "DEBUG: Force PC usage in linear transformations." << std::endl;
2794  #endif
2795 
2796  // Compute transformation matrices for non-uniform PC matrix
2797  if (!m_unity) {
2798 
2799  // Compute the pixel-to-image transformation matrix
2800  for (int i = 0, index = 0; i < m_naxis; ++i) {
2801  for (int j = 0; j < m_naxis; ++j, ++index) {
2802  m_piximg.push_back(m_cdelt[i] * m_pc[index]);
2803  }
2804  }
2805 
2806  // Compute the image-to-pixel transformation matrix
2808 
2809  }
2810 
2811  // Signal that linear transformation matrix has been set
2812  m_linset = true;
2813 
2814  // Return
2815  return;
2816 }
2817 
2818 
2819 /***********************************************************************//**
2820  * @brief Pixel-to-world linear transformation
2821  *
2822  * @param[in] ncoord Number of coordinates.
2823  * @param[in] nelem Vector length of each coordinate (>=m_naxis).
2824  * @param[in] pixcrd Array of pixel coordinates (ncoord*nelem).
2825  * @param[out] imgcrd Array of intermediate world coordinates (ncoord*nelem).
2826  *
2827  * Transforms pixel coordinates to intermediate world coordinates. This
2828  * method has been adapted from lin.c::linp2x(). The method performs distinct
2829  * computations depending of whether the PC matrix is unity or not, avoiding
2830  * a matrix multiplication when it is not necessary.
2831  ***************************************************************************/
2832 void GWcs::lin_p2x(int ncoord, int nelem, const double* pixcrd, double* imgcrd) const
2833 {
2834  // Initialize linear transformations if required
2835  if (!m_linset) {
2836  lin_set();
2837  }
2838 
2839  // Convert pixel coordinates to intermediate world coordinates
2840  const double* pix = pixcrd;
2841  double* img = imgcrd;
2842 
2843  // Case A: we have a unity PC matrix
2844  if (m_unity) {
2845 
2846  // Loop over all coordinates
2847  for (int k = 0; k < ncoord; ++k) {
2848 
2849  // Transform the first m_naxis elements
2850  for (int i = 0; i < m_naxis; ++i) {
2851  *(img++) = m_cdelt[i] * (*(pix++) - m_crpix[i]);
2852  }
2853 
2854  // Go to next coordinate
2855  pix += (nelem - m_naxis);
2856  img += (nelem - m_naxis);
2857 
2858  } // endfor: looped over all coordinates
2859 
2860  }
2861 
2862  // Case B: we have a non-unity PC matrix
2863  else {
2864 
2865  // Loop over all coordinates
2866  for (int k = 0; k < ncoord; ++k) {
2867 
2868  // Clear first m_naxis elements
2869  for (int i = 0; i < m_naxis; ++i) {
2870  img[i] = 0.0;
2871  }
2872 
2873  // Perform matrix multiplication (column-wise multiplication
2874  // allows this to be cached)
2875  for (int j = 0; j < m_naxis; ++j) {
2876  double temp = *(pix++) - m_crpix[j];
2877  for (int i = 0, ji = j; i < m_naxis; ++i, ji += m_naxis) {
2878  img[i] += m_piximg[ji] * temp;
2879  }
2880  }
2881 
2882  // Go to next coordinate
2883  pix += (nelem - m_naxis);
2884  img += nelem;
2885 
2886  } // endfor: looped over all coordinates
2887 
2888  } // endelse: PC matrix was not unity
2889 
2890  // Return
2891  return;
2892 }
2893 
2894 
2895 /***********************************************************************//**
2896  * @brief World-to-pixel linear transformation
2897  *
2898  * @param[in] ncoord Number of coordinates.
2899  * @param[in] nelem Vector length of each coordinate (>=m_naxis).
2900  * @param[in] imgcrd Array of intermediate world coordinates (ncoord*nelem).
2901  * @param[out] pixcrd Array of pixel coordinates (ncoord*nelem).
2902  *
2903  * Transforms intermediate world coordinates to pixel coordinates. This
2904  * method has been adapted from lin.c::linx2p(). The method performs distinct
2905  * computations depending of whether the PC matrix is unity or not, avoiding
2906  * a matrix multiplication when it is not necessary.
2907  ***************************************************************************/
2908 void GWcs::lin_x2p(int ncoord, int nelem, const double* imgcrd, double* pixcrd) const
2909 {
2910  // Initialize linear transformations if required
2911  if (!m_linset) {
2912  lin_set();
2913  }
2914 
2915  // Convert pixel coordinates to intermediate world coordinates
2916  const double* img = imgcrd;
2917  double* pix = pixcrd;
2918 
2919  // Case A: we have a unity PC matrix
2920  if (m_unity) {
2921 
2922  // Loop over all coordinates
2923  for (int k = 0; k < ncoord; ++k) {
2924 
2925  // Transform the first m_naxis elements
2926  for (int i = 0; i < m_naxis; ++i) {
2927  *(pix++) = (*(img++) / m_cdelt[i]) + m_crpix[i];
2928  }
2929 
2930  // Go to next coordinate
2931  pix += (nelem - m_naxis);
2932  img += (nelem - m_naxis);
2933 
2934  } // endfor: looped over all coordinates
2935 
2936  }
2937 
2938  // Case B: we have a non-unity PC matrix
2939  else {
2940 
2941  // Loop over all coordinates
2942  for (int k = 0; k < ncoord; ++k) {
2943 
2944  // Perform matrix multiplication
2945  for (int j = 0, ji = 0; j < m_naxis; ++j) {
2946  *pix = 0.0;
2947  for (int i = 0; i < m_naxis; ++i, ++ji) {
2948  *pix += m_imgpix[ji] * img[i];
2949  }
2950  *(pix++) += m_crpix[j];
2951  }
2952 
2953  // Go to next coordinate
2954  pix += (nelem - m_naxis);
2955  img += nelem;
2956 
2957  } // endfor: looped over all coordinates
2958 
2959  } // endelse: PC matrix was not unity
2960 
2961  // Return
2962  return;
2963 }
2964 
2965 
2966 /***********************************************************************//**
2967  * @brief Invert linear transformation matrix
2968  *
2969  * @param[in] mat Matrix
2970  * @param[out] inv Inverted matrix
2971  *
2972  * @exception GException::wcs_singular_matrix
2973  * Singular matrix encountered.
2974  *
2975  * Inverts a linear transformation matrix that is stored in a vector. The
2976  * matrix dimension is assumed to be m_naxis*m_naxis. This method has been
2977  * adapted from lin.c::matinv().
2978  ***************************************************************************/
2979 void GWcs::lin_matinv(const std::vector<double>& mat, std::vector<double>& inv) const
2980 {
2981  // Continue only if naxis is valid
2982  if (m_naxis > 0) {
2983 
2984  // Declare internal arrays
2985  std::vector<int> mxl;
2986  std::vector<int> lxm;
2987  std::vector<double> rowmax;
2988  std::vector<double> lu = mat;
2989 
2990  // Clear inverse matrix
2991  inv.clear();
2992  inv.assign(m_naxis*m_naxis, 0.0);
2993 
2994  // Initialize arrays
2995  for (int i = 0, ij = 0; i < m_naxis; ++i) {
2996 
2997  // Vector that records row interchanges
2998  mxl.push_back(i);
2999  lxm.push_back(0);
3000 
3001  // Maximum element value in row
3002  rowmax.push_back(0.0);
3003  for (int j = 0; j < m_naxis; ++j, ++ij) {
3004  double dtemp = std::abs(mat[ij]);
3005  if (dtemp > rowmax[i]) {
3006  rowmax[i] = dtemp;
3007  }
3008  }
3009 
3010  // Throw exception if matrix is singular
3011  if (rowmax[i] == 0.0) {
3012  throw GException::wcs_singular_matrix(G_LIN_MATINV, m_naxis, mat);
3013  }
3014 
3015  } // endfor: initialized array
3016 
3017  // Form the LU triangular factorization using scaled partial pivoting
3018  for (int k = 0; k < m_naxis; ++k) {
3019 
3020  // Decide whether to pivot
3021  double colmax = std::abs(lu[k*m_naxis+k]) / rowmax[k];
3022  int pivot = k;
3023  for (int i = k+1; i < m_naxis; ++i) {
3024  int ik = i*m_naxis + k;
3025  double dtemp = std::abs(lu[ik]) / rowmax[i];
3026  if (dtemp > colmax) {
3027  colmax = dtemp;
3028  pivot = i;
3029  }
3030  }
3031 
3032  // Do we need to pivot
3033  if (pivot > k) {
3034 
3035  // We must pivot, interchange the rows of the design matrix
3036  for (int j = 0, pj = pivot*m_naxis, kj = k*m_naxis; j < m_naxis; ++j, ++pj, ++kj) {
3037  double dtemp = lu[pj];
3038  lu[pj] = lu[kj];
3039  lu[kj] = dtemp;
3040  }
3041 
3042  // Amend the vector of row maxima
3043  double dtemp = rowmax[pivot];
3044  rowmax[pivot] = rowmax[k];
3045  rowmax[k] = dtemp;
3046 
3047  // Record the interchange for later use
3048  int itemp = mxl[pivot];
3049  mxl[pivot] = mxl[k];
3050  mxl[k] = itemp;
3051 
3052  } // endif: pivoting required
3053 
3054  // Gaussian elimination
3055  for (int i = k+1; i < m_naxis; ++i) {
3056 
3057  // Compute matrix index
3058  int ik = i*m_naxis + k;
3059 
3060  // Nothing to do if lu[ik] is zero
3061  if (lu[ik] != 0.0) {
3062 
3063  // Save the scaling factor
3064  lu[ik] /= lu[k*m_naxis+k];
3065 
3066  // Subtract rows
3067  for (int j = k+1; j < m_naxis; ++j) {
3068  lu[i*m_naxis+j] -= lu[ik]*lu[k*m_naxis+j];
3069  }
3070 
3071  } // endif: lu[ik] was not zero
3072 
3073  } // endfor: Gaussian elimination
3074 
3075  } // endfor: formed LU triangular factorization
3076 
3077  // mxl[i] records which row of mat corresponds to row i of lu
3078  // lxm[i] records which row of lu corresponds to row i of mat
3079  for (int i = 0; i < m_naxis; ++i) {
3080  lxm[mxl[i]] = i;
3081  }
3082 
3083  // Determine the inverse matrix
3084  for (int k = 0; k < m_naxis; ++k) {
3085 
3086  // ...
3087  inv[lxm[k]*m_naxis+k] = 1.0;
3088 
3089  // Forward substitution
3090  for (int i = lxm[k]+1; i < m_naxis; ++i) {
3091  for (int j = lxm[k]; j < i; ++j) {
3092  inv[i*m_naxis+k] -= lu[i*m_naxis+j]*inv[j*m_naxis+k];
3093  }
3094  }
3095 
3096  // Backward substitution
3097  for (int i = m_naxis-1; i >= 0; --i) {
3098  for (int j = i+1; j < m_naxis; ++j) {
3099  inv[i*m_naxis+k] -= lu[i*m_naxis+j]*inv[j*m_naxis+k];
3100  }
3101  inv[i*m_naxis+k] /= lu[i*m_naxis+i];
3102  }
3103 
3104  } // endfor: determined inverse matrix
3105 
3106  } // endif: naxis was valid
3107 
3108  // Return
3109  return;
3110 }
3111 
3112 
3113 /*==========================================================================
3114  = =
3115  = Projection methods adapted from wcslib =
3116  = =
3117  ==========================================================================*/
3118 
3119 /***********************************************************************//**
3120  * @brief Initialise projection parameters
3121  *
3122  * Code adapted from prj.c::prjini().
3123  ***************************************************************************/
3124 void GWcs::prj_ini(void) const
3125 {
3126  // Initialise parameters
3127  m_prjset = false;
3128  m_r0 = 0.0;
3129  m_pv[0] = 0.0;
3130  m_pv[1] = UNDEFINED;
3131  m_pv[2] = UNDEFINED;
3132  m_pv[3] = UNDEFINED;
3133  for (int k = 4; k < PVN; ++k) {
3134  m_pv[k] = 0.0;
3135  }
3136  m_bounds = true; // 1 in wcslib
3137  m_x0 = 0.0;
3138  m_y0 = 0.0;
3139  m_w.clear();
3140 
3141  // Return
3142  return;
3143 }
3144 
3145 
3146 /***********************************************************************//**
3147  * @brief Compute fiducial offset to force (x,y) = (0,0) at (phi0,theta0)
3148  *
3149  * @param[in] phi0 Fiducial longitude
3150  * @param[in] theta0 Fiducial latitude
3151  *
3152  * Code adapted from prj.c::prjoff().
3153  ***************************************************************************/
3154 void GWcs::prj_off(const double& phi0, const double& theta0) const
3155 {
3156  // Initialise fiducial offsets
3157  m_x0 = 0.0;
3158  m_y0 = 0.0;
3159 
3160  // Set both to the projection-specific default if either undefined
3161  if (undefined(m_phi0) || undefined(m_theta0)) {
3162  m_phi0 = phi0;
3163  m_theta0 = theta0;
3164  }
3165 
3166  // ... otherwise compute (x,y) for (phi0,theta0)
3167  else {
3168  // Get (x,y) at (phi0,theta0) from projection
3169  int stat;
3170  double x0;
3171  double y0;
3172  prj_s2x(1, 1, 1, 1, &m_phi0, &m_theta0, &x0, &y0, &stat);
3173  m_x0 = x0;
3174  m_y0 = y0;
3175  }
3176 
3177  // Return
3178  return;
3179 }
3180 
3181 
3182 /***********************************************************************//**
3183  * @brief Performs bounds checking on native spherical coordinates
3184  *
3185  * @param[in] tol Tolerance for the bounds check [deg]
3186  * @param[in] nphi Number of bins in phi.
3187  * @param[in] ntheta Number of bins in theta.
3188  * @param[in] spt Step size.
3189  * @param[in,out] phi Pointer to phi array [deg]
3190  * @param[in,out] theta Pointer to theta array [deg]
3191  * @param[in,out] stat Pointer to status array.
3192  *
3193  * Performs bounds checking on native spherical coordinates. As returned by
3194  * the deprojection (x2s) routines, native longitude is expected to lie in
3195  * the closed interval [-180,180], with latitude in [-90,90].
3196  *
3197  * A tolerance may be specified to provide a small allowance for numerical
3198  * imprecision. Values that lie outside the allowed range by not more than
3199  * the specified tolerance will be adjusted back into range.
3200  *
3201  * Code adapted from prj.c::prjbchk().
3202  ***************************************************************************/
3203 int GWcs::prj_bchk(const double& tol,
3204  const int& nphi,
3205  const int& ntheta,
3206  const int& spt,
3207  double* phi,
3208  double* theta,
3209  int* stat) const
3210 {
3211  // Initialize result
3212  int status = 0;
3213 
3214  // Continue only if bounds checking is requested
3215  if (m_bounds) {
3216 
3217  // Initialise array pointers
3218  register double* phip = phi;
3219  register double* thetap = theta;
3220  register int* statp = stat;
3221 
3222  // Loop over all theta bins
3223  for (int itheta = 0; itheta < ntheta; ++itheta) {
3224 
3225  // Loop over all phi bins
3226  for (int iphi = 0; iphi < nphi; ++iphi, phip += spt, thetap += spt, statp++) {
3227 
3228  // Skip values already marked as illegal
3229  if (*statp == 0) {
3230 
3231  // Check lower phi boundary
3232  if (*phip < -180.0) {
3233  if (*phip < -180.0-tol) {
3234  *statp = 1;
3235  status = 1;
3236  }
3237  else {
3238  *phip = -180.0;
3239  }
3240  }
3241 
3242  // Check upper phi boundary
3243  else if (180.0 < *phip) {
3244  if (180.0+tol < *phip) {
3245  *statp = 1;
3246  status = 1;
3247  }
3248  else {
3249  *phip = 180.0;
3250  }
3251  }
3252 
3253  // Check lower theta boundary
3254  if (*thetap < -90.0) {
3255  if (*thetap < -90.0-tol) {
3256  *statp = 1;
3257  status = 1;
3258  }
3259  else {
3260  *thetap = -90.0;
3261  }
3262  }
3263 
3264  // Check upper theta boundary
3265  else if (90.0 < *thetap) {
3266  if (90.0+tol < *thetap) {
3267  *statp = 1;
3268  status = 1;
3269  }
3270  else {
3271  *thetap = 90.0;
3272  }
3273  }
3274 
3275  } // endif: values were marked as valid
3276 
3277  } // endfor: looped over phi
3278 
3279  } // endfor: looped over theta
3280 
3281  } // endif: bounds checking was requested
3282 
3283  // Return status
3284  return status;
3285 }
#define G_WCS_SET_CTYPE
Definition: GWcs.cpp:45
bool m_linset
Linear transformation is set.
Definition: GWcs.hpp:190
void wcs_p2s(int ncoord, int nelem, const double *pixcrd, double *imgcrd, double *phi, double *theta, double *world, int *stat) const
Pixel-to-world transformation.
Definition: GWcs.cpp:1491
void init_members(void)
Initialise class members.
Definition: GWcs.cpp:949
GSkyPixel m_last_pix2dir_pix
Last pixel for pix2dir.
Definition: GWcs.hpp:227
void cel_ini(void) const
Initialise celestial transformation parameters.
Definition: GWcs.cpp:1826
std::vector< double > m_imgpix
Image to pixel transformation matrix.
Definition: GWcs.hpp:196
bool has_card(const int &cardno) const
Check existence of header card.
Definition: GFitsHDU.hpp:233
void set_members(const std::string &coords, const double &crval1, const double &crval2, const double &crpix1, const double &crpix2, const double &cdelt1, const double &cdelt2)
Set World Coordinate System parameters.
Definition: GWcs.cpp:1077
virtual GSkyDir pix2dir(const GSkyPixel &pixel) const
Returns sky direction of sky map pixel.
Definition: GWcs.cpp:682
double dec_deg(void) const
Returns Declination in degrees.
Definition: GSkyDir.hpp:250
void prj_off(const double &phi0, const double &theta0) const
Compute fiducial offset to force (x,y) = (0,0) at (phi0,theta0)
Definition: GWcs.cpp:3154
int m_naxis
Number of axes.
Definition: GWcs.hpp:172
void prj_ini(void) const
Initialise projection parameters.
Definition: GWcs.cpp:3124
void x(const double &x)
Set x value of sky map pixel.
Definition: GSkyPixel.hpp:204
std::vector< double > m_crpix
CRPIXja keyvalues for each pixel axis.
Definition: GWcs.hpp:192
const double pi
Definition: GMath.hpp:35
void set(const std::string &coords, const double &crval1, const double &crval2, const double &crpix1, const double &crpix2, const double &cdelt1, const double &cdelt2)
Set World Coordinate System parameters.
Definition: GWcs.cpp:849
std::vector< std::string > m_ctype
CTYPEia keyvalues for each coord axis.
Definition: GWcs.hpp:175
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1163
void lin_p2x(int ncoord, int nelem, const double *pixcrd, double *imgcrd) const
Pixel-to-world linear transformation.
Definition: GWcs.cpp:2832
std::vector< double > m_piximg
Pixel to image transformation matrix.
Definition: GWcs.hpp:195
Abstract FITS extension base class.
Definition: GFitsHDU.hpp:51
virtual void prj_set(void) const =0
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1100
bool undefined(const double &value) const
Check if value is undefined.
Definition: GWcs.hpp:262
int m_lng
Longitude axis.
Definition: GWcs.hpp:185
bool m_celset
Celestial transformation is set.
Definition: GWcs.hpp:199
int m_latpreq
LATPOLEa requirement.
Definition: GWcs.hpp:207
#define G_CRPIX
Definition: GWcs.cpp:43
#define G_LIN_MATINV
Definition: GWcs.cpp:51
void xy(const double &x, const double &y)
Set x and y value of sky map pixel.
Definition: GSkyPixel.hpp:238
void wcs_s2p(int ncoord, int nelem, const double *world, double *phi, double *theta, double *imgcrd, double *pixcrd, int *stat) const
World-to-pixel transformation.
Definition: GWcs.cpp:1557
int m_spec
Spectral axis.
Definition: GWcs.hpp:187
virtual GSkyProjection & operator=(const GSkyProjection &proj)
Assignment operator.
double acos(const double &arg)
Computes acos by avoiding NaN due to rounding errors.
Definition: GMath.cpp:69
void cel_x2s(int nx, int ny, int sxy, int sll, const double *x, const double *y, double *phi, double *theta, double *lng, double *lat, int *stat) const
Pixel-to-world celestial transformation.
Definition: GWcs.cpp:2151
void spc_ini(void)
Initialise spectral transformation parameters.
Definition: GWcs.cpp:2694
double crval(const int &inx) const
Return value of reference pixel.
Definition: GWcs.cpp:877
double m_ref[4]
Celestial coordinates of fiducial.
Definition: GWcs.hpp:203
void unset_lock(const int &lock_id=0) const
Releases a previously set OpenMP lock.
Definition: GWcs.hpp:305
void clear(void)
Clear instance.
Definition: GSkyPixel.cpp:259
virtual double solidangle(const GSkyPixel &pixel) const
Returns solid angle of pixel in units of steradians.
Definition: GWcs.cpp:477
Gammalib tools definition.
void y(const double &y)
Set y value of sky map pixel.
Definition: GSkyPixel.hpp:221
GSkyPixel m_last_dir2pix_pix
Last pixel for dir2pix.
Definition: GWcs.hpp:228
std::vector< double > m_cdelt
CDELTia keyvalues for each coord axis.
Definition: GWcs.hpp:194
void sph_s2x(int nlng, int nlat, int sll, int spt, const double *lng, const double *lat, double *phi, double *theta) const
Rotation in the pixel-to-world direction.
Definition: GWcs.cpp:2489
double m_equinox
EQUINOX keyvalue.
Definition: GWcs.hpp:182
void free_members(void)
Delete class members.
Definition: GWcs.cpp:1053
void cel_set(void) const
Setup of celestial transformation.
Definition: GWcs.cpp:1881
std::string strip_whitespace(const std::string &arg)
Strip leading and trailing whitespace from string.
Definition: GTools.cpp:73
virtual std::string code(void) const =0
void wcs_ini(int naxis)
Initialise World Coordinate System.
Definition: GWcs.cpp:1157
#define G_CRVAL
Definition: GWcs.cpp:42
double real(const std::string &keyname) const
Return card value as double precision.
Definition: GFitsHDU.hpp:423
std::vector< double > m_crval
CRVALia keyvalues for each coord axis.
Definition: GWcs.hpp:173
Sky map pixel class.
Definition: GSkyPixel.hpp:74
double m_latpole
LATPOLEa keyvalue.
Definition: GWcs.hpp:178
const double deg2rad
Definition: GMath.hpp:43
bool m_unity
Signals unity PC matrix.
Definition: GWcs.hpp:191
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1268
void init_lock(const int &lock_id=0) const
Initializes an OpenMP lock with a specific name.
Definition: GWcs.hpp:274
double crpix(const int &inx) const
Return reference pixel.
Definition: GWcs.cpp:900
#define G_WCS_P2S
Definition: GWcs.cpp:46
std::vector< double > m_crota
CROTAia keyvalues for each coord axis.
Definition: GWcs.hpp:184
void wcs_set(void) const
Setup of World Coordinate System.
Definition: GWcs.cpp:1249
virtual void prj_s2x(int nphi, int ntheta, int spt, int sxy, const double *phi, const double *theta, double *x, double *y, int *stat) const =0
double m_lonpole
LONPOLEa keyvalue.
Definition: GWcs.hpp:177
double acosd(const double &value)
Compute arc cosine in degrees.
Definition: GMath.cpp:218
double m_restfrq
RESTFRQa keyvalue.
Definition: GWcs.hpp:179
GVector tan(const GVector &vector)
Computes tangens of vector elements.
Definition: GVector.cpp:1289
std::vector< std::string > m_ctype_c
CTYPEia comments for each coord axis.
Definition: GWcs.hpp:176
GVector cross(const GVector &a, const GVector &b)
Vector cross product.
Definition: GVector.cpp:762
void celvector(const GVector &vector)
Set sky direction from 3D vector in celestial coordinates.
Definition: GSkyDir.cpp:291
#define G_READ
Definition: GWcs.cpp:39
void sincosd(const double &angle, double *s, double *c)
Compute sine and cosine of angle in degrees.
Definition: GMath.cpp:342
double cdelt(const int &inx) const
Return pixel size.
Definition: GWcs.cpp:923
virtual GWcs & operator=(const GWcs &wcs)
Assignment operator.
Definition: GWcs.cpp:180
void wcs_set_ctype(void) const
Set CTYPEa keywords.
Definition: GWcs.cpp:1378
virtual void write(GFitsHDU &hdu) const
Write WCS definition into FITS HDU header.
Definition: GWcs.cpp:332
virtual GSkyPixel dir2pix(const GSkyDir &dir) const
Returns sky map pixel of sky direction.
Definition: GWcs.cpp:766
int m_coordsys
0=EQU, 1=GAL
std::vector< double > m_cd
CDi_ja linear transformation matrix.
Definition: GWcs.hpp:183
double ra_deg(void) const
Returns Right Ascension in degrees.
Definition: GSkyDir.hpp:223
virtual void read(const GFitsHDU &hdu)
Read WCS definition from FITS header.
Definition: GWcs.cpp:220
double m_pv[PVN]
Projection parameters.
Definition: GWcs.hpp:213
void sph_x2s(int nphi, int ntheta, int spt, int sll, const double *phi, const double *theta, double *lng, double *lat) const
Rotation in the pixel-to-world direction.
Definition: GWcs.cpp:2253
virtual bool compare(const GSkyProjection &proj) const
Compares sky projection.
Definition: GWcs.cpp:1116
#define G_CDELT
Definition: GWcs.cpp:44
double m_x0
Fiducial x offset.
Definition: GWcs.hpp:215
GChatter
Definition: GTypemaps.hpp:33
GSkyDir m_last_pix2dir_dir
Last sky direction for pix2dir.
Definition: GWcs.hpp:225
double m_restwav
RESTWAVa keyvalue.
Definition: GWcs.hpp:180
std::string wcs_print(const GChatter &chatter=NORMAL) const
Print WCS information.
Definition: GWcs.cpp:1602
GSkyDir m_last_dir2pix_dir
Last sky direction for dir2pix.
Definition: GWcs.hpp:226
double m_theta0
Native zenith angle of fiducial point.
Definition: GWcs.hpp:202
void radec_deg(const double &ra, const double &dec)
Set equatorial sky direction (degrees)
Definition: GSkyDir.cpp:202
void lin_ini(int naxis)
Initialise linear transformation parameters.
Definition: GWcs.cpp:2712
std::string wcs_print_value(const double &value) const
Helper function for value printing.
Definition: GWcs.cpp:1797
void wcs_set_radesys(void) const
Set radesys and equinox members.
Definition: GWcs.cpp:1315
virtual std::string coordsys(void) const
Returns coordinate system.
bool m_isolat
True if |latitude| is preserved.
Definition: GWcs.hpp:208
virtual ~GWcs(void)
Destructor.
Definition: GWcs.cpp:158
double m_phi0
Native azimuth angle of fiducial point.
Definition: GWcs.hpp:201
std::vector< double > m_w
Intermediate values.
Definition: GWcs.hpp:217
bool m_has_pix2dir_cache
Has valid pix2dir cache value.
Definition: GWcs.hpp:223
bool m_offset
Force (x,y) = (0,0) at (phi_0,theta_0)
Definition: GWcs.hpp:200
std::string m_radesys
RADESYS keyvalue.
Definition: GWcs.hpp:181
Abstract world coordinate system base class.
Definition: GWcs.hpp:51
int prj_bchk(const double &tol, const int &nphi, const int &ntheta, const int &spt, double *phi, double *theta, int *stat) const
Performs bounds checking on native spherical coordinates.
Definition: GWcs.cpp:3203
void clear(void)
Clear sky direction.
Definition: GSkyDir.cpp:142
GWcs(void)
Void constructor.
Definition: GWcs.cpp:80
double m_y0
Fiducial y offset.
Definition: GWcs.hpp:216
double atan2d(const double &y, const double &x)
Compute arc tangens in degrees.
Definition: GMath.cpp:308
void copy_members(const GWcs &wcs)
Copy class members.
Definition: GWcs.cpp:979
double cosd(const double &angle)
Compute cosine of angle in degrees.
Definition: GMath.cpp:134
void set_lock(const int &lock_id=0) const
Sets an OpenMP lock with a specific name.
Definition: GWcs.hpp:291
void lin_matinv(const std::vector< double > &mat, std::vector< double > &inv) const
Invert linear transformation matrix.
Definition: GWcs.cpp:2979
void lb_deg(const double &l, const double &b)
Set galactic sky direction (degrees)
Definition: GSkyDir.cpp:256
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1226
void cel_s2x(int nlng, int nlat, int sll, int sxy, const double *lng, const double *lat, double *phi, double *theta, double *x, double *y, int *stat) const
World-to-pixel celestial transformation.
Definition: GWcs.cpp:2197
std::string string(const std::string &keyname) const
Return card value as string.
Definition: GFitsHDU.hpp:410
Exception handler interface definition.
std::vector< double > m_pc
PCi_ja linear transformation matrix.
Definition: GWcs.hpp:193
double asind(const double &value)
Compute arc sine in degrees.
Definition: GMath.cpp:250
bool m_prjset
Projection is set.
Definition: GWcs.hpp:211
double dist(const GSkyDir &dir) const
Compute angular distance between sky directions in radians.
Definition: GSkyDir.hpp:265
virtual void clear(void)=0
Clear object.
const double twopi
Definition: GMath.hpp:36
Vector class.
Definition: GVector.hpp:46
double l_deg(void) const
Return galactic longitude in degrees.
Definition: GSkyDir.hpp:169
double b_deg(void) const
Returns galactic latitude in degrees.
Definition: GSkyDir.hpp:196
Abstract sky projection base class.
void lin_set(void) const
Initialise linear transformation parameters.
Definition: GWcs.cpp:2762
static const double UNDEFINED
Definition: GWcs.hpp:95
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
int m_lat
Latitude axis.
Definition: GWcs.hpp:186
bool m_wcsset
WCS information is set.
Definition: GWcs.hpp:171
double m_euler[5]
Euler angles and functions thereof.
Definition: GWcs.hpp:206
static const int PVN
Definition: GWcs.hpp:94
GFitsHeaderCard & card(const int &cardno)
Return header card.
Definition: GFitsHDU.hpp:259
std::vector< std::string > m_cunit
CUNITia keyvalues for each coord axis.
Definition: GWcs.hpp:174
Sky direction class.
Definition: GSkyDir.hpp:62
virtual void prj_x2s(int nx, int ny, int sxy, int spt, const double *x, const double *y, double *phi, double *theta, int *stat) const =0
GVector atan(const GVector &vector)
Computes arctan of vector elements.
Definition: GVector.cpp:1058
double sind(const double &angle)
Compute sine of angle in degrees.
Definition: GMath.cpp:163
#define G_WCS_S2P
Definition: GWcs.cpp:48
virtual std::string name(void) const =0
Abstract world coordinate system base class definition.
double m_r0
Radius of the generating sphere.
Definition: GWcs.hpp:212
Mathematical function definitions.
void lin_x2p(int ncoord, int nelem, const double *imgcrd, double *pixcrd) const
World-to-pixel linear transformation.
Definition: GWcs.cpp:2908
bool m_bounds
Enable strict bounds checking.
Definition: GWcs.hpp:214
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413
bool m_has_dir2pix_cache
Has valid dir2pix cache value.
Definition: GWcs.hpp:224
#define G_CEL_SET
Definition: GWcs.cpp:50