Dmitry Malyshev / gammalib

Commit 5db80affa598f1703da833e8fcd2c2431113a780

Authored by knodlseder
1 parent 0775ee7d

update for new wcs design (wcslib adaption for all non-Healpix classes)

Inline Side-by-side
Showing 3 changed files with 581 additions and 396 deletions
include/GWcslib.hpp
  View file @5db80af
... ... @@ -41,14 +41,10 @@
41 41 /***********************************************************************//**
42 42 * @class GWcslib
43 43 *
44   - * @brief GWcslib virtual base class interface defintion
  44 + * @brief Virtual base class for wcslib based World Coordinate System
45 45 ***************************************************************************/
46 46 class GWcslib : public GWcs {
47 47  
48   - // Operator friends
49   - friend bool operator== (const GWcslib &a, const GWcslib &b);
50   - friend bool operator!= (const GWcslib &a, const GWcslib &b);
51   -
52 48 public:
53 49 // Constructors and destructors
54 50 GWcslib(void);
... ... @@ -56,6 +52,7 @@ public:
56 52 const double& crval1, const double& crval2,
57 53 const double& crpix1, const double& crpix2,
58 54 const double& cdelt1, const double& cdelt2);
  55 + explicit GWcslib(const GFitsHDU* hdu);
59 56 GWcslib(const GWcslib& wcs);
60 57 virtual ~GWcslib(void);
61 58  
... ... @@ -65,14 +62,25 @@ public:
65 62 // Pure virtual methods (not implemented)
66 63 virtual void clear(void) = 0;
67 64 virtual GWcslib* clone(void) const = 0;
  65 + virtual std::string code(void) const = 0;
  66 + virtual std::string name(void) const = 0;
68 67 virtual std::string print(void) const = 0;
69 68  
70   - // Implemented methods
71   - void read(const GFitsHDU* hdu);
72   - void write(GFitsHDU* hdu) const;
73   - double omega(const GSkyPixel& pix) const;
74   - GSkyDir xy2dir(const GSkyPixel& pix) const;
75   - GSkyPixel dir2xy(GSkyDir dir) const;
  69 + // Implemented virtual methods
  70 + virtual void read(const GFitsHDU* hdu);
  71 + virtual void write(GFitsHDU* hdu) const;
  72 + virtual double omega(const int& pix) const;
  73 + virtual double omega(const GSkyPixel& pix) const;
  74 + virtual GSkyDir pix2dir(const int& pix) const;
  75 + virtual int dir2pix(const GSkyDir& dir) const;
  76 + virtual GSkyDir xy2dir(const GSkyPixel& pix) const;
  77 + virtual GSkyPixel dir2xy(const GSkyDir& dir) const;
  78 +
  79 + // Other methods
  80 + void set(const std::string& coords,
  81 + const double& crval1, const double& crval2,
  82 + const double& crpix1, const double& crpix2,
  83 + const double& cdelt1, const double& cdelt2);
76 84  
77 85 private:
78 86 // Constants
... ... @@ -81,24 +89,27 @@ private:
81 89  
82 90 protected:
83 91 // Protected methods
84   - void init_members(void);
85   - void copy_members(const GWcslib& wcs);
86   - void free_members(void);
87   - bool undefined(const double& value) const { return (value==UNDEFINED); }
88   - void set_members(const std::string& coords,
89   - const double& crval1, const double& crval2,
90   - const double& crpix1, const double& crpix2,
91   - const double& cdelt1, const double& cdelt2);
  92 + void init_members(void);
  93 + void copy_members(const GWcslib& wcs);
  94 + void free_members(void);
  95 + void set_members(const std::string& coords,
  96 + const double& crval1, const double& crval2,
  97 + const double& crpix1, const double& crpix2,
  98 + const double& cdelt1, const double& cdelt2);
  99 + virtual bool compare(const GWcs& wcs) const;
  100 + bool undefined(const double& value) const { return (value==UNDEFINED); }
92 101  
93 102 // Methods adapted from wcslib::wcs.c
94 103 void wcs_ini(int naxis);
95 104 void wcs_set(void) const;
  105 + void wcs_set_ctype(void) const;
96 106 void wcs_p2s(int ncoord, int nelem, const double* pixcrd, double* imgcrd,
97 107 double* phi, double* theta, double* world, int* stat) const;
98 108 void wcs_s2p(int ncoord, int nelem, const double* world,
99 109 double* phi, double* theta, double* imgcrd,
100 110 double* pixcrd, int* stat) const;
101 111 std::string wcs_print(void) const;
  112 + std::string wcs_print_value(const double& value) const;
102 113  
103 114 // Methods adapted from wcslib::cel.c
104 115 void cel_ini(void) const;
... ... @@ -152,51 +163,53 @@ protected:
152 163 void sincosd(const double& angle, double *s, double *c) const;
153 164  
154 165 // World Coordinate System parameters
155   - mutable bool m_wcsset; //!< WCS information is set
156   - int m_naxis; //!< Number of axes
157   - std::vector<double> m_Crval; //!< CRVALia keyvalues for each coord axis
158   - std::vector<std::string> m_cunit; //!< CUNITia keyvalues for each coord axis
159   - std::vector<std::string> m_ctype; //!< CTYPEia keyvalues for each coord axis
160   - mutable double m_lonpole; //!< LONPOLEa keyvalue
161   - mutable double m_latpole; //!< LATPOLEa keyvalue
162   - double m_restfrq; //!< RESTFRQa keyvalue
163   - double m_restwav; //!< RESTWAVa keyvalue
164   - std::vector<double> m_Cd; //!< CDi_ja linear transformation matrix
165   - std::vector<double> m_crota; //!< CROTAia keyvalues for each coord axis
166   - int m_lng; //!< Longitude axis
167   - int m_lat; //!< Latitude axis
168   - int m_spec; //!< Spectral axis
  166 + mutable bool m_wcsset; //!< WCS information is set
  167 + int m_naxis; //!< Number of axes
  168 + std::vector<double> m_crval; //!< CRVALia keyvalues for each coord axis
  169 + std::vector<std::string> m_cunit; //!< CUNITia keyvalues for each coord axis
  170 + mutable std::vector<std::string> m_ctype; //!< CTYPEia keyvalues for each coord axis
  171 + mutable std::vector<std::string> m_ctype_c; //!< CTYPEia comments for each coord axis
  172 + mutable double m_lonpole; //!< LONPOLEa keyvalue
  173 + mutable double m_latpole; //!< LATPOLEa keyvalue
  174 + double m_restfrq; //!< RESTFRQa keyvalue
  175 + double m_restwav; //!< RESTWAVa keyvalue
  176 + std::string m_radesys; //!< RADESYS keyvalue
  177 + double m_equinox; //!< EQUINOX keyvalue
  178 + std::vector<double> m_cd; //!< CDi_ja linear transformation matrix
  179 + std::vector<double> m_crota; //!< CROTAia keyvalues for each coord axis
  180 + int m_lng; //!< Longitude axis
  181 + int m_lat; //!< Latitude axis
  182 + int m_spec; //!< Spectral axis
169 183  
170 184 // Linear transformation parameters
171   - mutable bool m_linset; //!< Linear transformation is set
172   - mutable bool m_unity; //!< Signals unity PC matrix
173   - std::vector<double> m_Crpix; //!< CRPIXja keyvalues for each pixel axis
174   - std::vector<double> m_pc; //!< PCi_ja linear transformation matrix
175   - std::vector<double> m_Cdelt; //!< CDELTia keyvalues for each coord axis
176   - mutable std::vector<double> m_piximg; //!< Pixel to image transformation matrix
177   - mutable std::vector<double> m_imgpix; //!< Image to pixel transformation matrix
  185 + mutable bool m_linset; //!< Linear transformation is set
  186 + mutable bool m_unity; //!< Signals unity PC matrix
  187 + std::vector<double> m_crpix; //!< CRPIXja keyvalues for each pixel axis
  188 + std::vector<double> m_pc; //!< PCi_ja linear transformation matrix
  189 + std::vector<double> m_cdelt; //!< CDELTia keyvalues for each coord axis
  190 + mutable std::vector<double> m_piximg; //!< Pixel to image transformation matrix
  191 + mutable std::vector<double> m_imgpix; //!< Image to pixel transformation matrix
