Johan Bregeon / gammalib

Commit c1344de7157c596b359efb628479dfafac3d475c

Authored by knodlseder
1 parent c8fc0afc

update for new wcslib code

Inline Side-by-side
Showing 7 changed files with 657 additions and 170 deletions
include/GWcs.hpp
  View file @c1344de
1 1 /***************************************************************************
2 2 * GWcs.hpp - World Coordinate System virtual base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010 by Jurgen Knodlseder *
  4 + * copyright (C) 2010-2011 by Jurgen Knodlseder *
5 5 * ----------------------------------------------------------------------- *
6 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 2 of the License, or *
10   - * (at your option) any later version. *
  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/>. *
11 19 * *
12 20 ***************************************************************************/
13 21 /**
... ... @@ -51,10 +59,10 @@ class GWcs {
51 59 public:
52 60 // Constructors and destructors
53 61 GWcs(void);
54   - explicit GWcs(const std::string& coords,
55   - const double& crval1, const double& crval2,
56   - const double& crpix1, const double& crpix2,
57   - const double& cdelt1, const double& cdelt2);
  62 + //explicit GWcs(const std::string& coords,
  63 + // const double& crval1, const double& crval2,
  64 + // const double& crpix1, const double& crpix2,
  65 + // const double& cdelt1, const double& cdelt2);
58 66 GWcs(const GWcs& wcs);
59 67 virtual ~GWcs(void);
60 68  
... ... @@ -93,15 +101,11 @@ protected:
93 101 void wcs_init(const double& theta0);
94 102 GVector wcs_getpole(const double& theta0);
95 103 GMatrix wcs_get_rot(void);
96   - void wcs_set(const std::string& coords,
97   - const double& crval1, const double& crval2,
98   - const double& crpix1, const double& crpix2,
99   - const double& cdelt1, const double& cdelt2);
100   - void wcs_read(const GFitsHDU* hdu);
101   - void wcs_write(GFitsHDU* hdu) const;
  104 + //void wcs_read(const GFitsHDU* hdu);
  105 + //void wcs_write(GFitsHDU* hdu) const;
102 106 std::string wcs_crval1(void) const;
103 107 std::string wcs_crval2(void) const;
104   - std::string wcs_dump(void) const;
  108 + //std::string wcs_dump(void) const;
105 109  
106 110 // Protected members (astr structure)
107 111 std::string m_type; //!< WCS type
... ... @@ -115,7 +119,7 @@ protected:
115 119 GVector m_pv2; //!< Projection parameters (up to 21)
116 120  
117 121 // Derived members
118   - double m_theta0; //!< Native latitude of the fiducial point
  122 + //double m_theta0; //!< Native latitude of the fiducial point
119 123 GVector m_refval; //!< Ordered value of reference point
120 124 GVector m_refpix; //!< Pixel of reference point (starting from 0)
121 125 GMatrix m_invcd; //!< Inverse of CD matrix
... ...
include/GWcsCAR.hpp
  View file @c1344de
1 1 /***************************************************************************
2 2 * GWcsCAR.hpp - Cartesian projection class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010 by Jurgen Knodlseder *
  4 + * copyright (C) 2010-2011 by Jurgen Knodlseder *
5 5 * ----------------------------------------------------------------------- *
6 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 2 of the License, or *
10   - * (at your option) any later version. *
  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/>. *
11 19 * *
12 20 ***************************************************************************/
13 21 /**
... ... @@ -20,7 +28,7 @@
20 28 #define GWCSCAR_HPP
21 29  
22 30 /* __ Includes ___________________________________________________________ */
23   -#include "GWcs.hpp"
  31 +#include "GWcslib.hpp"
