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