178 192  
179 193 // Celestial transformation parameters
180   - mutable bool m_celset; //!< Celestial transformation is set
181   - mutable bool m_offset; //!< Force (x,y) = (0,0) at (phi_0,theta_0)
182   - mutable double m_phi0; //!< Native azimuth angle of fiducial point
183   - mutable double m_theta0; //!< Native zenith angle of fiducial point
184   - mutable double m_ref[4]; //!< Celestial coordinates of fiducial
185   - // point and native coordinates of celestial
186   - // pole
187   - mutable double m_euler[5];//!< Euler angles and functions thereof
188   - mutable int m_latpreq; //!< LATPOLEa requirement
189   - mutable bool m_isolat; //!< True if |latitude| is preserved
  194 + mutable bool m_celset; //!< Celestial transformation is set
  195 + mutable bool m_offset; //!< Force (x,y) = (0,0) at (phi_0,theta_0)
  196 + mutable double m_phi0; //!< Native azimuth angle of fiducial point
  197 + mutable double m_theta0; //!< Native zenith angle of fiducial point
  198 + mutable double m_ref[4]; //!< Celestial coordinates of fiducial
  199 + // point and native coordinates of celestial
  200 + // pole
  201 + mutable double m_euler[5];//!< Euler angles and functions thereof
  202 + mutable int m_latpreq; //!< LATPOLEa requirement
  203 + mutable bool m_isolat; //!< True if |latitude| is preserved
190 204  
191 205 // Projection parameters
192   - mutable bool m_prjset; //!< Projection is set
193   - mutable std::string m_code; //!< Three-letter projection code
194   - mutable std::string m_name; //!< Projection name
195   - mutable double m_r0; //!< Radius of the generating sphere
196   - mutable double m_pv[PVN]; //!< Projection parameters
197   - mutable bool m_bounds; //!< Enable strict bounds checking
198   - mutable double m_x0; //!< Fiducial x offset
199   - mutable double m_y0; //!< Fiducial y offset
  206 + mutable bool m_prjset; //!< Projection is set
  207 + mutable double m_r0; //!< Radius of the generating sphere
  208 + mutable double m_pv[PVN]; //!< Projection parameters
  209 + mutable bool m_bounds; //!< Enable strict bounds checking
  210 + mutable double m_x0; //!< Fiducial x offset
  211 + mutable double m_y0; //!< Fiducial y offset
  212 + mutable std::vector<double> m_w; //!< Intermediate values
200 213  
201 214 // Spectral transformation parameters
202 215 //struct spcprm spc
... ...
src/sky/GWcslib.cpp
  View file @5db80af
... ... @@ -35,9 +35,10 @@
35 35  
36 36 /* __ Method name definitions ____________________________________________ */
37 37 #define G_READ "GWcslib::read(GFitsHDU*)"
38   -#define G_WRITE "GWcslib::write(GFitsHDU*)"
39   -#define G_WCS_SET "GWcslib::GWcslib(std::string,double,double,double,double,"\
40   - "double,double)"
  38 +#define G_OMEGA "GWcslib::omega(int&)"
  39 +#define G_PIX2DIR "GWcslib::pix2dir(int&)"
  40 +#define G_DIR2PIX "GWcslib::dir2pix(GSkyDir&)"
  41 +#define G_WCS_SET_CTYPE "GWcslib::wcs_set_ctype()"
41 42 #define G_WCS_P2S "GWcslib::wcs_s2p(int,int,double*,double*,double*,double*,"\
42 43 "double*,int*)"
43 44 #define G_WCS_S2P "GWcslib::wcs_s2p(int,int,double*,double*,double*,double*,"\
... ... @@ -70,7 +71,7 @@
70 71 /***********************************************************************//**
71 72 * @brief Void constructor
72 73 ***************************************************************************/
73   -GWcslib::GWcslib(void)
  74 +GWcslib::GWcslib(void) : GWcs()
74 75 {
75 76 // Initialise class members
76 77 init_members();
... ... @@ -96,14 +97,11 @@ GWcslib::GWcslib(void)
96 97 GWcslib::GWcslib(const std::string& coords,
97 98 const double& crval1, const double& crval2,
98 99 const double& crpix1, const double& crpix2,
99   - const double& cdelt1, const double& cdelt2)
  100 + const double& cdelt1, const double& cdelt2) : GWcs()