24 32 #include "GFitsHDU.hpp"
25 33 #include "GSkyDir.hpp"
26 34 #include "GSkyPixel.hpp"
... ... @@ -31,10 +39,10 @@
31 39 *
32 40 * @brief GWcsCAR class interface defintion
33 41 ***************************************************************************/
34   -class GWcsCAR : public GWcs {
  42 +class GWcsCAR : public GWcslib {
35 43  
36 44 // I/O friends
37   - friend std::ostream& operator<< (std::ostream& os, const GWcsCAR& wcs);
  45 + //friend std::ostream& operator<< (std::ostream& os, const GWcsCAR& wcs);
38 46  
39 47 public:
40 48 // Constructors and destructors
... ... @@ -54,20 +62,33 @@ public:
54 62 // Implemented pure virtual base class methods
55 63 void clear(void);
56 64 GWcsCAR* clone(void) const;
57   - void read(const GFitsHDU* hdu);
58   - void write(GFitsHDU* hdu) const;
  65 + //void read(const GFitsHDU* hdu);
  66 + //void write(GFitsHDU* hdu) const;
59 67 std::string print(void) const;
60 68  
61 69 // Overloaded base class methods
62   - double omega(const GSkyPixel& pix) const;
  70 + //double omega(const GSkyPixel& pix) const;
63 71  
64 72 private:
65 73 // Private methods
66   - void init_members(void);
67   - void copy_members(const GWcsCAR& wcs);
68   - void free_members(void);
69   - void std2nat(GVector *coord) const;
70   - void nat2std(GVector *coord) const;
  74 + void init_members(void);
  75 + void copy_members(const GWcsCAR& wcs);
  76 + void free_members(void);
  77 + void std2nat(GVector *coord) const;
  78 + void nat2std(GVector *coord) const;
  79 +
  80 + // NEW VERSION
  81 + void prj_set(void) const;
  82 + void prj_x2s(int nx, int ny, int sxy, int spt,
  83 + const double* x, const double* y,
  84 + double* phi, double* theta, int* stat) const;
  85 + void prj_s2x(int nphi, int ntheta, int spt, int sxy,
  86 + const double* phi, const double* theta,
  87 + double* x, double* y, int* stat) const;
  88 +
  89 + // Protected members
  90 + mutable double m_w0;
  91 + mutable double m_w1;
71 92 };
72 93  
73 94 #endif /* GWCSCAR_HPP */
... ...
include/GWcsHPX.hpp
  View file @c1344de
1 1 /***************************************************************************
2 2 * GWcsHPX.hpp - Healpix projection class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010 by Jurgen Knodlseder *
  4 + * copyright (C) 2010-2011 by Jurgen Knodlseder *
5 5 * ----------------------------------------------------------------------- *
6 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 2 of the License, or *
10   - * (at your option) any later version. *
  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/>. *
11 19 * *
12 20 ***************************************************************************/
13 21 /**
... ... @@ -79,6 +87,15 @@ private:
79 87 int ang2pix_z_phi_nest(double z, double phi) const;
80 88 unsigned int isqrt(unsigned int arg) const;
81 89  
  90 + // NEW VERSION
  91 + void prj_set(void);
  92 + int prj_x2s(int nx, int ny, int sxy, int spt,
  93 + const double* x, const double* y,
  94 + double* phi, double* theta, int* stat) const;
  95 + int prj_s2x(int nphi, int ntheta, int spt, int sxy,
  96 + const double* phi, const double* theta,
  97 + double* x, double* y, int* stat) const;
  98 +
82 99 // Private data area
83 100 int m_nside; //!< Number of divisions of each base pixel (1-8192)
84 101 int m_npface; //!<
... ...
src/sky/GWcs.cpp
  View file @c1344de
... ... @@ -3,11 +3,18 @@
3 3 * ----------------------------------------------------------------------- *
4 4 * copyright (C) 2010-2011 by Jurgen Knodlseder *
5 5 * ----------------------------------------------------------------------- *
  6 + * This program is free software: you can redistribute it and/or modify *
  7 + * it under the terms of the GNU General Public License as published by *
  8 + * the Free Software Foundation, either version 3 of the License, or *
  9 + * (at your option) any later version. *
6 10 * *
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 2 of the License, or *
10   - * (at your option) any later version. *
  11 + * This program is distributed in the hope that it will be useful, *
  12 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  13 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  14 + * GNU General Public License for more details. *
  15 + * *
  16 + * You should have received a copy of the GNU General Public License *
  17 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
11 18 * *
12 19 ***************************************************************************/
13 20 /**
... ... @@ -20,6 +27,7 @@
20 27 #ifdef HAVE_CONFIG_H
21 28 #include <config.h>
22 29 #endif
  30 +#include <cmath>
23 31 #include "GException.hpp"
24 32 #include "GTools.hpp"
25 33 #include "GWcs.hpp"
... ... @@ -84,6 +92,7 @@ GWcs::GWcs(void)
84 92 * Construct standard WCS sky map from standard definition parameters. This
85 93 * method
86 94 ***************************************************************************/
  95 +/*
87 96 GWcs::GWcs(const std::string& coords,
88 97 const double& crval1, const double& crval2,
89 98 const double& crpix1, const double& crpix2,
... ... @@ -97,12 +106,12 @@ GWcs::GWcs(const std::string&amp; coords,
97 106 m_reverse = 0;
98 107  
99 108 // Set standard parameters
100   - wcs_set(coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2);
  109 + set_members(coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2);
101 110  
102 111 // Return
103 112 return;
104 113 }
105   -
  114 +*/
106 115  
107 116 /***********************************************************************//**
108 117 * @brief Copy constructor
... ... @@ -371,13 +380,13 @@ GSkyDir GWcs::xy2dir(const GSkyPixel&amp; pix) const
371 380 if (m_reverse)
372 381 xy = perm(xy, __permutation);
373 382  
374   - // Perform CAR map projection
  383 + // Perform map projection
375 384 GVector native = xy;
376 385 nat2std(&native);
377 386  
378 387 // Get sky direction for native coordinates
379 388 GSkyDir dir = wcs_native2dir(native);
380   -
  389 +
381 390 // Debug: Dump transformation steps
382 391 #if defined(G_XY2DIR_DEBUG)
383 392 std::cout << "xy2dir: pixel=" << offset << " xy=" << xy
... ... @@ -404,7 +413,7 @@ GSkyPixel GWcs::dir2xy(GSkyDir dir) const
404 413 // Get native coordinates for sky direction in radians
405 414 GVector native = wcs_dir2native(dir);
406 415  
407   - // Perform CAR map projection
  416 + // Perform map projection
408 417 GVector xy = native;
409 418 std2nat(&xy);
410 419  
... ... @@ -526,11 +535,22 @@ GSkyDir GWcs::wcs_native2dir(GVector native) const
526 535 ***************************************************************************/
527 536 void GWcs::wcs_init(const double& theta0)
528 537 {
  538 + // Initialise World Coordinate System. This has to be done here and
  539 + // not in wcs_set(...) as this method will be called by the derived
  540 + // class constructor. We're not allowed to call pure virtual functions
  541 + // in base class constructors !!!
  542 + //wcs_set();
  543 +
529 544 // Set native latitude of the fiducial point
530   - m_theta0 = theta0;
  545 + //m_theta0 = theta0;
  546 +
  547 + // Set pole
  548 + //m_npole[0] = (m_refval[1] < theta0) ? 180.0 : 0.0;
  549 + //m_npole[1] = 90.0;
531 550  
532 551 // Compute the coordinates of native pole for non-polar projection
533   - m_native_pole = wcs_getpole(m_theta0);
  552 + //m_native_pole = wcs_getpole(m_theta0);
  553 + m_native_pole = wcs_getpole(theta0);
534 554  
535 555 // Get matrix for rotation between sky and native coordinates
536 556 m_rot = wcs_get_rot();
... ... @@ -556,8 +576,9 @@ void GWcs::wcs_init(const double&amp; theta0)
556 576 * http://www.aoc.nrao.edu/~egreisen.
557 577 * This method returns the coordinates of the native pole for non-polar
558 578 * projections in spherical coordinates. Units are radians.
559   - * The code has largely be inspired by the IDL routine wcs_getpole.pro
560   - * (update of May 1998).
  579 + * The code has largely been inspired by the IDL routine wcs_getpole.pro
  580 + * (update of May 1998). Is is equivalent to the code found in the
  581 + * wcslib routine cel.c::celset().
561 582 ***************************************************************************/
562 583 GVector GWcs::wcs_getpole(const double& theta0)
563 584 {
... ... @@ -566,37 +587,39 @@ GVector GWcs::wcs_getpole(const double&amp; theta0)
566 587 double delta_p;
567 588  
568 589 // Compute crval in radians
569   - double alpha_0 = m_refval[0] * deg2rad;
570   - double delta_0 = m_refval[1] * deg2rad;
  590 + double alpha_0 = m_refval[0] * deg2rad; // lng0
  591 + double delta_0 = m_refval[1] * deg2rad; // lat0
571 592  
572 593 // Get longitude and latitude of pole
573   - double lonpole = m_npole[0];
574   - double latpole = m_npole[1];
  594 + double lonpole = m_npole[0]; // phip
  595 + double latpole = m_npole[1]; // latp
575 596  
576 597 // If theta0=90 then the coordinate of the native pole are given by crval
  598 + // (fiducial point at the native pole)
577 599 if (theta0 == 90.0) {
578   - alpha_p = alpha_0;
579   - delta_p = delta_0;
  600 + alpha_p = alpha_0; // lngp = lng0
  601 + delta_p = delta_0; // latp = lat0
580 602 }
581 603  
582 604 // ... otherwise perform computations
583 605 else {
584 606 // Compute some angles
585   - double phi_p = lonpole*deg2rad;
586   - double sp = sin(phi_p);
587   - double cp = cos(phi_p);
588   - double sd = sin(delta_0);
589   - double cd = sin(delta_0);
590   - double tand = sin(delta_0);
  607 + double phi_p = lonpole * deg2rad;
  608 + double theta_p = latpole * deg2rad;
  609 + double sp = sin(phi_p);
  610 + double cp = cos(phi_p);
  611 + double sd = sin(delta_0);
  612 + double cd = cos(delta_0);
  613 + double tand = tan(delta_0);
591 614  
592 615 // Case A: theta0 = 0
593 616 if (theta0 == 0.0) {
594   - if (delta_0 == 0.0 && lonpole == 90.0)
595   - delta_p = latpole;
  617 + if (delta_0 == 0.0 && std::abs(lonpole) == 90.0)
  618 + delta_p = theta_p;
596 619 else
597 620 delta_p = acos(sd/cp);
598 621 if (latpole != 90.0) {
599   - if (fabs(latpole + delta_p) < fabs(latpole - delta_p))
  622 + if (std::abs(theta_p + delta_p) < std::abs(theta_p - delta_p))
600 623 delta_p = -delta_p;
601 624 }
602 625 if (lonpole == 180.0 || cd == 0.0)
... ... @@ -625,7 +648,7 @@ GVector GWcs::wcs_getpole(const double&amp; theta0)
625 648 else { // Two valid solutions
626 649 delta_p1 = (term1+term2)*rad2deg;
627 650 delta_p2 = (term1-term2)*rad2deg;
628   - if (fabs(latpole-delta_p1) < fabs(latpole-delta_p2))
  651 + if (std::abs(latpole-delta_p1) < std::abs(latpole-delta_p2))
629 652 delta_p = term1 + term2;
630 653 else
631 654 delta_p = term1 - term2;
... ... @@ -662,7 +685,7 @@ GVector GWcs::wcs_getpole(const double&amp; theta0)
662 685 /***********************************************************************//**
663 686 * @brief Get matrix for rotation between sky and native coordinates
664 687 *
665   - * Requires m_npole and m_native_pole to be set.
  688 + * Requires m_npole[0] and m_native_pole to be set.
666 689 ***************************************************************************/
667 690 GMatrix GWcs::wcs_get_rot(void)
668 691 {
... ... @@ -694,66 +717,6 @@ GMatrix GWcs::wcs_get_rot(void)
694 717  
695 718  
696 719 /***********************************************************************//**
697   - * @brief Set standard WCS parameters
698   - *
699   - * @param[in] coords Coordinate system.
700   - * @param[in] crval1 X value of reference pixel.
701   - * @param[in] crval2 Y value of reference pixel.
702   - * @param[in] crpix1 X index of reference pixel (first pixel is 1).
703   - * @param[in] crpix2 Y index of reference pixel (first pixel is 1).
704   - * @param[in] cdelt1 Increment in x direction at reference pixel (deg).
705   - * @param[in] cdelt2 Increment in y direction at reference pixel (deg).
706   - *
707   - * @exception GException::wcs
708   - * Unable to invert CD matrix.
709   - *
710   - * This method sets the WCS standard parameters as class members. It does
711   - * however not set the rotation matrices. A call to wcs_init is mandatory
712   - * before the projection can be used.
713   - *
714   - * @todo Implement parameter check
715   - ***************************************************************************/
716   -void GWcs::wcs_set(const std::string& coords,
717   - const double& crval1, const double& crval2,
718   - const double& crpix1, const double& crpix2,
719   - const double& cdelt1, const double& cdelt2)
720   -
721   -{
722   - //TODO: Check parameters
723   -
724   - // Set coordinate system
725   - coordsys(coords);
726   -
727   - // Set parameters
728   - m_crval = GVector(crval1, crval2);
729   - m_crpix = GVector(crpix1, crpix2);
730   - m_cdelt = GVector(cdelt1, cdelt2);
731   - m_refval = m_crval;
732   - m_refpix = m_crpix - GVector(1.0,1.0);
733   -
734   - // Set CD matrix corresponding to no rotation
735   - m_cd(0,0) = m_cdelt[0];
736   - m_cd(1,0) = 0.0;
737   - m_cd(0,1) = 0.0;
738   - m_cd(1,1) = m_cdelt[1];
739   -
740   - // Compute inverse CD matrix
741   - double delta = m_cd(0,0)*m_cd(1,1) - m_cd(0,1)*m_cd(1,0);
742   - if (delta != 0.0) {
743   - m_invcd(0,0) = m_cd(1,1) / delta;
744   - m_invcd(0,1) = -m_cd(0,1) / delta;
745   - m_invcd(1,0) = -m_cd(1,0) / delta;
746   - m_invcd(1,1) = m_cd(0,0) / delta;
747   - }
748   - else
749   - throw GException::wcs(G_WCS_SET, "Unable to invert CD matrix.");
750   -
751   - // Return
752   - return;
753   -}
754   -
755   -
756   -/***********************************************************************//**
757 720 * @brief Read WCS definiton from FITS header.
758 721 *
759 722 * @param[in] hdu FITS HDU containing the WCS definition.
... ... @@ -764,11 +727,12 @@ void GWcs::wcs_set(const std::string&amp; coords,
764 727 * Coordinate system is not of valid type.
765 728 *
766 729 * This method reads the WCS definition from the FITS header. It calls the
767   - * wcs_set() method
  730 + * set_members() method.
768 731 *
769 732 * @todo Projection ketwords are not yet extracted from the header. Assume
770 733 * theta0=0.0 for the moment.
771 734 ***************************************************************************/
  735 +/*
772 736 void GWcs::wcs_read(const GFitsHDU* hdu)
773 737 {
774 738 // Continue only if HDU is valid
... ... @@ -808,7 +772,7 @@ void GWcs::wcs_read(const GFitsHDU* hdu)
808 772 throw GException::wcs_bad_coords(G_WCS_READ, coordsys());
809 773  
810 774 // Set standard parameters
811   - wcs_set(coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2);
  775 + set_members(coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2);
812 776  
813 777 // Get projection keywords
814 778 //TODO
... ... @@ -821,7 +785,7 @@ void GWcs::wcs_read(const GFitsHDU* hdu)
821 785 // Return
822 786 return;
823 787 }
824   -
  788 +*/
825 789  
826 790 /***********************************************************************//**
827 791 * @brief Write WCS definiton into FITS HDU
... ... @@ -831,6 +795,7 @@ void GWcs::wcs_read(const GFitsHDU* hdu)
831 795 * @exception GException::wcs_bad_coords
832 796 * Coordinate system is not of valid type.
833 797 ***************************************************************************/
  798 +/*
834 799 void GWcs::wcs_write(GFitsHDU* hdu) const
835 800 {
836 801 // Continue only if HDU is valid
... ... @@ -842,8 +807,6 @@ void GWcs::wcs_write(GFitsHDU* hdu) const
842 807 std::string c_crval2;
843 808 std::string c_cdelt1;
844 809 std::string c_cdelt2;
845   - std::string c_lonpole;
846   - std::string c_latpole;
847 810 if (m_coordsys == 0) {
848 811 c_equinox = "Equinox for celestial coordinate system";
849 812 c_crval1 = "[deg] Right Ascension at reference pixel";
... ... @@ -852,8 +815,6 @@ void GWcs::wcs_write(GFitsHDU* hdu) const
852 815 " reference pixel";
853 816 c_cdelt2 = "[deg/pixel] Declination increment at"
854 817 " reference pixel";
855   - c_lonpole = "[deg] Native longitude of Celestial North Pole";
856   - c_latpole = "[deg] Native latitude of Celestial North Pole";
857 818 }
858 819 else if (m_coordsys == 1) {
859 820 c_equinox = "Equinox for celestial coordinate system";
... ... @@ -863,8 +824,6 @@ void GWcs::wcs_write(GFitsHDU* hdu) const
863 824 " reference pixel";
864 825 c_cdelt2 = "[deg/pixel] Galactic Latitude increment at"
865 826 " reference pixel";
866   - c_lonpole = "[deg] Native longitude of Galactic North Pole";
867   - c_latpole = "[deg] Native latitude of Galactic North Pole";
868 827 }
869 828 else
870 829 throw GException::wcs_bad_coords(G_WCS_WRITE, coordsys());
... ... @@ -882,17 +841,17 @@ void GWcs::wcs_write(GFitsHDU* hdu) const
882 841 hdu->card("CRVAL1", m_crval[0], c_crval1);
883 842 hdu->card("CRVAL2", m_crval[1], c_crval2);
884 843 hdu->card("CROTA2", 0.0, "[deg] Rotation Angle");
885   - //hdu->card("LONPOLE", m_npole[0], c_lonpole);
886   - //hdu->card("LATPOLE", m_npole[1], c_latpole);
887   - hdu->card("PV2_1", 0.0, "Projection parameter 1");
888   - hdu->card("PV2_2", 0.0, "Projection parameter 2");
  844 + hdu->card("LONPOLE", m_lonpole, "[deg] Native longitude of celestial pole");
  845 + hdu->card("LATPOLE", m_latpole, "[deg] Native latitude of celestial pole");
  846 + //hdu->card("PV2_1", 0.0, "Projection parameter 1");
  847 + //hdu->card("PV2_2", 0.0, "Projection parameter 2");
889 848  
890 849 } // endif: HDU was valid
891 850  
892 851 // Return
893 852 return;
894 853 }
895   -
  854 +*/
896 855  
897 856 /***********************************************************************//**
898 857 * @brief Return crval1 string
... ... @@ -959,6 +918,7 @@ std::string GWcs::wcs_crval2(void) const
959 918 /***********************************************************************//**
960 919 * @brief Print WCS information
961 920 ***************************************************************************/
  921 +/*
962 922 std::string GWcs::wcs_dump(void) const
963 923 {
964 924 // Initialise result string
... ... @@ -987,21 +947,116 @@ std::string GWcs::wcs_dump(void) const
987 947 result.append(" "+str(m_invcd(1,1))+"]\n");
988 948  
989 949 // Append coordinate of native pole
990   - result.append(parformat("Coordinate of native pole")+native.print());
  950 + result.append(parformat("Coordinate of native pole")+native.print()+"\n");
  951 +
  952 + // NEW VERSION
  953 +
  954 + // Append World Coordinate parameters
  955 + result.append(parformat("Number of axes")+str(m_naxis)+"\n");
  956 + result.append(parformat("Longitude axis")+str(m_lng)+"\n");
  957 + result.append(parformat("Latitude axis")+str(m_lat)+"\n");
  958 + result.append(parformat("Spectral axis")+str(m_spec)+"\n");
  959 +
  960 + // Append coordinates
  961 + result.append(parformat("Reference coordinate")+"(");
  962 + for (int i = 0; i < m_Crval.size(); ++i) {
  963 + if (i > 0)
  964 + result.append(", ");
  965 + result.append(str(m_Crval[i]));
  966 + if (m_cunit.size() > i)
  967 + result.append(" "+m_cunit[i]);
  968 +
  969 + }
  970 + result.append(")\n");
  971 + result.append(parformat("Reference pixel")+"(");
  972 + for (int i = 0; i < m_Crpix.size(); ++i) {
  973 + if (i > 0)
  974 + result.append(", ");
  975 + result.append(str(m_Crpix[i]));
  976 + }
  977 + result.append(")\n");
  978 + result.append(parformat("Increment at reference")+"(");
  979 + for (int i = 0; i < m_Cdelt.size(); ++i) {
  980 + if (i > 0)
  981 + result.append(", ");
  982 + result.append(str(m_Cdelt[i]));
  983 + if (m_cunit.size() > i)
  984 + result.append(" "+m_cunit[i]);
  985 + }
  986 + result.append(")\n");
  987 +
  988 + // Append celestial transformation parameters
  989 + result.append(parformat("(Phi0, Theta0)")+"(");
  990 + if (undefined(m_phi0))
  991 + result.append("UNDEFINED, ");
  992 + else
  993 + result.append(str(m_phi0)+", ");
  994 + if (undefined(m_theta0))
  995 + result.append("UNDEFINED)\n");
  996 + else
  997 + result.append(str(m_theta0)+") deg\n");
  998 + result.append(parformat("(Lonpole, Latpole)")+"(");
  999 + if (undefined(m_lonpole))
  1000 + result.append("UNDEFINED, ");
  1001 + else
  1002 + result.append(str(m_lonpole)+", ");
  1003 + result.append(str(m_latpole)+") deg\n");
  1004 + result.append(parformat("Ref")+"(");
  1005 + for (int k = 0; k < 4; ++k) {
  1006 + if (k > 0)
  1007 + result.append(", ");
  1008 + if (undefined(m_ref[k]))
  1009 + result.append("UNDEFINED");
  1010 + else
  1011 + result.append(str(m_ref[k]));
  1012 + }
  1013 + result.append(") deg\n");
  1014 + result.append(parformat("Euler")+"(");
  1015 + for (int k = 0; k < 5; ++k) {
  1016 + if (k > 0)
  1017 + result.append(", ");
  1018 + result.append(str(m_euler[k]));
  1019 + }
  1020 + result.append(") deg\n");
  1021 + result.append(parformat("Latpole requirement")+str(m_latpreq)+"\n");
  1022 + if (m_isolat)
  1023 + result.append(parformat("Latitude preserved")+"True\n");
  1024 + else
  1025 + result.append(parformat("Latitude preserved")+"False\n");
  1026 +
  1027 + // Append projection parameters
  1028 + result.append(parformat("Projection code")+m_code+"\n");
  1029 + result.append(parformat("Radius of the gen. sphere")+str(m_r0)+"\n");
  1030 + if (m_bounds)
  1031 + result.append(parformat("Strict bounds checking")+"True\n");
  1032 + else
  1033 + result.append(parformat("Strict bounds checking")+"False\n");
  1034 + if (m_offset)
  1035 + result.append(parformat("Offset")+"True\n");
  1036 + else
  1037 + result.append(parformat("Offset")+"False\n");
  1038 + result.append(parformat("Fiducial offset")+str(m_x0)+", "+str(m_y0)+"\n");
  1039 +
  1040 +
  1041 + // Append spectral transformation parameters
  1042 + result.append(parformat("Rest frequency")+str(m_restfrq)+"\n");
  1043 + result.append(parformat("Rest wavelength")+str(m_restwav)+"\n");
991 1044  
992 1045 // Return
993 1046 return result;
994 1047 }
995   -
  1048 +*/
996 1049  
997 1050 /*==========================================================================
998 1051 = =
999   - = GWcs private methods =
  1052 + = GWcs protected methods =
1000 1053 = =
1001 1054 ==========================================================================*/
1002 1055  
1003 1056 /***********************************************************************//**
1004 1057 * @brief Initialise class members
  1058 + *
  1059 + * Code adapted from wcs.c::wcsini().
1005 1060 ***************************************************************************/
1006 1061 void GWcs::init_members(void)
1007 1062 {
... ... @@ -1017,14 +1072,14 @@ void GWcs::init_members(void)
1017 1072 m_pv2 = GVector(21);
1018 1073  
1019 1074 // Initialise derived members
1020   - m_theta0 = 0.0;
  1075 + //m_theta0 = 0.0;
1021 1076 m_refval = GVector(2);
1022 1077 m_refpix = GVector(2);
1023 1078 m_invcd = GMatrix(2,2);
1024 1079 m_native_pole = GVector(2);
1025 1080 m_rot = GMatrix(3,3);
1026 1081 m_trot = GMatrix(3,3);
1027   -
  1082 +
1028 1083 // Return
1029 1084 return;
1030 1085 }
... ... @@ -1049,7 +1104,7 @@ void GWcs::copy_members(const GWcs&amp; wcs)
1049 1104 m_pv2 = wcs.m_pv2;
1050 1105  
1051 1106 // Copy derived attributes
1052   - m_theta0 = wcs.m_theta0;
  1107 + //m_theta0 = wcs.m_theta0;
1053 1108 m_refval = wcs.m_refval;
1054 1109 m_refpix = wcs.m_refpix;
1055 1110 m_invcd = wcs.m_invcd;
... ...
src/sky/GWcsCAR.cpp
  View file @c1344de
... ... @@ -4,10 +4,18 @@
4 4 * copyright (C) 2010-2011 by Jurgen Knodlseder *
5 5 * ----------------------------------------------------------------------- *
6 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 2 of the License, or *
10   - * (at your option) any later version. *
  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/>. *
11 19 * *
12 20 ***************************************************************************/
13 21 /**
... ... @@ -50,7 +58,7 @@
50 58 /***********************************************************************//**
51 59 * @brief Void constructor
52 60 ***************************************************************************/
53   -GWcsCAR::GWcsCAR(void) : GWcs()
  61 +GWcsCAR::GWcsCAR(void) : GWcslib()
54 62 {
55 63 // Initialise class members
56 64 init_members();
... ... @@ -78,12 +86,18 @@ GWcsCAR::GWcsCAR(const std::string&amp; coords,
78 86 const double& crpix1, const double& crpix2,
79 87 const double& cdelt1, const double& cdelt2,
80 88 const GMatrix& cd, const GVector& pv2) :
81   - GWcs(coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2)
  89 + GWcslib(coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2)
82 90  
83 91 {
84 92 // Initialise class members
85 93 init_members();
86 94  
  95 + // Initialise World Coordinate System. This has to be done here and
  96 + // not in wcs_set(...) as this method will be called by the derived
  97 + // class constructor. We're not allowed to call pure virtual functions
  98 + // in base class constructors !!!
  99 + wcs_set();
  100 +
87 101 // Initialise derived projection parameters
88 102 wcs_init(0.0); // theta0 = 0.0
89 103  
... ... @@ -97,7 +111,7 @@ GWcsCAR::GWcsCAR(const std::string&amp; coords,
97 111 *
98 112 * @param[in] hdu FITS HDU.
99 113 ***************************************************************************/
100   -GWcsCAR::GWcsCAR(const GFitsHDU* hdu) : GWcs()
  114 +GWcsCAR::GWcsCAR(const GFitsHDU* hdu) : GWcslib()
101 115 {
102 116 // Initialise class members
103 117 init_members();
... ... @@ -115,7 +129,7 @@ GWcsCAR::GWcsCAR(const GFitsHDU* hdu) : GWcs()
115 129 *
116 130 * @param[in] wcs World Coordinate System.
117 131 ***************************************************************************/
118   -GWcsCAR::GWcsCAR(const GWcsCAR& wcs) : GWcs(wcs)
  132 +GWcsCAR::GWcsCAR(const GWcsCAR& wcs) : GWcslib(wcs)
119 133 {
120 134 // Initialise class members for clean destruction
121 135 init_members();
... ... @@ -158,7 +172,7 @@ GWcsCAR&amp; GWcsCAR::operator= (const GWcsCAR&amp; wcs)
158 172 if (this != &wcs) {
159 173  
160 174 // Copy base class members
161   - this->GWcs::operator=(wcs);
  175 + this->GWcslib::operator=(wcs);
162 176  
163 177 // Free members
164 178 free_members();
... ... @@ -191,10 +205,12 @@ void GWcsCAR::clear(void)
191 205 {
192 206 // Free class members (base and derived classes, derived class first)
193 207 free_members();
  208 + this->GWcslib::free_members();
194 209 this->GWcs::free_members();
195 210  
196 211 // Initialise members
197 212 this->GWcs::init_members();
  213 + this->GWcslib::init_members();
198 214 init_members();
199 215  
200 216 // Return
... ... @@ -216,6 +232,7 @@ GWcsCAR* GWcsCAR::clone(void) const
216 232 *
217 233 * @param[in] hdu FITS HDU.
218 234 ***************************************************************************/
  235 +/*
219 236 void GWcsCAR::read(const GFitsHDU* hdu)
220 237 {
221 238 // Clear object
... ... @@ -232,13 +249,14 @@ void GWcsCAR::read(const GFitsHDU* hdu)
232 249 // Return
233 250 return;
234 251 }
235   -
  252 +*/
236 253  
237 254 /***********************************************************************//**
238 255 * @brief Write Healpix definiton into FITS HDU
239 256 *
240 257 * @param[in] hdu FITS HDU.
241 258 ***************************************************************************/
  259 +/*
242 260 void GWcsCAR::write(GFitsHDU* hdu) const
243 261 {
244 262 // Continue only if HDU pointer is valid
... ... @@ -252,7 +270,7 @@ void GWcsCAR::write(GFitsHDU* hdu) const
252 270 // Return
253 271 return;
254 272 }
255   -
  273 +*/
256 274  
257 275 /***********************************************************************//**
258 276 * @brief Returns solid angle of pixel in units of steradians
... ... @@ -273,6 +291,7 @@ void GWcsCAR::write(GFitsHDU* hdu) const
273 291 * @todo Implement accurate solid angle computation (so far only brute force
274 292 * estimation)
275 293 ***************************************************************************/
  294 +/*
276 295 double GWcsCAR::omega(const GSkyPixel& pix) const
277 296 {
278 297 // Bypass const correctness
... ... @@ -297,7 +316,7 @@ double GWcsCAR::omega(const GSkyPixel&amp; pix) const
297 316 // Return solid angle
298 317 return omega;
299 318 }
300   -
  319 +*/
301 320  
302 321 /***********************************************************************//**
303 322 * @brief Print WCS information
... ... @@ -309,7 +328,7 @@ std::string GWcsCAR::print(void) const
309 328  
310 329 // Append header
311 330 result.append("=== GWcsCAR ===\n");
312   - result.append(wcs_dump());
  331 + result.append(wcs_print());
313 332  
314 333 // Return result
315 334 return result;
... ... @@ -324,12 +343,30 @@ std::string GWcsCAR::print(void) const
324 343  
325 344 /***********************************************************************//**
326 345 * @brief Initialise class members
  346 + *
  347 + * Given and/or returned:
  348 + * m_r0 Reset to 180/pi if 0.
  349 + * m_phi0 Reset to 0.0 if undefined.
  350 + * m_theta0 Reset to 0.0 if undefined.
  351 + *
  352 + * Returned:
  353 + * m_code "CAR"
  354 + * m_name "plate caree"
  355 + * m_x0 Fiducial offset in x.
  356 + * m_y0 Fiducial offset in y.
  357 + * m_w0 r0*(pi/180)
  358 + * m_w1 (180/pi)/r0
327 359 ***************************************************************************/
328 360 void GWcsCAR::init_members(void)
329 361 {
330 362 // Initialise members
331 363 m_type = "CAR";
332 364  
  365 + // NEW VERSION
  366 +
  367 + // Initialise projection parameters
  368 + prj_set();
  369 +
333 370 // Return
334 371 return;
335 372 }
... ... @@ -393,6 +430,184 @@ void GWcsCAR::nat2std(GVector *coord) const
393 430 }
394 431  
395 432  
  433 +/***********************************************************************//**
  434 + * @brief Setup of projection
  435 + *
  436 + * This method sets up the projection information. The method has been
  437 + * adapted from the wcslib function prj.c::carset.
  438 + ***************************************************************************/
  439 +void GWcsCAR::prj_set(void) const
  440 +{
  441 + // Initialise projection parameters
  442 + m_code = "CAR";
  443 + m_name = "plate caree";
  444 +
  445 + // Precompute
  446 + if (m_r0 == 0.0) {
  447 + m_r0 = rad2deg;
  448 + m_w0 = 1.0;
  449 + m_w1 = 1.0;
  450 + } else {
  451 + m_w0 = m_r0 * deg2rad;
  452 + m_w1 = 1.0/m_w0;
  453 + }
  454 +
  455 + // Compute fiducial offset
  456 + prj_off(0.0, 0.0);
  457 +
  458 + // Signal that projection has been set
  459 + m_prjset = true;
  460 +
  461 + // Return
  462 + return;
  463 +}
  464 +
  465 +
  466 +/***********************************************************************//**
  467 + * @brief Cartesian-to-spherical deprojection
  468 + *
  469 + * @param[in] nx X vector length.
  470 + * @param[in] ny Y vector length (0=no replication).
  471 + * @param[in] sxy Input vector step.
  472 + * @param[in] spt Output vector step.
  473 + * @param[in] x Vector of projected x coordinates.
  474 + * @param[in] y Vector of projected y coordinates.
  475 + * @param[out] phi Longitude of the projected point in native spherical
  476 + * coordinates [deg].
  477 + * @param[out] theta Latitude of the projected point in native spherical
  478 + * coordinates [deg].
  479 + * @param[out] stat Status return value for each vector element (always 0)
  480 + *
  481 + * Deproject Cartesian (x,y) coordinates in the plane of projection to native
  482 + * spherical coordinates (phi,theta).
  483 + *
  484 + * This method has been adapted from the wcslib function prj.c::carx2s().
  485 + * The interface follows very closely that of wcslib. In contrast to the
  486 + * wcslib routine, however, the method assumes that the projection has been
  487 + * setup previsouly (as this will be done by the constructor).
  488 + ***************************************************************************/
  489 +void GWcsCAR::prj_x2s(int nx, int ny, int sxy, int spt,
  490 + const double* x, const double* y,
  491 + double* phi, double* theta, int* stat) const
  492 +{
  493 + // Initialize projection if required
  494 + if (!m_prjset)
  495 + prj_set();
  496 +
  497 + // Set value replication length mx,my
  498 + int mx;
  499 + int my;
  500 + if (ny > 0) {
  501 + mx = nx;
  502 + my = ny;
  503 + }
  504 + else {
  505 + mx = 1;
  506 + my = 1;
  507 + ny = nx;
  508 + }
  509 +
  510 + // Do x dependence
  511 + const double* xp = x;
  512 + int rowoff = 0;
  513 + int rowlen = nx * spt;
  514 + for (int ix = 0; ix < nx; ++ix, rowoff += spt, xp += sxy) {
  515 + double s = m_w1 * (*xp + m_x0);
  516 + double* phip = phi + rowoff;
  517 + for (int iy = 0; iy < my; ++iy, phip += rowlen)
  518 + *phip = s;
  519 + }
  520 +
  521 + // Do y dependence
  522 + const double* yp = y;
  523 + double* thetap = theta;
  524 + int* statp = stat;
  525 + for (int iy = 0; iy < ny; ++iy, yp += sxy) {
  526 + double t = m_w1 * (*yp + m_y0);
  527 + for (int ix = 0; ix < mx; ++ix, thetap += spt) {
  528 + *thetap = t;
  529 + *(statp++) = 0;
  530 + }
  531 + }
  532 +
  533 + // Return
  534 + return;
  535 +}
  536 +
  537 +
  538 +/***********************************************************************//**
  539 + * @brief Generic spherical-to-Cartesian projection
  540 + *
  541 + * @param[in] nphi Longitude vector length.
  542 + * @param[in] ntheta Latitude vector length (0=no replication).
  543 + * @param[in] spt Input vector step.
  544 + * @param[in] sxy Output vector step.
  545 + * @param[in] phi Longitude vector of the projected point in native spherical
  546 + * coordinates [deg].
  547 + * @param[in] theta Latitude vector of the projected point in native spherical
  548 + * coordinates [deg].
  549 + * @param[out] x Vector of projected x coordinates.
  550 + * @param[out] y Vector of projected y coordinates.
  551 + * @param[out] stat Status return value for each vector element (always 0)
  552 + *
  553 + * Project native spherical coordinates (phi,theta) to Cartesian (x,y)
  554 + * coordinates in the plane of projection.
  555 + *
  556 + * This method has been adapted from the wcslib function prj.c::cars2x().
  557 + * The interface follows very closely that of wcslib. In contrast to the
  558 + * wcslib routine, however, the method assumes that the projection has been
  559 + * setup previsouly (as this will be done by the constructor).
  560 + ***************************************************************************/
  561 +void GWcsCAR::prj_s2x(int nphi, int ntheta, int spt, int sxy,
  562 + const double* phi, const double* theta,
  563 + double* x, double* y, int* stat) const
  564 +{
  565 + // Initialize projection if required
  566 + if (!m_prjset)
  567 + prj_set();
  568 +
  569 + // Set value replication length mphi,mtheta
  570 + int mphi;
  571 + int mtheta;
  572 + if (ntheta > 0) {
  573 + mphi = nphi;
  574 + mtheta = ntheta;
  575 + }
  576 + else {
  577 + mphi = 1;
  578 + mtheta = 1;
  579 + ntheta = nphi;
  580 + }
  581 +
  582 + // Do phi dependence
  583 + const double* phip = phi;
  584 + int rowoff = 0;
  585 + int rowlen = nphi * sxy;
  586 + for (int iphi = 0; iphi < nphi; ++iphi, rowoff += sxy, phip += spt) {
  587 + double xi = m_w0 * (*phip) - m_x0;
  588 + double* xp = x + rowoff;
  589 + for (int itheta = 0; itheta < mtheta; ++itheta, xp += rowlen)
  590 + *xp = xi;
  591 + }
  592 +
  593 +
  594 + // Do theta dependence
  595 + const double* thetap = theta;
  596 + double* yp = y;
  597 + int* statp = stat;
  598 + for (int itheta = 0; itheta < ntheta; ++itheta, thetap += spt) {
  599 + double eta = m_w0 * (*thetap) - m_y0;
  600 + for (int iphi = 0; iphi < mphi; ++iphi, yp += sxy) {
  601 + *yp = eta;
  602 + *(statp++) = 0;
  603 + }
  604 + }
  605 +
  606 + // Return
  607 + return;
  608 +}
  609 +
  610 +
396 611 /*==========================================================================
397 612 = =
398 613 = Friends =
... ...
src/sky/GWcsHPX.cpp
  View file @c1344de
... ... @@ -900,6 +900,46 @@ unsigned int GWcsHPX::isqrt(unsigned int arg) const
900 900 }
901 901  
902 902  
  903 +/***********************************************************************//**
  904 + * @brief Setup of projection
  905 + *
  906 + * Dummy (only implemented for non HPX classes).
  907 + ***************************************************************************/
  908 +void GWcsHPX::prj_set(void)
  909 +{
  910 + // Return
  911 + return;
  912 +}
  913 +
  914 +
  915 +/***********************************************************************//**
  916 + * @brief Cartesian-to-spherical deprojection
  917 + *
  918 + * Dummy (only implemented for non HPX classes).
  919 + ***************************************************************************/
  920 +int GWcsHPX::prj_x2s(int nx, int ny, int sxy, int spt,
  921 + const double* x, const double* y,
  922 + double* phi, double* theta, int* stat) const
  923 +{
  924 + // Return
  925 + return 0;
  926 +}
  927 +
  928 +
  929 +/***********************************************************************//**
  930 + * @brief Generic spherical-to-Cartesian projection
  931 + *
  932 + * Dummy (only implemented for non HPX classes).
  933 + ***************************************************************************/
  934 +int GWcsHPX::prj_s2x(int nphi, int ntheta, int spt, int sxy,
  935 + const double* phi, const double* theta,
  936 + double* x, double* y, int* stat) const
  937 +{
  938 + // Return
  939 + return 0;
  940 +}
  941 +
  942 +
903 943 /*==========================================================================
904 944 = =
905 945 = Friends =
... ...
test/test_GSky.cpp
  View file @c1344de
1 1 /***************************************************************************
2 2 * test_GSky.cpp - test GSky classes *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010 by Jurgen Knodlseder *
  4 + * copyright (C) 2010-2011 by Jurgen Knodlseder *
5 5 * ----------------------------------------------------------------------- *
6 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 2 of the License, or *
10   - * (at your option) any later version. *
  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/>. *
11 19 * *
12 20 ***************************************************************************/
13 21 /**
... ... @@ -31,6 +39,132 @@
31 39  
32 40  
33 41 /***************************************************************************
  42 + * Test consistency of forward and background transformations *
  43 + ***************************************************************************/
  44 +double wcs_forth_back_pixel(GWcs& wcs, double& crpix1, double& crpix2)
  45 +{
  46 + // Initialise maximal distance
  47 + double dist_max = 0.0;
  48 +
  49 + // Initialise number of pixels
  50 + int nx = 100;
  51 + int ny = 100;
  52 +
  53 + // Loop over x pixels
  54 + for (int ix = -nx; ix <= nx; ++ix) {
  55 +
  56 + // Set x value
  57 + double x = double(ix) + crpix1;
  58 +
  59 + // Loop over y pixels
  60 + for (int iy = -ny; iy <= ny; ++iy) {
  61 +
  62 + // Set y value
  63 + double y = double(iy) + crpix2;
  64 +
  65 + // Set sky pixel
  66 + GSkyPixel pix_in(x,y);
  67 +
  68 + // Forth: Transform to world
  69 + GSkyDir dir = wcs.xy2dir(pix_in);
  70 +
  71 + // Back: Transform to pixel
  72 + GSkyPixel pix_out = wcs.dir2xy(dir);
  73 +
  74 + // Compute distance
  75 + double dx = pix_in.x()-pix_out.x();
  76 + double dy = pix_in.y()-pix_out.y();
  77 + double dist = sqrt(dx*dx+dy*dy);
  78 + //std::cout << dist << std::endl;
  79 +
  80 + // Store maximum distance
  81 + if (dist > dist_max)
  82 + dist_max = dist;
  83 +
  84 + } // endfor: y pixels
  85 + } // endfor: x pixels
  86 +
  87 + // Return
  88 + return dist_max;
  89 +}
  90 +
  91 +
  92 +/***************************************************************************
  93 + * Test: GWcsCAR *
  94 + ***************************************************************************/
  95 +void test_GWcsCAR(void)
  96 +{
  97 + // Dump header
  98 + std::cout << "Test GWcsCAR: ";
  99 +
  100 + // Allocate parameters for testing
  101 + std::string coords = "CEL";
  102 + double crval1 = 83.63;
  103 + double crval2 = 22.01;
  104 + double crpix1 = 100.5;
  105 + double crpix2 = 100.5;
  106 + double cdelt1 = 0.02;
  107 + double cdelt2 = 0.02;
  108 + GMatrix cd(2,2);
  109 + GVector pv2(21);
  110 +
  111 + // Test void constructor
  112 + try {
  113 + GWcsCAR wcs;
  114 + }
  115 + catch (std::exception &e) {
  116 + std::cout << std::endl
  117 + << "TEST ERROR: Unable to construct empty GWcsCAR."
  118 + << std::endl;
  119 + std::cout << e.what() << std::endl;
  120 + throw;
  121 + }
  122 + std::cout << ".";
  123 +
  124 + // Test constructor
  125 + try {
  126 + GWcsCAR wcs(coords, crval1, crval2, crpix1, crpix2,
  127 + cdelt1, cdelt2, cd, pv2);
  128 + std::cout << wcs << std::endl;
  129 + }
  130 + catch (std::exception &e) {
  131 + std::cout << std::endl
  132 + << "TEST ERROR: Unable to construct GWcsCAR."
  133 + << std::endl;
  134 + std::cout << e.what() << std::endl;
  135 + throw;
  136 + }
  137 + std::cout << ".";
  138 +
  139 + // Build now skymap for conversion testing
  140 + GWcsCAR car(coords, crval1, crval2, crpix1, crpix2,
  141 + cdelt1, cdelt2, cd, pv2);
  142 +
  143 + // Test (x,y) to (lon,lat)
  144 + try {
  145 + double tol = 0.0;
  146 + if ((tol = wcs_forth_back_pixel(car, crpix1, crpix2)) > 1.0e-10) {
  147 + std::cout << std::endl
  148 + << "TEST ERROR: GWcsCAR transformation tolerance "
  149 + << tol << std::endl;
  150 + throw;
  151 + }
  152 + }
  153 + catch (std::exception &e) {
  154 + std::cout << std::endl
  155 + << "TEST ERROR: Error while testing sky conversion."
  156 + << std::endl;
  157 + std::cout << e.what() << std::endl;
  158 + throw;
  159 + }
  160 + std::cout << ".";
  161 +
  162 + // Exit test
  163 + return;
  164 +}
  165 +
  166 +
  167 +/***************************************************************************
34 168 * Test: GSkymap_healpix_construct *
35 169 ***************************************************************************/
36 170 void test_GSkymap_healpix_construct(void)
... ... @@ -468,10 +602,11 @@ int main(void)
468 602 std::cout << "************************" << std::endl;
469 603  
470 604 // Execute Healpix tests
471   - test_GSkymap_healpix_construct();
472   - test_GSkymap_healpix_io();
473   - test_GSkymap_wcs_construct();
474   - test_GSkymap_wcs_io();
  605 + test_GWcsCAR();
  606 + //test_GSkymap_healpix_construct();
  607 + //test_GSkymap_healpix_io();
  608 + //test_GSkymap_wcs_construct();
  609 + //test_GSkymap_wcs_io();
475 610  
476 611 // Return
477 612 return 0;
... ...