GammaLib 2.0.0
Loading...
Searching...
No Matches
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 __________________________________________________________ */
68const double GWcs::UNDEFINED = 987654321.0e99;
69
70
71/*==========================================================================
72 = =
73 = Constructors/destructors =
74 = =
75 ==========================================================================*/
76
77/***********************************************************************//**
78 * @brief Void constructor
79 ***************************************************************************/
80GWcs::GWcs(void) : GSkyProjection()
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 ***************************************************************************/
103GWcs::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 ***************************************************************************/
124GWcs::GWcs(const GFitsHDU& hdu) : GSkyProjection()
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 ***************************************************************************/
142GWcs::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 ***************************************************************************/
158GWcs::~GWcs(void)
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 ***************************************************************************/
180GWcs& GWcs::operator=(const GWcs& wcs)
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 ***************************************************************************/
220void GWcs::read(const GFitsHDU& hdu)
221{
222 // Free class members (base and derived classes, derived class first)
223 free_members();
224 this->GSkyProjection::free_members();
225
226 // Initialise members
227 this->GSkyProjection::init_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 ***************************************************************************/
342void 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 ***************************************************************************/
487double 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) -
670 gammalib::twopi;
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 ***************************************************************************/
692GSkyDir 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 ***************************************************************************/
776GSkyPixel 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 ***************************************************************************/
859void 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 ***************************************************************************/
887double 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 ***************************************************************************/
911double 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 ***************************************************************************/
935double 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 ***************************************************************************/
962void GWcs::init_members(void)
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;
973 m_last_pix2dir_dir.clear();
974 m_last_dir2pix_dir.clear();
975 m_last_pix2dir_pix.clear();
976 m_last_dir2pix_pix.clear();
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 ***************************************************************************/
992void 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
1051 m_has_pix2dir_cache = wcs.m_has_pix2dir_cache;
1052 m_has_dir2pix_cache = wcs.m_has_dir2pix_cache;
1053 m_last_pix2dir_dir = wcs.m_last_pix2dir_dir;
1054 m_last_dir2pix_dir = wcs.m_last_dir2pix_dir;
1055 m_last_pix2dir_pix = wcs.m_last_pix2dir_pix;
1056 m_last_dir2pix_pix = wcs.m_last_dir2pix_pix;
1057
1058 // Return
1059 return;
1060}
1061
1062
1063/***********************************************************************//**
1064 * @brief Delete class members
1065 ***************************************************************************/
1066void GWcs::free_members(void)
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 ***************************************************************************/
1090void 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 ***************************************************************************/
1129bool 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 ***************************************************************************/
1170void 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 ***************************************************************************/
1262void 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 ***************************************************************************/
1328void 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 ***************************************************************************/
1391void 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.";
1398 throw GException::invalid_value(G_WCS_SET_CTYPE, msg);
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).";
1430 throw GException::invalid_value(G_WCS_SET_CTYPE, msg);
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).";
1467 throw GException::invalid_value(G_WCS_SET_CTYPE, msg);
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 ***************************************************************************/
1514void 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.";
1533 throw GException::invalid_argument(G_WCS_P2S, msg);
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.";
1540 throw GException::invalid_argument(G_WCS_P2S, msg);
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 ***************************************************************************/
1589void 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.";
1608 throw GException::invalid_argument(G_WCS_S2P, msg);
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.";
1615 throw GException::invalid_argument(G_WCS_S2P, msg);
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 ***************************************************************************/
1642std::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 ***************************************************************************/
1837std::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 ***************************************************************************/
1866void 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 ***************************************************************************/
1921void 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.";
2008 throw GException::runtime_error(G_CEL_SET, msg);
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.";
2038 throw GException::runtime_error(G_CEL_SET, msg);
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.";
2117 throw GException::runtime_error(G_CEL_SET, msg);
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.";
2161 throw GException::runtime_error(G_CEL_SET, msg);
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 ***************************************************************************/
2194void 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 ***************************************************************************/
2247void 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 ***************************************************************************/
2310void 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 ***************************************************************************/
2551void 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 ***************************************************************************/
2761void 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 ***************************************************************************/
2779void 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 ***************************************************************************/
2829void 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
2874 lin_matinv(m_piximg, m_imgpix);
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 ***************************************************************************/
2899void 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 ***************************************************************************/
2978void 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 ***************************************************************************/
3052void 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.";
3089 throw GException::runtime_error(G_LIN_MATINV, msg);
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 ***************************************************************************/
3202void 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 ***************************************************************************/
3232void 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 ***************************************************************************/
3281int 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}
G_READ
#define G_READ
Definition GApplicationPars.cpp:59
GException.hpp
Exception handler interface definition.
GMath.hpp
Mathematical function definitions.
GTools.hpp
Gammalib tools definition.
GChatter
GChatter
Definition GTypemaps.hpp:33
EXPLICIT
@ EXPLICIT
Definition GTypemaps.hpp:37
SILENT
@ SILENT
Definition GTypemaps.hpp:34
sqrt
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition GVector.cpp:1358
cross
GVector cross(const GVector &a, const GVector &b)
Vector cross product.
Definition GVector.cpp:793
G_LIN_MATINV
#define G_LIN_MATINV
Definition GWcs.cpp:51
G_WCS_S2P
#define G_WCS_S2P
Definition GWcs.cpp:48
G_CRVAL
#define G_CRVAL
Definition GWcs.cpp:42
G_CEL_SET
#define G_CEL_SET
Definition GWcs.cpp:50
G_CRPIX
#define G_CRPIX
Definition GWcs.cpp:43
G_CDELT
#define G_CDELT
Definition GWcs.cpp:44
G_WCS_SET_CTYPE
#define G_WCS_SET_CTYPE
Definition GWcs.cpp:45
G_WCS_P2S
#define G_WCS_P2S
Definition GWcs.cpp:46
GWcs.hpp
Abstract world coordinate system base class definition.
GException::invalid_argument
Definition GException.hpp:88
GException::invalid_value
Definition GException.hpp:81
GException::out_of_range
Definition GException.hpp:105
GException::runtime_error
Definition GException.hpp:126
GFitsHDU
Abstract FITS extension base class.
Definition GFitsHDU.hpp:51
GFitsHDU::has_card
bool has_card(const int &cardno) const
Check existence of header card.
Definition GFitsHDU.hpp:233
GFitsHDU::real
double real(const std::string &keyname) const
Return card value as double precision.
Definition GFitsHDU.hpp:423
GFitsHDU::string
std::string string(const std::string &keyname) const
Return card value as string.
Definition GFitsHDU.hpp:410
GFitsHDU::card
GFitsHeaderCard & card(const int &cardno)
Return header card.
Definition GFitsHDU.hpp:259
GSkyDir
Sky direction class.
Definition GSkyDir.hpp:62
GSkyDir::lb_deg
void lb_deg(const double &l, const double &b)
Set galactic sky direction (degrees)
Definition GSkyDir.cpp:278
GSkyDir::dec_deg
double dec_deg(void) const
Returns Declination in degrees.
Definition GSkyDir.hpp:256
GSkyDir::ra_deg
double ra_deg(void) const
Returns Right Ascension in degrees.
Definition GSkyDir.hpp:229
GSkyDir::radec_deg
void radec_deg(const double &ra, const double &dec)
Set equatorial sky direction (degrees)
Definition GSkyDir.cpp:224
GSkyDir::b_deg
double b_deg(void) const
Returns galactic latitude in degrees.
Definition GSkyDir.hpp:202
GSkyDir::clear
void clear(void)
Clear sky direction.
Definition GSkyDir.cpp:164
GSkyDir::dist
double dist(const GSkyDir &dir) const
Compute angular distance between sky directions in radians.
Definition GSkyDir.hpp:271
GSkyDir::celvector
void celvector(const GVector &vector)
Set sky direction from 3D vector in celestial coordinates.
Definition GSkyDir.cpp:313
GSkyDir::l_deg
double l_deg(void) const
Return galactic longitude in degrees.
Definition GSkyDir.hpp:175
GSkyPixel
Sky map pixel class.
Definition GSkyPixel.hpp:74
GSkyPixel::clear
void clear(void)
Clear instance.
Definition GSkyPixel.cpp:259
GSkyPixel::y
void y(const double &y)
Set y value of sky map pixel.
Definition GSkyPixel.hpp:221
GSkyPixel::x
void x(const double &x)
Set x value of sky map pixel.
Definition GSkyPixel.hpp:204
GSkyPixel::xy
void xy(const double &x, const double &y)
Set x and y value of sky map pixel.
Definition GSkyPixel.hpp:238
GSkyProjection
Abstract sky projection base class.
Definition GSkyProjection.hpp:47
GSkyProjection::free_members
void free_members(void)
Delete class members.
Definition GSkyProjection.cpp:242
GSkyProjection::coordsys
virtual std::string coordsys(void) const
Returns coordinate system.
Definition GSkyProjection.cpp:147
GSkyProjection::m_coordsys
int m_coordsys
0=CEL, 1=GAL
Definition GSkyProjection.hpp:88
GSkyProjection::init_members
void init_members(void)
Initialise class members.
Definition GSkyProjection.cpp:214
GSkyProjection::operator=
virtual GSkyProjection & operator=(const GSkyProjection &proj)
Assignment operator.
Definition GSkyProjection.cpp:111
GVector
Vector class.
Definition GVector.hpp:46
GWcs
Abstract world coordinate system base class.
Definition GWcs.hpp:51
GWcs::cel_set
void cel_set(void) const
Setup of celestial transformation.
Definition GWcs.cpp:1921
GWcs::cdelt
double cdelt(const int &inx) const
Return pixel size.
Definition GWcs.cpp:935
GWcs::free_members
void free_members(void)
Delete class members.
Definition GWcs.cpp:1066
GWcs::wcs_ini
void wcs_ini(int naxis)
Initialise World Coordinate System.
Definition GWcs.cpp:1170
GWcs::copy_members
void copy_members(const GWcs &wcs)
Copy class members.
Definition GWcs.cpp:992
GWcs::read
virtual void read(const GFitsHDU &hdu)
Read WCS definition from FITS header.
Definition GWcs.cpp:220
GWcs::sph_x2s
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
GWcs::lin_p2x
void lin_p2x(int ncoord, int nelem, const double *pixcrd, double *imgcrd) const
Pixel-to-world linear transformation.
Definition GWcs.cpp:2899
GWcs::m_radesys
std::string m_radesys
RADESYS keyvalue.
Definition GWcs.hpp:181
GWcs::PVN
static const int PVN
Definition GWcs.hpp:94
GWcs::m_lat
int m_lat
Latitude axis.
Definition GWcs.hpp:186
GWcs::pix2dir
virtual GSkyDir pix2dir(const GSkyPixel &pixel) const
Returns sky direction of sky map pixel.
Definition GWcs.cpp:692
GWcs::~GWcs
virtual ~GWcs(void)
Destructor.
Definition GWcs.cpp:158
GWcs::crpix
double crpix(const int &inx) const
Return reference pixel.
Definition GWcs.cpp:911
GWcs::prj_set
virtual void prj_set(void) const =0
GWcs::m_restfrq
double m_restfrq
RESTFRQa keyvalue.
Definition GWcs.hpp:179
GWcs::prj_s2x
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
GWcs::m_crval
std::vector< double > m_crval
CRVALia keyvalues for each coord axis.
Definition GWcs.hpp:173
GWcs::m_euler
double m_euler[5]
Euler angles and functions thereof.
Definition GWcs.hpp:206
GWcs::spc_ini
void spc_ini(void)
Initialise spectral transformation parameters.
Definition GWcs.cpp:2761
GWcs::m_last_dir2pix_dir
GSkyDir m_last_dir2pix_dir
Last sky direction for dir2pix.
Definition GWcs.hpp:226
GWcs::sph_s2x
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
GWcs::set
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
GWcs::prj_x2s
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
GWcs::m_unity
bool m_unity
Signals unity PC matrix.
Definition GWcs.hpp:191
GWcs::lin_ini
void lin_ini(int naxis)
Initialise linear transformation parameters.
Definition GWcs.cpp:2779
GWcs::lin_x2p
void lin_x2p(int ncoord, int nelem, const double *imgcrd, double *pixcrd) const
World-to-pixel linear transformation.
Definition GWcs.cpp:2978
GWcs::m_has_pix2dir_cache
bool m_has_pix2dir_cache
Has valid pix2dir cache value.
Definition GWcs.hpp:223
GWcs::write
virtual void write(GFitsHDU &hdu) const
Write WCS definition into FITS HDU header.
Definition GWcs.cpp:342
GWcs::m_piximg
std::vector< double > m_piximg
Pixel to image transformation matrix.
Definition GWcs.hpp:195
GWcs::prj_bchk
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
GWcs::code
virtual std::string code(void) const =0
GWcs::m_has_dir2pix_cache
bool m_has_dir2pix_cache
Has valid dir2pix cache value.
Definition GWcs.hpp:224
GWcs::prj_off
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
GWcs::m_wcsset
bool m_wcsset
WCS information is set.
Definition GWcs.hpp:171
GWcs::wcs_print_value
std::string wcs_print_value(const double &value) const
Helper function for value printing.
Definition GWcs.cpp:1837
GWcs::UNDEFINED
static const double UNDEFINED
Definition GWcs.hpp:95
GWcs::compare
virtual bool compare(const GSkyProjection &proj) const
Compares sky projection.
Definition GWcs.cpp:1129
GWcs::m_linset
bool m_linset
Linear transformation is set.
Definition GWcs.hpp:190
GWcs::m_imgpix
std::vector< double > m_imgpix
Image to pixel transformation matrix.
Definition GWcs.hpp:196
GWcs::wcs_set_ctype
void wcs_set_ctype(void) const
Set CTYPEa keywords.
Definition GWcs.cpp:1391
GWcs::prj_ini
void prj_ini(void) const
Initialise projection parameters.
Definition GWcs.cpp:3202
GWcs::cel_x2s
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
GWcs::solidangle
virtual double solidangle(const GSkyPixel &pixel) const
Returns solid angle of pixel in units of steradians.
Definition GWcs.cpp:487
GWcs::m_naxis
int m_naxis
Number of axes.
Definition GWcs.hpp:172
GWcs::operator=
virtual GWcs & operator=(const GWcs &wcs)
Assignment operator.
Definition GWcs.cpp:180
GWcs::m_lng
int m_lng
Longitude axis.
Definition GWcs.hpp:185
GWcs::m_latpole
double m_latpole
LATPOLEa keyvalue.
Definition GWcs.hpp:178
GWcs::m_lonpole
double m_lonpole
LONPOLEa keyvalue.
Definition GWcs.hpp:177
GWcs::wcs_print
std::string wcs_print(const GChatter &chatter=NORMAL) const
Print WCS information.
Definition GWcs.cpp:1642
GWcs::wcs_s2p
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
GWcs::init_members
void init_members(void)
Initialise class members.
Definition GWcs.cpp:962
GWcs::lin_set
void lin_set(void) const
Initialise linear transformation parameters.
Definition GWcs.cpp:2829
GWcs::m_crpix
std::vector< double > m_crpix
CRPIXja keyvalues for each pixel axis.
Definition GWcs.hpp:192
GWcs::m_pv
double m_pv[PVN]
Projection parameters.
Definition GWcs.hpp:213
GWcs::GWcs
GWcs(void)
Void constructor.
Definition GWcs.cpp:80
GWcs::m_ctype_c
std::vector< std::string > m_ctype_c
CTYPEia comments for each coord axis.
Definition GWcs.hpp:176
GWcs::m_phi0
double m_phi0
Native azimuth angle of fiducial point.
Definition GWcs.hpp:201
GWcs::m_cunit
std::vector< std::string > m_cunit
CUNITia keyvalues for each coord axis.
Definition GWcs.hpp:174
GWcs::cel_s2x
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
GWcs::unset_lock
void unset_lock(const int &lock_id=0) const
Releases a previously set OpenMP lock.
Definition GWcs.hpp:305
GWcs::m_cdelt
std::vector< double > m_cdelt
CDELTia keyvalues for each coord axis.
Definition GWcs.hpp:194
GWcs::set_members
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
GWcs::wcs_set
void wcs_set(void) const
Setup of World Coordinate System.
Definition GWcs.cpp:1262
GWcs::m_crota
std::vector< double > m_crota
CROTAia keyvalues for each coord axis.
Definition GWcs.hpp:184
GWcs::cel_ini
void cel_ini(void) const
Initialise celestial transformation parameters.
Definition GWcs.cpp:1866
GWcs::lin_matinv
void lin_matinv(const std::vector< double > &mat, std::vector< double > &inv) const
Invert linear transformation matrix.
Definition GWcs.cpp:3052
GWcs::crval
double crval(const int &inx) const
Return value of reference pixel.
Definition GWcs.cpp:887
GWcs::m_latpreq
int m_latpreq
LATPOLEa requirement.
Definition GWcs.hpp:207
GWcs::m_bounds
bool m_bounds
Enable strict bounds checking.
Definition GWcs.hpp:214
GWcs::m_ref
double m_ref[4]
Celestial coordinates of fiducial.
Definition GWcs.hpp:203
GWcs::m_isolat
bool m_isolat
True if |latitude| is preserved.
Definition GWcs.hpp:208
GWcs::clear
virtual void clear(void)=0
Clear object.
GWcs::m_celset
bool m_celset
Celestial transformation is set.
Definition GWcs.hpp:199
GWcs::m_ctype
std::vector< std::string > m_ctype
CTYPEia keyvalues for each coord axis.
Definition GWcs.hpp:175
GWcs::m_last_pix2dir_dir
GSkyDir m_last_pix2dir_dir
Last sky direction for pix2dir.
Definition GWcs.hpp:225
GWcs::m_theta0
double m_theta0
Native zenith angle of fiducial point.
Definition GWcs.hpp:202
GWcs::m_restwav
double m_restwav
RESTWAVa keyvalue.
Definition GWcs.hpp:180
GWcs::m_equinox
double m_equinox
EQUINOX keyvalue.
Definition GWcs.hpp:182
GWcs::wcs_p2s
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
GWcs::name
virtual std::string name(void) const =0
GWcs::init_lock
void init_lock(const int &lock_id=0) const
Initializes an OpenMP lock with a specific name.
Definition GWcs.hpp:274
GWcs::m_pc
std::vector< double > m_pc
PCi_ja linear transformation matrix.
Definition GWcs.hpp:193
GWcs::m_last_dir2pix_pix
GSkyPixel m_last_dir2pix_pix
Last pixel for dir2pix.
Definition GWcs.hpp:228
GWcs::m_spec
int m_spec
Spectral axis.
Definition GWcs.hpp:187
GWcs::m_prjset
bool m_prjset
Projection is set.
Definition GWcs.hpp:211
GWcs::m_w
std::vector< double > m_w
Intermediate values.
Definition GWcs.hpp:217
GWcs::m_offset
bool m_offset
Force (x,y) = (0,0) at (phi_0,theta_0)
Definition GWcs.hpp:200
GWcs::dir2pix
virtual GSkyPixel dir2pix(const GSkyDir &dir) const
Returns sky map pixel of sky direction.
Definition GWcs.cpp:776
GWcs::set_lock
void set_lock(const int &lock_id=0) const
Sets an OpenMP lock with a specific name.
Definition GWcs.hpp:291
GWcs::wcs_set_radesys
void wcs_set_radesys(void) const
Set radesys and equinox members.
Definition GWcs.cpp:1328
GWcs::m_y0
double m_y0
Fiducial y offset.
Definition GWcs.hpp:216
GWcs::m_last_pix2dir_pix
GSkyPixel m_last_pix2dir_pix
Last pixel for pix2dir.
Definition GWcs.hpp:227
GWcs::m_x0
double m_x0
Fiducial x offset.
Definition GWcs.hpp:215
GWcs::m_r0
double m_r0
Radius of the generating sphere.
Definition GWcs.hpp:212
GWcs::m_cd
std::vector< double > m_cd
CDi_ja linear transformation matrix.
Definition GWcs.hpp:183
GWcs::undefined
bool undefined(const double &value) const
Check if value is undefined.
Definition GWcs.hpp:262
gammalib::parformat
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition GTools.cpp:1143
gammalib::sincosd
void sincosd(const double &angle, double *s, double *c)
Compute sine and cosine of angle in degrees.
Definition GMath.cpp:342
gammalib::str
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition GTools.cpp:489
gammalib::pi
const double pi
Definition GMath.hpp:35
gammalib::asind
double asind(const double &value)
Compute arc sine in degrees.
Definition GMath.cpp:250
gammalib::strip_whitespace
std::string strip_whitespace(const std::string &arg)
Strip leading and trailing whitespace from string.
Definition GTools.cpp:80
gammalib::deg2rad
const double deg2rad
Definition GMath.hpp:43
gammalib::cosd
double cosd(const double &angle)
Compute cosine of angle in degrees.
Definition GMath.cpp:134
gammalib::sind
double sind(const double &angle)
Compute sine of angle in degrees.
Definition GMath.cpp:163
gammalib::acos
double acos(const double &arg)
Computes acos by avoiding NaN due to rounding errors.
Definition GMath.cpp:69
gammalib::atan2d
double atan2d(const double &y, const double &x)
Compute arc tangens in degrees.
Definition GMath.cpp:308
gammalib::acosd
double acosd(const double &value)
Compute arc cosine in degrees.
Definition GMath.cpp:218
gammalib::twopi
const double twopi
Definition GMath.hpp:36