100 101 {
101 102 // Initialise class members
102 103 init_members();
103   -
104   - // Initialise reverse flag
105   - m_reverse = 0;
106   -
  104 +
107 105 // Set standard parameters
108 106 set_members(coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2);
109 107  
... ... @@ -113,6 +111,24 @@ GWcslib::GWcslib(const std::string&amp; coords,
113 111  
114 112  
115 113 /***********************************************************************//**
  114 + * @brief Construct from FITS HDU table
  115 + *
  116 + * @param[in] hdu FITS HDU.
  117 + ***************************************************************************/
  118 +GWcslib::GWcslib(const GFitsHDU* hdu) : GWcs()
  119 +{
  120 + // Initialise class members
  121 + init_members();
  122 +
  123 + // Read WCS definition from FITS HDU
  124 + read(hdu);
  125 +
  126 + // Return
  127 + return;
  128 +}
  129 +
  130 +
  131 +/***********************************************************************//**
116 132 * @brief Copy constructor
117 133 *
118 134 * @param wcs World Coordinate System.
... ... @@ -189,16 +205,10 @@ GWcslib&amp; GWcslib::operator= (const GWcslib&amp; wcs)
189 205 *
190 206 * @param[in] hdu FITS HDU containing the WCS definition.
191 207 *
192   - * @exception GException::fits_key_not_found
193   - * Unable to find required FITS header keyword.
194 208 * @exception GException::wcs_bad_coords
195 209 * Coordinate system is not of valid type.
196 210 *
197   - * This method reads the WCS definition from the FITS header. It calls the
198   - * set_members() method.
199   - *
200   - * @todo Projection ketwords are not yet extracted from the header. Assume
201   - * theta0=0.0 for the moment.
  211 + * This method reads the WCS definition from the FITS header.
202 212 ***************************************************************************/
203 213 void GWcslib::read(const GFitsHDU* hdu)
204 214 {
... ... @@ -223,20 +233,54 @@ void GWcslib::read(const GFitsHDU* hdu)
223 233 std::string xcoord = ctype1.substr(0,4);
224 234 std::string ycoord = ctype2.substr(0,4);
225 235 if (xcoord == "RA--" && ycoord == "DEC-") {
226   - m_reverse = 0;
227   - coords = "EQU";
  236 + m_lng = 0;
  237 + m_lat = 1;
  238 + coords = "EQU";
228 239 }
229 240 else if (xcoord == "DEC-" && ycoord == "RA--") {
230   - m_reverse = 1;
231   - coords = "EQU";
  241 + m_lng = 1;
  242 + m_lat = 0;
  243 + coords = "EQU";
232 244 }
233 245 else if (xcoord == "GLON" && ycoord == "GLAT") {
234   - m_reverse = 0;
235   - coords = "GAL";
  246 + m_lng = 0;
  247 + m_lat = 1;
  248 + coords = "GAL";
236 249 }
237 250 else if (xcoord == "GLAT" && ycoord == "GLON") {
238   - m_reverse = 1;
239   - coords = "GAL";
  251 + m_lng = 1;
  252 + m_lat = 0;
  253 + coords = "GAL";
  254 + }
  255 + else if (xcoord == "ELON" && ycoord == "ELAT") {
  256 + m_lng = 0;
  257 + m_lat = 1;
  258 + coords = "ECL";
  259 + }
  260 + else if (xcoord == "ELAT" && ycoord == "ELON") {
  261 + m_lng = 1;
  262 + m_lat = 0;
  263 + coords = "ECL";
  264 + }
  265 + else if (xcoord == "HLON" && ycoord == "HLAT") {
  266 + m_lng = 0;
  267 + m_lat = 1;
  268 + coords = "HEL";
  269 + }
  270 + else if (xcoord == "HLAT" && ycoord == "HLON") {
  271 + m_lng = 1;
  272 + m_lat = 0;
  273 + coords = "HEL";
  274 + }
  275 + else if (xcoord == "SLON" && ycoord == "SLAT") {
  276 + m_lng = 0;
  277 + m_lat = 1;
  278 + coords = "SGL";
  279 + }
  280 + else if (xcoord == "SLAT" && ycoord == "SLON") {
  281 + m_lng = 1;
  282 + m_lat = 0;
  283 + coords = "SGL";
240 284 }
241 285 else
242 286 throw GException::wcs_bad_coords(G_READ, coordsys());
... ... @@ -244,12 +288,6 @@ void GWcslib::read(const GFitsHDU* hdu)
244 288 // Set standard parameters
245 289 set_members(coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2);
246 290  
247   - //TODO: Get projection keywords
248   -
249   - // Initialise derived projection parameters
250   - wcs_set();
251   - wcs_init(0.0); // theta0 = 0.0 <= THIS DISAPPEARS LATER
252   -
253 291 } // endif: HDU was valid
254 292  
255 293 // Return
... ... @@ -258,62 +296,119 @@ void GWcslib::read(const GFitsHDU* hdu)
258 296  
259 297  
260 298 /***********************************************************************//**
261   - * @brief Write WCS definiton into FITS HDU
  299 + * @brief Write WCS definiton into FITS HDU header
262 300 *
263   - * @param[in] hdu FITS HDU into which WCS definition will be written.
  301 + * @param[in] hdu FITS HDU.
264 302 *
265   - * @exception GException::wcs_bad_coords
266   - * Coordinate system is not of valid type.
  303 + * This method writes the World Coordinate System definition into the FITS
  304 + * HDU header.
  305 + *
  306 + * The method does nothing if the HDU is not valid.
  307 + *
  308 + * This method has been adapted from wcshdr.c::wcshdo().
267 309 ***************************************************************************/
268 310 void GWcslib::write(GFitsHDU* hdu) const
269 311 {
270 312 // Continue only if HDU is valid
271 313 if (hdu != NULL) {
272 314  
273   - // Set coordinate system dependent strings
274   - std::string c_equinox;
275   - std::string c_crval1;
276   - std::string c_crval2;
277   - std::string c_cdelt1;
278   - std::string c_cdelt2;
279   - if (m_coordsys == 0) {
280   - c_equinox = "Equinox for celestial coordinate system";
281   - c_crval1 = "[deg] Right Ascension at reference pixel";
282   - c_crval2 = "[deg] Declination at reference pixel";
283   - c_cdelt1 = "[deg/pixel] Right Ascension increment at"
284   - " reference pixel";
285   - c_cdelt2 = "[deg/pixel] Declination increment at"
286   - " reference pixel";
  315 + // Write reference pixel coordinates
  316 + for (int i = 0; i < m_naxis; ++i) {
  317 + std::string keyname = "CRPIX"+str(i+1);
  318 + std::string comment = "Pixel coordinate of reference point (starting from 1)";
  319 + hdu->card(keyname, m_crpix.at(i), comment);
287 320 }
288   - else if (m_coordsys == 1) {
289   - c_equinox = "Equinox for celestial coordinate system";
290   - c_crval1 = "[deg] Galactic Longitude at reference pixel";
291   - c_crval2 = "[deg] Galactic Latitude at reference pixel";
292   - c_cdelt1 = "[deg/pixel] Galactic Longitude increment at"
293   - " reference pixel";
294   - c_cdelt2 = "[deg/pixel] Galactic Latitude increment at"
295   - " reference pixel";
  321 +
  322 + //TODO: Write linear transformation matrix
  323 +
  324 + // Write coordinate increment at reference point
  325 + for (int i = 0; i < m_naxis; ++i) {
  326 + std::string keyname = "CDELT"+str(i+1);
  327 + std::string comment;
  328 + if (m_cunit.at(i).length() > 0)
  329 + comment += "["+strip_whitespace(m_cunit.at(i))+"] ";
  330 + comment += "Coordinate increment at reference point";
  331 + hdu->card(keyname, m_cdelt.at(i), comment);
296 332 }
297   - else
298   - throw GException::wcs_bad_coords(G_WRITE, coordsys());
299   -
300   - // Set keywords
301   - hdu->card("EQUINOX", 2000.0, c_equinox);
302   - hdu->card("CTYPE1", wcs_crval1(), "Projection Type");
303   - hdu->card("CTYPE2", wcs_crval2(), "Projection Type");
304   - hdu->card("CDELT1", m_cdelt[0], c_cdelt1);
305   - hdu->card("CDELT2", m_cdelt[1], c_cdelt2);
306   - hdu->card("CRPIX1", m_crpix[0],
307   - "X index of reference pixel (starting from 1)");
308   - hdu->card("CRPIX2", m_crpix[1],
309   - "Y index of reference pixel (starting from 1)");
310   - hdu->card("CRVAL1", m_crval[0], c_crval1);
311   - hdu->card("CRVAL2", m_crval[1], c_crval2);
312   - hdu->card("CROTA2", 0.0, "[deg] Rotation Angle");
313   - hdu->card("LONPOLE", m_lonpole, "[deg] Native longitude of celestial pole");
314   - hdu->card("LATPOLE", m_latpole, "[deg] Native latitude of celestial pole");
  333 +
  334 + // Write units of coordinate increment and reference value
  335 + for (int i = 0; i < m_naxis; ++i) {
  336 + if (m_cunit.at(i).length() > 0) {
  337 + std::string keyname = "CUNIT"+str(i+1);
  338 + std::string comment = "Units of coordinate increment and value";
  339 + hdu->card(keyname, m_cunit.at(0), comment);
  340 + }
  341 + }
  342 +
  343 + // Write coordinate type
  344 + wcs_set_ctype();
  345 + for (int i = 0; i < m_naxis; ++i) {
  346 + if (i == m_lng) {
  347 + std::string keyname = "CTYPE"+str(i+1);
  348 + hdu->card(keyname, m_ctype.at(i), m_ctype_c.at(i));
  349 + }
  350 + if (i == m_lat) {
  351 + std::string keyname = "CTYPE"+str(i+1);
  352 + hdu->card(keyname, m_ctype.at(i), m_ctype_c.at(i));
  353 + }
  354 + if (i == m_spec) {
  355 + std::string keyname = "CTYPE"+str(i+1);
  356 + hdu->card(keyname, m_ctype.at(i), m_ctype_c.at(i));
  357 + }
  358 + }
  359 +
  360 + // Write coordinate value at reference point
  361 + for (int i = 0; i < m_naxis; ++i) {
  362 + std::string keyname = "CRVAL"+str(i+1);
  363 + std::string comment;
  364 + if (m_cunit.at(i).length() > 0)
  365 + comment += "["+strip_whitespace(m_cunit.at(i))+"] ";
  366 + comment += "Coordinate value at reference point";
  367 + hdu->card(keyname, m_crval.at(i), comment);
  368 + }
  369 +
  370 + //TODO: Parameter values
  371 + hdu->card("CROTA2", 0.0, "[deg] Rotation Angle"); // Old style, use PV instead
315 372 //hdu->card("PV2_1", 0.0, "Projection parameter 1");
316 373 //hdu->card("PV2_2", 0.0, "Projection parameter 2");
  374 +
  375 + // Celestial and spectral transformation parameters
  376 + if (!undefined(m_lonpole))
  377 + hdu->card("LONPOLE", m_lonpole, "[deg] Native longitude of celestial pole");
  378 + if (!undefined(m_latpole))
  379 + hdu->card("LATPOLE", m_latpole, "[deg] Native latitude of celestial pole");
  380 + if (!undefined(m_restfrq))
  381 + hdu->card("RESTFRQ", m_restfrq, "[Hz] Line rest frequency");
  382 + if (!undefined(m_restwav))
  383 + hdu->card("RESTWAV", m_restwav, "[Hz] Line rest wavelength");
  384 +
  385 + // Equatorial coordinate system type
  386 + if (m_radesys.length() > 0)
  387 + hdu->card("RADESYS", m_radesys, "Equatorial coordinate system");
  388 +
  389 + // Equinox of equatorial coordinate system
  390 + if (!undefined(m_equinox))
  391 + hdu->card("EQUINOX", m_equinox, "[yr] Equinox of equatorial coordinates");
  392 +
  393 + //TODO: Reference frame of spectral coordinates
  394 +
  395 + //TODO: Reference frame of spectral observation
  396 +
  397 + //TODO: Observer's velocity towards source
  398 +
  399 + //TODO: Reference frame of source redshift
  400 +
  401 + //TODO: Redshift of the source
  402 +
  403 + //TODO: Observatory coordinates
  404 +
  405 + //TODO: MJD of observation
  406 +
  407 + //TODO: MJD mid-observation time
  408 +
  409 + //TODO: ISO-8601 date corresponding to MJD-OBS
  410 +
  411 + //TODO: ISO-8601 date corresponding to MJD-AVG
317 412  
318 413 } // endif: HDU was valid
319 414  
... ... @@ -325,14 +420,40 @@ void GWcslib::write(GFitsHDU* hdu) const
325 420 /***********************************************************************//**
326 421 * @brief Returns solid angle of pixel in units of steradians
327 422 *
  423 + * @param[in] pix Pixel index.
  424 + *
  425 + * This dummy method throws an error when called as pixel indices are not
  426 + * implemented for a HealPix grid. Maybe we should drop this method and
  427 + * implement a general GSkyPixel class tha also handles HealPix?
  428 + *
  429 + * @todo Think about implementation of sky pixel for HealPix data. We may
  430 + * use only the first argument, and even carry a usage flag in
  431 + * GSkyPixel
  432 + ***************************************************************************/
  433 +double GWcslib::omega(const int& pix) const
  434 +{
  435 + // Set error message
  436 + std::string message = "Method not defined for general WCS projection.";
  437 +
  438 + // Throw error
  439 + throw GException::wcs(G_OMEGA, message);
  440 +
  441 + // Return
  442 + return 0.0;
  443 +}
  444 +
  445 +
  446 +/***********************************************************************//**
  447 + * @brief Returns solid angle of pixel in units of steradians
  448 + *
328 449 * @param[in] pix Pixel index (x,y)
329 450 *
330 451 * Estimate solid angles of pixel by compuing the coordinates in the 4 pixel
331 452 * corners. The surface is computed using a cartesian approximation:
332 453 * a
333   - * 1-----------2 a+b
334   - * h / \ where A = h ---
335   - * 4---------------3 2
  454 + * 1-----5-----2 a+b
  455 + * / | h \ where A = h ---
  456 + * 4-------6-------3 2
336 457 * b
337 458 * This is a brute force technique that works sufficiently well for non-
338 459 * rotated sky maps. Something more intelligent should be implemented in
... ... @@ -340,16 +461,13 @@ void GWcslib::write(GFitsHDU* hdu) const
340 461 ***************************************************************************/
341 462 double GWcslib::omega(const GSkyPixel& pix) const
342 463 {
343   - // Bypass const correctness
344   - GWcslib* ptr = (GWcslib*)this;
345   -
346   - // Get the sky directions of the 4 corners
347   - GSkyDir dir1 = ptr->xy2dir(GSkyPixel(pix.x()-0.5, pix.y()-0.5));
348   - GSkyDir dir2 = ptr->xy2dir(GSkyPixel(pix.x()+0.5, pix.y()-0.5));
349   - GSkyDir dir3 = ptr->xy2dir(GSkyPixel(pix.x()+0.5, pix.y()+0.5));
350   - GSkyDir dir4 = ptr->xy2dir(GSkyPixel(pix.x()-0.5, pix.y()+0.5));
351   - GSkyDir dir5 = ptr->xy2dir(GSkyPixel(pix.x(), pix.y()-0.5));
352   - GSkyDir dir6 = ptr->xy2dir(GSkyPixel(pix.x(), pix.y()+0.5));
  464 + // Get the sky directions of the 6 points
  465 + GSkyDir dir1 = xy2dir(GSkyPixel(pix.x()-0.5, pix.y()-0.5));
  466 + GSkyDir dir2 = xy2dir(GSkyPixel(pix.x()+0.5, pix.y()-0.5));
  467 + GSkyDir dir3 = xy2dir(GSkyPixel(pix.x()+0.5, pix.y()+0.5));
  468 + GSkyDir dir4 = xy2dir(GSkyPixel(pix.x()-0.5, pix.y()+0.5));
  469 + GSkyDir dir5 = xy2dir(GSkyPixel(pix.x(), pix.y()-0.5));
  470 + GSkyDir dir6 = xy2dir(GSkyPixel(pix.x(), pix.y()+0.5));
353 471  
354 472 // Compute distances between sky directions
355 473 double a = dir1.dist(dir2);
... ... @@ -365,36 +483,65 @@ double GWcslib::omega(const GSkyPixel&amp; pix) const
365 483  
366 484  
367 485 /***********************************************************************//**
368   - * @brief Returns sky direction of pixel
  486 + * @brief Returns sky direction of pixel (DUMMY METHOD)
369 487 *
370   - * @param[in] pix Sky pixel.
  488 + * @param[in] pix Pixel number (0,1,...,m_num_pixels).
371 489 *
372   - * Note that pixel indices start from 0.
  490 + * @exception GException::wcs
  491 + * Method not defined for projection.
  492 + *
  493 + * This method should be overloaded by the appropriate method of the derived
  494 + * class. It thus should never get called.
373 495 ***************************************************************************/
374   -GSkyDir GWcslib::xy2dir(const GSkyPixel& pix) const
  496 +GSkyDir GWcslib::pix2dir(const int& pix) const
375 497 {
376   - // Set constants
377   - const int __permutation[2] = {1,0};
  498 + // Set error message
  499 + std::string message = "Method not defined for general WCS projection.";
  500 +
  501 + // Throw error
  502 + throw GException::wcs(G_PIX2DIR, message);
378 503  
379   - // Determine offset
380   - GVector offset = GVector(pix.x(), pix.y());
  504 + // Define direction
  505 + GSkyDir dir;
  506 +
  507 + // Return
  508 + return dir;
  509 +}
381 510  
382   - // Determine xy
383   - GVector xy = m_cd * (offset - m_refpix);
384 511  
385   - // Swap result if coordinates are reversed
386   - if (m_reverse)
387   - xy = perm(xy, __permutation);
  512 +/***********************************************************************//**
  513 + * @brief Returns pixel for sky direction (DUMMY METHOD)
  514 + *
  515 + * @param[in] dir Sky direction.
  516 + *
  517 + * @exception GException::wcs
  518 + * Method not defined for projection.
  519 + *
  520 + * This method should be overloaded by the appropriate method of the derived
  521 + * class. It thus should never get called.
  522 + ***************************************************************************/
  523 +int GWcslib::dir2pix(const GSkyDir& dir) const
  524 +{
  525 + // Set error message
  526 + std::string message = "Method not defined for general WCS projection.";
388 527  
389   - // Perform map projection
390   - GVector native = xy;
391   - nat2std(&native);
  528 + // Throw error
  529 + throw GException::wcs(G_DIR2PIX, message);
392 530  
393   - // Get sky direction for native coordinates
394   - //GSkyDir dir = wcs_native2dir(native);
395   -
396   - // NEW METHOD -- DELETE ALL OF ABOVE
  531 + // Return
  532 + return 0;
  533 +}
397 534  
  535 +
  536 +/***********************************************************************//**
  537 + * @brief Returns sky direction of pixel
  538 + *
  539 + * @param[in] pix Sky pixel.
  540 + *
  541 + * Note that pixel indices start from 0.
  542 + ***************************************************************************/
  543 +GSkyDir GWcslib::xy2dir(const GSkyPixel& pix) const
  544 +{
398 545 // Allocate memory for transformation
399 546 double pixcrd[2];
400 547 double imgcrd[2];
... ... @@ -420,8 +567,7 @@ GSkyDir GWcslib::xy2dir(const GSkyPixel&amp; pix) const
420 567  
421 568 // Debug: Dump transformation steps
422 569 #if defined(G_XY2DIR_DEBUG)
423   - std::cout << "xy2dir: pixel=" << offset << " xy=" << xy
424   - << " native=" << native << " dir=" << wcs_native2dir(native)
  570 + std::cout << "xy2dir: pixel=" << pix
425 571 << " (x,y)=(" << pixcrd[0] << "," << pixcrd[1] << ")"
426 572 << " (phi,theta)=(" << phi << "," << theta << ")"
427 573 << " (lng,lat)=(" << world[0] << "," << world[1] << ")" << std::endl;
... ... @@ -440,30 +586,8 @@ GSkyDir GWcslib::xy2dir(const GSkyPixel&amp; pix) const
440 586 * Note that GSkyPixel indices start from 0 while the WCS pixel reference
441 587 * starts from 1.
442 588 ***************************************************************************/
443   -GSkyPixel GWcslib::dir2xy(GSkyDir dir) const
  589 +GSkyPixel GWcslib::dir2xy(const GSkyDir& dir) const
444 590 {
445   - // Set constants
446   - const int __permutation[2] = {1,0};
447   -
448   - // Get native coordinates for sky direction in radians
449   - GVector native = wcs_dir2native(dir);
450   -
451   - // Perform map projection
452   - GVector xy = native;
453   - std2nat(&xy);
454   -
455   - // Swap result if coordinates are reversed
456   - if (m_reverse)
457   - xy = perm(xy, __permutation);
458   -
459   - // Determine pixel offset
460   - GVector offset = m_invcd * xy + m_refpix;
461   -
462   - // Set sky pixel
463   - //GSkyPixel pixel(offset[0], offset[1]);
464   -
465   - // NEW METHOD -- DELETE ALL OF ABOVE
466   -
467 591 // Allocate memory for transformation
468 592 double pixcrd[2];
469 593 double imgcrd[2];
... ... @@ -491,8 +615,7 @@ GSkyPixel GWcslib::dir2xy(GSkyDir dir) const
491 615  
492 616 // Debug: Dump transformation steps
493 617 #if defined(G_DIR2XY_DEBUG)
494   - std::cout << "dir2xy: dir=" << dir << " native=" << native
495   - << " xy=" << xy << " pixel=" << offset
  618 + std::cout << "dir2xy: dir=" << dir
496 619 << " (lng,lat)=(" << world[0] << "," << world[1] << ")"
497 620 << " (phi,theta)=(" << phi << "," << theta << ")"
498 621 << " (x,y)=" << pixel << std::endl;
... ... @@ -503,6 +626,39 @@ GSkyPixel GWcslib::dir2xy(GSkyDir dir) const
503 626 }
504 627  
505 628  
  629 +/***********************************************************************//**
  630 + * @brief Set World Coordinate System parameters
  631 + *
  632 + * @param[in] coords Coordinate system.
  633 + * @param[in] crval1 X value of reference pixel.
  634 + * @param[in] crval2 Y value of reference pixel.
  635 + * @param[in] crpix1 X index of reference pixel (first pixel is 1).
  636 + * @param[in] crpix2 Y index of reference pixel (first pixel is 1).
  637 + * @param[in] cdelt1 Increment in x direction at reference pixel (deg).
  638 + * @param[in] cdelt2 Increment in y direction at reference pixel (deg).
  639 + *
  640 + * This method sets the WCS parameters.
  641 + ***************************************************************************/
  642 +void GWcslib::set(const std::string& coords,
  643 + const double& crval1, const double& crval2,
  644 + const double& crpix1, const double& crpix2,
  645 + const double& cdelt1, const double& cdelt2)
  646 +
  647 +{
  648 + // Clear any existing information
  649 + clear();
  650 +
  651 + // Set standard parameters
  652 + set_members(coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2);
  653 +
  654 + // Setup WCS derived parameters
  655 + wcs_set();
  656 +
  657 + // Return
  658 + return;
  659 +}
  660 +
  661 +
506 662 /*==========================================================================
507 663 = =
508 664 = Protected methods =
... ... @@ -534,14 +690,17 @@ void GWcslib::copy_members(const GWcslib&amp; wcs)
534 690 // Copy WCS attributes
535 691 m_wcsset = wcs.m_wcsset;
536 692 m_naxis = wcs.m_naxis;
537   - m_Crval = wcs.m_Crval;
  693 + m_crval = wcs.m_crval;
538 694 m_cunit = wcs.m_cunit;
539 695 m_ctype = wcs.m_ctype;
  696 + m_ctype_c = wcs.m_ctype_c;
540 697 m_lonpole = wcs.m_lonpole;
541 698 m_latpole = wcs.m_latpole;
542 699 m_restfrq = wcs.m_restfrq;
543 700 m_restwav = wcs.m_restwav;
544   - m_Cd = wcs.m_Cd;
  701 + m_radesys = wcs.m_radesys;
  702 + m_equinox = wcs.m_equinox;
  703 + m_cd = wcs.m_cd;
545 704 m_crota = wcs.m_crota;
546 705 m_lng = wcs.m_lng;
547 706 m_lat = wcs.m_lat;
... ... @@ -550,9 +709,9 @@ void GWcslib::copy_members(const GWcslib&amp; wcs)
550 709 // Copy linear transformation parameters
551 710 m_linset = wcs.m_linset;
552 711 m_unity = wcs.m_unity;
553   - m_Crpix = wcs.m_Crpix;
  712 + m_crpix = wcs.m_crpix;
554 713 m_pc = wcs.m_pc;
555   - m_Cdelt = wcs.m_Cdelt;
  714 + m_cdelt = wcs.m_cdelt;
556 715 m_piximg = wcs.m_piximg;
557 716 m_imgpix = wcs.m_imgpix;
558 717  
... ... @@ -575,12 +734,11 @@ void GWcslib::copy_members(const GWcslib&amp; wcs)
575 734  
576 735 // Copy projection parameters
577 736 m_prjset = wcs.m_prjset;
578   - m_code = wcs.m_code;
579   - m_name = wcs.m_name;
580 737 m_r0 = wcs.m_r0;
581 738 m_bounds = wcs.m_bounds;
582 739 m_x0 = wcs.m_x0;
583 740 m_y0 = wcs.m_y0;
  741 + m_w = wcs.m_w;
584 742 for (int i = 0; i < PVN; ++i)
585 743 m_pv[i] = wcs.m_pv[i];
586 744  
... ... @@ -600,7 +758,7 @@ void GWcslib::free_members(void)
600 758  
601 759  
602 760 /***********************************************************************//**
603   - * @brief Set standard WCS parameters
  761 + * @brief Set World Coordinate System parameters
604 762 *
605 763 * @param[in] coords Coordinate system.
606 764 * @param[in] crval1 X value of reference pixel.
... ... @@ -610,18 +768,11 @@ void GWcslib::free_members(void)
610 768 * @param[in] cdelt1 Increment in x direction at reference pixel (deg).
611 769 * @param[in] cdelt2 Increment in y direction at reference pixel (deg).
612 770 *
613   - * @exception GException::wcs
614   - * Unable to invert CD matrix.
615   - *
616   - * This method sets the WCS standard parameters as class members. It does
617   - * however not set the rotation matrices. A call to wcs_init is mandatory
618   - * before the projection can be used.
  771 + * This method sets the WCS parameters. It does not call wcs_set(), however,
  772 + * as wcs_set() may depend on the availability of derived class methods which
  773 + * can not be used in a base class constructor.
619 774 *
620   - * @todo Implement parameter check
621   - *
622   - * Note: This method will be called by the base class constructors, hence
623   - * it cannot contain any call to pure virtual functions!!! The method
624   - * should definitely be renamed in the new version!!!
  775 + * @todo Implement parameter validity check
625 776 ***************************************************************************/
626 777 void GWcslib::set_members(const std::string& coords,
627 778 const double& crval1, const double& crval2,
... ... @@ -631,58 +782,64 @@ void GWcslib::set_members(const std::string&amp; coords,
631 782 {
632 783 //TODO: Check parameters
633 784  
  785 + // Initialise WCS
  786 + wcs_ini(2);
  787 +
634 788 // Set coordinate system
635 789 coordsys(coords);
636 790  
637   - // Set parameters
638   - m_crval = GVector(crval1, crval2);
639   - m_crpix = GVector(crpix1, crpix2);
640   - m_cdelt = GVector(cdelt1, cdelt2);
641   - m_refval = m_crval;
642   - m_refpix = m_crpix - GVector(1.0,1.0); //!< Needed as skypixel starts with 0
643   -
644   - // Set CD matrix corresponding to no rotation
645   - m_cd(0,0) = m_cdelt[0];
646   - m_cd(1,0) = 0.0;
647   - m_cd(0,1) = 0.0;
648   - m_cd(1,1) = m_cdelt[1];
649   -
650   - // Compute inverse CD matrix
651   - double delta = m_cd(0,0)*m_cd(1,1) - m_cd(0,1)*m_cd(1,0);
652   - if (delta != 0.0) {
653   - m_invcd(0,0) = m_cd(1,1) / delta;
654   - m_invcd(0,1) = -m_cd(0,1) / delta;
655   - m_invcd(1,0) = -m_cd(1,0) / delta;
656   - m_invcd(1,1) = m_cd(0,0) / delta;
657   - }
658   - else
659   - throw GException::wcs(G_WCS_SET, "Unable to invert CD matrix.");
660   -
661   - // NEW VERSION
662   -
663 791 // Set World Coordinate parameters
664   - wcs_ini(2);
665 792 m_lng = 0;
666 793 m_lat = 1;
667   - m_Crpix[0] = crpix1;
668   - m_Crpix[1] = crpix2;
669   - m_Cdelt[0] = cdelt1;
670   - m_Cdelt[1] = cdelt2;
671   - m_Crval[0] = crval1;
672   - m_Crval[1] = crval2;
673   - //m_Cd.push_back(cdelt1);
674   - //m_Cd.push_back(0.0);
675   - //m_Cd.push_back(0.0);
676   - //m_Cd.push_back(cdelt2);
  794 + m_crpix[0] = crpix1;
  795 + m_crpix[1] = crpix2;
  796 + m_cdelt[0] = cdelt1;
  797 + m_cdelt[1] = cdelt2;
  798 + m_crval[0] = crval1;
  799 + m_crval[1] = crval2;
677 800  
678 801 // Return
679 802 return;
680 803 }
681 804  
682 805  
  806 +/***********************************************************************//**
  807 + * @brief Returns true if argument is identical
  808 + *
  809 + * @param[in] wcs Pointer to World Coordinate System
  810 + *
  811 + * This method is a helper for the World Coordinate Comparison friends. It
  812 + * does not compare derived quantities as those may not have been initialised
  813 + * in both objects.
  814 + ***************************************************************************/
  815 +bool GWcslib::compare(const GWcs& wcs) const
  816 +{
  817 + // Initialise result
  818 + bool result = false;
  819 +
  820 + // Continue only we compare to a GWcslib object
  821 + const GWcslib* ptr = dynamic_cast<const GWcslib*>(&wcs);
  822 + if (ptr != NULL) {
  823 +
  824 + // Perform comparion of standard (non derived) parameters
  825 + result = ((code() == ptr->code()) &&
  826 + (m_coordsys == ptr->m_coordsys) &&
  827 + (m_lng == ptr->m_lng) &&
  828 + (m_lat == ptr->m_lat) &&
  829 + (m_crval == ptr->m_crval) &&
  830 + (m_crpix == ptr->m_crpix) &&
  831 + (m_cdelt == ptr->m_cdelt));
  832 +
  833 + }
  834 +
  835 + // Return result
  836 + return result;
  837 +}
  838 +
  839 +
683 840 /*==========================================================================
684 841 = =
685   - = World Coordinate Methods adapted from wcslib =
  842 + = World Coordinate Methods adapted from wcslib =
686 843 = =
687 844 ==========================================================================*/
688 845  
... ... @@ -703,13 +860,14 @@ void GWcslib::wcs_ini(int naxis)
703 860 m_wcsset = false;
704 861  
705 862 // Clear vectors
706   - m_Crpix.clear();
  863 + m_crpix.clear();
707 864 m_pc.clear();
708   - m_Cdelt.clear();
709   - m_Crval.clear();
  865 + m_cdelt.clear();
  866 + m_crval.clear();
710 867 m_cunit.clear();
711 868 m_ctype.clear();
712   - m_Cd.clear();
  869 + m_ctype_c.clear();
  870 + m_cd.clear();
713 871 m_crota.clear();
714 872  
715 873 // Set number of axes
... ... @@ -721,9 +879,10 @@ void GWcslib::wcs_ini(int naxis)
721 879 // Set default axes parameters
722 880 if (naxis > 0) {
723 881 for (int i = 0; i < naxis; ++i) {
724   - m_Crval.push_back(0.0);
  882 + m_crval.push_back(0.0);
725 883 m_cunit.push_back("");
726 884 m_ctype.push_back("");
  885 + m_ctype_c.push_back("");
727 886 }
728 887 }
729 888  
... ... @@ -732,8 +891,8 @@ void GWcslib::wcs_ini(int naxis)
732 891 m_latpole = +90.0;
733 892  
734 893 // Set defaults for the spectral transformation parameters
735   - m_restfrq = 0.0;
736   - m_restwav = 0.0;
  894 + m_restfrq = UNDEFINED; // 0.0 in wcslib
  895 + m_restwav = UNDEFINED; // 0.0 in wcslib
737 896  
738 897 //TODO: Default parameter values
739 898  
... ... @@ -741,12 +900,14 @@ void GWcslib::wcs_ini(int naxis)
741 900 if (naxis > 0) {
742 901 for (int i = 0; i < naxis; ++i) {
743 902 for (int j = 0; j < naxis; ++j)
744   - m_Cd.push_back(0.0);
  903 + m_cd.push_back(0.0);
745 904 m_crota.push_back(0.0);
746 905 }
747 906 }
748 907  
749 908 //TODO: Defaults for auxiliary coordinate system information
  909 + m_radesys.clear();
  910 + m_equinox = UNDEFINED;
750 911  
751 912 // Reset derived values
752 913 m_lng = -1;
... ... @@ -796,8 +957,8 @@ void GWcslib::wcs_set(void) const
796 957 cel_ini();
797 958  
798 959 // Set CRVALia, LONPOLEa, and LATPOLEa keyvalues
799   - m_ref[0] = m_Crval[m_lng]; // Longitude reference value
800   - m_ref[1] = m_Crval[m_lat]; // Latitude reference value
  960 + m_ref[0] = m_crval[m_lng]; // Longitude reference value
  961 + m_ref[1] = m_crval[m_lat]; // Latitude reference value
801 962 m_ref[2] = m_lonpole; // LONPOLE
802 963 m_ref[3] = m_latpole; // LATPOLE
803 964  
... ... @@ -833,6 +994,98 @@ void GWcslib::wcs_set(void) const
833 994  
834 995  
835 996 /***********************************************************************//**
  997 + * @brief Set CTYPEa keywords
  998 + *
  999 + * @exception GException::wcs_invalid
  1000 + * WCS projection not valid.
  1001 + * @exception GException::wcs_bad_coords
  1002 + * Coordinate system is not of valid type.
  1003 + *
  1004 + * This method has been inspired by code from wcshdr.c::wcshdo.
  1005 + ***************************************************************************/
  1006 +void GWcslib::wcs_set_ctype(void) const
  1007 +{
  1008 + // Check on correct type length
  1009 + if (code().length() != 3)
  1010 + throw GException::wcs_invalid(G_WCS_SET_CTYPE, code(),
  1011 + "3-character type required.");
  1012 +
  1013 + // Set longitude keyword
  1014 + if (m_lng >= 0) {
  1015 +
  1016 + // Set coordinate system
  1017 + if (m_coordsys == 0) {
  1018 + m_ctype.at(m_lng) = "RA---" + code();
  1019 + m_ctype_c.at(m_lng) = "Right ascension, ";
  1020 + }
  1021 + else if (m_coordsys == 1) {
  1022 + m_ctype.at(m_lng) = "GLON-" + code();
  1023 + m_ctype_c.at(m_lng) = "Galactic longitude, ";
  1024 + }
  1025 + else if (m_coordsys == 2) {
  1026 + m_ctype.at(m_lng) = "ELON-" + code();
  1027 + m_ctype_c.at(m_lng) = "Ecliptic longitude, ";
  1028 + }
  1029 + else if (m_coordsys == 3) {
  1030 + m_ctype.at(m_lng) = "HLON-" + code();
  1031 + m_ctype_c.at(m_lng) = "Helioecliptic longitude, ";
  1032 + }
  1033 + else if (m_coordsys == 4) {
  1034 + m_ctype.at(m_lng) = "SLON-" + code();
  1035 + m_ctype_c.at(m_lng) = "Supergalactic longitude, ";
  1036 + }
  1037 + else
  1038 + throw GException::wcs_bad_coords(G_WCS_SET_CTYPE, coordsys());
  1039 +
  1040 + // Add projection name to comment
  1041 + m_ctype_c.at(m_lng).append(name());
  1042 + m_ctype_c.at(m_lng).append(" projection");
  1043 + }
  1044 +
  1045 + // Set latitude keyword
  1046 + if (m_lat >= 0) {
  1047 +
  1048 + // Set coordinate system
  1049 + if (m_coordsys == 0) {
  1050 + m_ctype.at(m_lat) = "DEC--" + code();
  1051 + m_ctype_c.at(m_lat) = "Declination, ";
  1052 + }
  1053 + else if (m_coordsys == 1) {
  1054 + m_ctype.at(m_lat) = "GLAT-" + code();
  1055 + m_ctype_c.at(m_lat) = "Galactic latitude, ";
  1056 + }
  1057 + else if (m_coordsys == 2) {
  1058 + m_ctype.at(m_lat) = "ELAT-" + code();
  1059 + m_ctype_c.at(m_lat) = "Ecliptic latitude, ";
  1060 + }
  1061 + else if (m_coordsys == 3) {
  1062 + m_ctype.at(m_lat) = "HLAT-" + code();
  1063 + m_ctype_c.at(m_lat) = "Helioecliptic latitude, ";
  1064 + }
  1065 + else if (m_coordsys == 4) {
  1066 + m_ctype.at(m_lat) = "SLAT-" + code();
  1067 + m_ctype_c.at(m_lat) = "Supergalactic latitude, ";
  1068 + }
  1069 + else
  1070 + throw GException::wcs_bad_coords(G_WCS_SET_CTYPE, coordsys());
  1071 +
  1072 + // Add projection name to comment
  1073 + m_ctype_c.at(m_lng).append(name());
  1074 + m_ctype_c.at(m_lng).append(" projection");
  1075 + }
  1076 +
  1077 + //TODO: Set spectral keyword
  1078 + if (m_spec >= 0) {
  1079 + m_ctype.at(m_spec) = "NONE";
  1080 + m_ctype_c.at(m_spec) = "Not yet implemented";
  1081 + }
  1082 +
  1083 + // Return
  1084 + return;
  1085 +}
  1086 +
  1087 +
  1088 +/***********************************************************************//**
836 1089 * @brief Pixel-to-world transformation
837 1090 *
838 1091 * @param[in] ncoord Number of coordinates.
... ... @@ -967,33 +1220,6 @@ std::string GWcslib::wcs_print(void) const
967 1220 // Initialise result string
968 1221 std::string result;
969 1222  
970   - // Append information
971   - GVector native = m_native_pole*rad2deg;
972   - result.append(parformat("Coordinate system")+coordsys()+"\n");
973   - result.append(parformat("Reference coordinate")+m_crval.print()+" deg\n");
974   - result.append(parformat("Reference pixel")+m_crpix.print()+"\n");
975   - result.append(parformat("Increment at reference")+m_cdelt.print()+" deg\n");
976   - result.append(parformat("Coordinate of North Pole")+m_npole.print()+" deg\n");
977   -
978   - // Append CD matrix
979   - result.append(parformat("CD matrix"));
980   - result.append("["+str(m_cd(0,0)));
981   - result.append(" "+str(m_cd(0,1)));
982   - result.append("]["+str(m_cd(1,0)));
983   - result.append(" "+str(m_cd(1,1))+"]\n");
984   -
985   - // Append inverse CD matrix
986   - result.append(parformat("Inverse CD matrix"));
987   - result.append("["+str(m_invcd(0,0)));
988   - result.append(" "+str(m_invcd(0,1)));
989   - result.append("]["+str(m_invcd(1,0)));
990   - result.append(" "+str(m_invcd(1,1))+"]\n");
991   -
992   - // Append coordinate of native pole
993   - result.append(parformat("Coordinate of native pole")+native.print()+"\n");
994   -
995   - // NEW VERSION
996   -
997 1223 // Append World Coordinate parameters
998 1224 result.append(parformat("Number of axes")+str(m_naxis)+"\n");
999 1225 result.append(parformat("Longitude axis")+str(m_lng)+"\n");
... ... @@ -1002,27 +1228,27 @@ std::string GWcslib::wcs_print(void) const
1002 1228  
1003 1229 // Append coordinates
1004 1230 result.append(parformat("Reference coordinate")+"(");
1005   - for (int i = 0; i < m_Crval.size(); ++i) {
  1231 + for (int i = 0; i < m_crval.size(); ++i) {
1006 1232 if (i > 0)
1007 1233 result.append(", ");
1008   - result.append(str(m_Crval[i]));
  1234 + result.append(str(m_crval[i]));
1009 1235 if (m_cunit[i].length() > 0)
1010 1236 result.append(" "+m_cunit[i]);
1011 1237  
1012 1238 }
1013 1239 result.append(")\n");
1014 1240 result.append(parformat("Reference pixel")+"(");
1015   - for (int i = 0; i < m_Crpix.size(); ++i) {
  1241 + for (int i = 0; i < m_crpix.size(); ++i) {
1016 1242 if (i > 0)
1017 1243 result.append(", ");
1018   - result.append(str(m_Crpix[i]));
  1244 + result.append(str(m_crpix[i]));
1019 1245 }
1020 1246 result.append(")\n");
1021 1247 result.append(parformat("Increment at reference")+"(");
1022   - for (int i = 0; i < m_Cdelt.size(); ++i) {
  1248 + for (int i = 0; i < m_cdelt.size(); ++i) {
1023 1249 if (i > 0)
1024 1250 result.append(", ");
1025   - result.append(str(m_Cdelt[i]));
  1251 + result.append(str(m_cdelt[i]));
1026 1252 if (m_cunit[i].length() > 0)
1027 1253 result.append(" "+m_cunit[i]);
1028 1254 }
... ... @@ -1030,22 +1256,13 @@ std::string GWcslib::wcs_print(void) const
1030 1256  
1031 1257 // Append origin
1032 1258 result.append(parformat("(Phi_0, Theta_0)")+"(");
1033   - if (undefined(m_phi0))
1034   - result.append("UNDEFINED, ");
1035   - else
1036   - result.append(str(m_phi0)+", ");
1037   - if (undefined(m_theta0))
1038   - result.append("UNDEFINED)\n");
1039   - else
1040   - result.append(str(m_theta0)+") deg\n");
  1259 + result.append(wcs_print_value(m_phi0)+", ");
  1260 + result.append(wcs_print_value(m_theta0)+") deg\n");
1041 1261  
1042 1262 // Append native pole
1043 1263 result.append(parformat("(Phi_p, Theta_p)")+"(");
1044   - if (undefined(m_lonpole))
1045   - result.append("UNDEFINED, ");
1046   - else
1047   - result.append(str(m_lonpole)+", ");
1048   - result.append(str(m_latpole)+") deg\n");
  1264 + result.append(wcs_print_value(m_lonpole)+", ");
  1265 + result.append(wcs_print_value(m_latpole)+") deg\n");
1049 1266  
1050 1267 // Append LATPOLEa keyword usage
1051 1268 result.append(parformat("LATPOLE keyword usage"));
... ... @@ -1061,6 +1278,9 @@ std::string GWcslib::wcs_print(void) const
1061 1278 case 2:
1062 1279 result.append("Theta_p was specified solely by LATPOLE.\n");
1063 1280 break;
  1281 + default:
  1282 + result.append("UNDEFINED\n");
  1283 + break;
1064 1284 }
1065 1285  
1066 1286 // Append celestial transformation parameters
... ... @@ -1068,10 +1288,7 @@ std::string GWcslib::wcs_print(void) const
1068 1288 for (int k = 0; k < 4; ++k) {
1069 1289 if (k > 0)
1070 1290 result.append(", ");
1071   - if (undefined(m_ref[k]))
1072   - result.append("UNDEFINED");
1073   - else
1074   - result.append(str(m_ref[k]));
  1291 + result.append(wcs_print_value(m_ref[k]));
1075 1292 }
1076 1293 result.append(") deg\n");
1077 1294  
... ... @@ -1112,9 +1329,12 @@ std::string GWcslib::wcs_print(void) const
1112 1329 }
1113 1330 result.append(")\n");
1114 1331  
  1332 + // Append coordinate system
  1333 + result.append(parformat("Coodinate system")+coordsys()+"\n");
  1334 +
1115 1335 // Append projection parameters
1116   - result.append(parformat("Projection code")+m_code+"\n");
1117   - result.append(parformat("Projection name")+m_name+"\n");
  1336 + result.append(parformat("Projection code")+code()+"\n");
  1337 + result.append(parformat("Projection name")+name()+"\n");
1118 1338 result.append(parformat("Radius of the gen. sphere")+str(m_r0)+" deg\n");
1119 1339  
1120 1340 // Append boundary checking information
... ... @@ -1131,14 +1351,38 @@ std::string GWcslib::wcs_print(void) const
1131 1351 result.append(parformat("Fiducial offset")+"("+str(m_x0)+", "+str(m_y0)+")\n");
1132 1352  
1133 1353 // Append spectral transformation parameters
1134   - result.append(parformat("Rest frequency")+str(m_restfrq)+"\n");
1135   - result.append(parformat("Rest wavelength")+str(m_restwav)+"\n");
  1354 + result.append(parformat("Rest frequency")+wcs_print_value(m_restfrq)+"\n");
  1355 + result.append(parformat("Rest wavelength")+wcs_print_value(m_restwav)+"\n");
1136 1356  
1137 1357 // Return
1138 1358 return result;
1139 1359 }
1140 1360  
1141 1361  
  1362 +/***********************************************************************//**
  1363 + * @brief Helper function for value printing
  1364 + *
  1365 + * @param[in] value Double precision value
  1366 + *
  1367 + * This helper function either prints a double precision value or the
  1368 + * word 'UNDEFINED', depending on whether the value is defined or not.
  1369 + ***************************************************************************/
  1370 +std::string GWcslib::wcs_print_value(const double& value) const
  1371 +{
  1372 + // Initialise result
  1373 + std::string result;
  1374 +
  1375 + // Depend value dependent string
  1376 + if (undefined(value))
  1377 + result.append("UNDEFINED");
  1378 + else
  1379 + result.append(str(value));
  1380 +
  1381 + // Return result
  1382 + return result;
  1383 +}
  1384 +
  1385 +
1142 1386 /*==========================================================================
1143 1387 = =
1144 1388 = Celestial transformation methods adapted from wcslib =
... ... @@ -1978,9 +2222,9 @@ void GWcslib::lin_ini(int naxis)
1978 2222 m_unity = false;
1979 2223  
1980 2224 // Clear vectors
1981   - m_Crpix.clear();
  2225 + m_crpix.clear();
1982 2226 m_pc.clear();
1983   - m_Cdelt.clear();
  2227 + m_cdelt.clear();
1984 2228 m_piximg.clear();
1985 2229 m_imgpix.clear();
1986 2230  
... ... @@ -1991,7 +2235,7 @@ void GWcslib::lin_ini(int naxis)
1991 2235 for (int i = 0; i < naxis; ++i) {
1992 2236  
1993 2237 // CRPIXja defaults to 0.0
1994   - m_Crpix.push_back(0.0);
  2238 + m_crpix.push_back(0.0);
1995 2239  
1996 2240 // PCi_ja defaults to the unit matrix
1997 2241 for (int j = 0; j < naxis; ++j) {
... ... @@ -2002,7 +2246,7 @@ void GWcslib::lin_ini(int naxis)
2002 2246 }
2003 2247  
2004 2248 // CDELTia defaults to 1.0
2005   - m_Cdelt.push_back(1.0);
  2249 + m_cdelt.push_back(1.0);
2006 2250  
2007 2251 } // endfor: looped over axes
2008 2252  
... ... @@ -2058,7 +2302,7 @@ void GWcslib::lin_set(void) const
2058 2302 // Compute the pixel-to-image transformation matrix
2059 2303 for (int i = 0, index = 0; i < m_naxis; ++i) {
2060 2304 for (int j = 0; j < m_naxis; ++j, ++index)
2061   - m_piximg.push_back(m_Cdelt[i] * m_pc[index]);
  2305 + m_piximg.push_back(m_cdelt[i] * m_pc[index]);
2062 2306 }
2063 2307  
2064 2308 // Compute the image-to-pixel transformation matrix
... ... @@ -2105,7 +2349,7 @@ void GWcslib::lin_p2x(int ncoord, int nelem, const double* pixcrd, double* imgcr
2105 2349  
2106 2350 // Transform the first m_naxis elements
2107 2351 for (int i = 0; i < m_naxis; ++i)
2108   - *(img++) = m_Cdelt[i] * (*(pix++) - m_Crpix[i]);
  2352 + *(img++) = m_cdelt[i] * (*(pix++) - m_crpix[i]);
2109 2353  
2110 2354 // Go to next coordinate
2111 2355 pix += (nelem - m_naxis);
... ... @@ -2128,7 +2372,7 @@ void GWcslib::lin_p2x(int ncoord, int nelem, const double* pixcrd, double* imgcr
2128 2372 // Perform matrix multiplication (column-wise multiplication
2129 2373 // allows this to be cached)
2130 2374 for (int j = 0; j < m_naxis; ++j) {
2131   - double temp = *(pix++) - m_Crpix[j];
  2375 + double temp = *(pix++) - m_crpix[j];
2132 2376 for (int i = 0, ji = j; i < m_naxis; ++i, ji += m_naxis)
2133 2377 img[i] += m_piximg[ji] * temp;
2134 2378 }
... ... @@ -2177,7 +2421,7 @@ void GWcslib::lin_x2p(int ncoord, int nelem, const double* imgcrd, double* pixcr
2177 2421  
2178 2422 // Transform the first m_naxis elements
2179 2423 for (int i = 0; i < m_naxis; ++i)
2180   - *(pix++) = (*(img++) / m_Cdelt[i]) + m_Crpix[i];
  2424 + *(pix++) = (*(img++) / m_cdelt[i]) + m_crpix[i];
2181 2425  
2182 2426 // Go to next coordinate
2183 2427 pix += (nelem - m_naxis);
... ... @@ -2198,7 +2442,7 @@ void GWcslib::lin_x2p(int ncoord, int nelem, const double* imgcrd, double* pixcr
2198 2442 *pix = 0.0;
2199 2443 for (int i = 0; i < m_naxis; ++i, ++ji)
2200 2444 *pix += m_imgpix[ji] * img[i];
2201   - *(pix++) += m_Crpix[j];
  2445 + *(pix++) += m_crpix[j];
2202 2446 }
2203 2447  
2204 2448 // Go to next coordinate
... ... @@ -2372,8 +2616,6 @@ void GWcslib::prj_ini(void) const
2372 2616 {
2373 2617 // Initialise parameters
2374 2618 m_prjset = false;
2375   - m_code = "";
2376   - m_name = "";
2377 2619 m_r0 = 0.0;
2378 2620 m_pv[0] = 0.0;
2379 2621 m_pv[1] = UNDEFINED;
... ... @@ -2384,6 +2626,7 @@ void GWcslib::prj_ini(void) const
2384 2626 m_bounds = true; // 1 in wcslib
2385 2627 m_x0 = 0.0;
2386 2628 m_y0 = 0.0;
  2629 + m_w.clear();
2387 2630  
2388 2631 // Return
2389 2632 return;
... ... @@ -2669,70 +2912,3 @@ void GWcslib::sincosd(const double&amp; angle, double *s, double *c) const
2669 2912 = Friends =
2670 2913 = =
2671 2914 ==========================================================================*/
2672   -
2673   -/***********************************************************************//**
2674   - * @brief Equality operator
2675   - *
2676   - * @param[in] a First WCS.
2677   - * @param[in] b Second WCS.
2678   - ***************************************************************************/
2679   -bool operator== (const GWcslib &a, const GWcslib &b)
2680   -{
2681   - // Set result
2682   - bool result = ((a.m_type == b.m_type) &&
2683   - (a.m_coordsys == b.m_coordsys) &&
2684   - (a.m_reverse == b.m_reverse) &&
2685   - (a.m_crval == b.m_crval) &&
2686   - (a.m_crpix == b.m_crpix) &&
2687   - (a.m_cdelt == b.m_cdelt) &&
2688   - (a.m_npole == b.m_npole) &&
2689   - (a.m_cd == b.m_cd));
2690   -
2691   - // Return result
2692   - return result;
2693   -}
2694   -
2695   -
2696   -/***********************************************************************//**
2697   - * @brief Non-equality operator
2698   - *
2699   - * @param[in] a First WCS.
2700   - * @param[in] b Second WCS.
2701   - ***************************************************************************/
2702   -bool operator!= (const GWcslib &a, const GWcslib &b)
2703   -{
2704   - // Return result
2705   - return !(a == b);
2706   -}
2707   -
2708   -
2709   -/***********************************************************************//**
2710   - * @brief Output operator
2711   - *
2712   - * @param[in] os Output stream.
2713   - * @param[in] wcs WCS.
2714   - ***************************************************************************/
2715   -std::ostream& operator<< (std::ostream& os, const GWcslib& wcs)
2716   -{
2717   - // Write WCS in output stream
2718   - os << wcs.print();
2719   -
2720   - // Return output stream
2721   - return os;
2722   -}
2723   -
2724   -
2725   -/***********************************************************************//**
2726   - * @brief Log operator
2727   - *
2728   - * @param[in] log Logger.
2729   - * @param[in] wcs WCS.
2730   - ***************************************************************************/
2731   -GLog& operator<< (GLog& log, const GWcslib& wcs)
2732   -{
2733   - // Write WCS into logger
2734   - log << wcs.print();
2735   -
2736   - // Return logger
2737   - return log;
2738   -}
... ...
test/test_GSky.cpp
  View file @5db80af
... ... @@ -105,8 +105,6 @@ void test_GWcsCAR(void)
105 105 double crpix2 = 100.5;
106 106 double cdelt1 = 0.02;
107 107 double cdelt2 = 0.02;
108   - GMatrix cd(2,2);
109   - GVector pv2(21);
110 108  
111 109 // Test void constructor
112 110 try {
... ... @@ -123,8 +121,7 @@ void test_GWcsCAR(void)
123 121  
124 122 // Test constructor
125 123 try {
126   - GWcsCAR wcs(coords, crval1, crval2, crpix1, crpix2,
127   - cdelt1, cdelt2, cd, pv2);
  124 + GWcsCAR wcs(coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2);
128 125 std::cout << wcs << std::endl;
129 126 }
130 127 catch (std::exception &e) {
... ... @@ -137,8 +134,7 @@ void test_GWcsCAR(void)
137 134 std::cout << ".";
138 135  
139 136 // Build now skymap for conversion testing
140   - GWcsCAR car(coords, crval1, crval2, crpix1, crpix2,
141   - cdelt1, cdelt2, cd, pv2);
  137 + GWcsCAR car(coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2);
142 138  
143 139 // Test (x,y) to (lon,lat)
144 140 try {
... ...