GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GFits.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GFits.cpp - FITS file access class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2008-2021 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GFits.cpp
23  * @brief FITS file access class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <cstdio>
32 #include "GException.hpp"
33 #include "GTools.hpp"
34 #include "GFitsCfitsio.hpp"
35 #include "GFits.hpp"
36 #include "GFitsImageByte.hpp"
37 #include "GFitsImageSByte.hpp"
38 #include "GFitsImageUShort.hpp"
39 #include "GFitsImageShort.hpp"
40 #include "GFitsImageULong.hpp"
41 #include "GFitsImageLong.hpp"
42 #include "GFitsImageLongLong.hpp"
43 #include "GFitsImageFloat.hpp"
44 #include "GFitsImageDouble.hpp"
45 #include "GFitsAsciiTable.hpp"
46 #include "GFitsBinTable.hpp"
47 #include "GVOClient.hpp"
48 #include "GVOTable.hpp"
49 
50 /* __ Method name definitions ____________________________________________ */
51 #define G_AT1 "GFits::at(int&)"
52 #define G_AT2 "GFits::at(std::string&)"
53 #define G_IMAGE1 "GFits::image(int&)"
54 #define G_IMAGE2 "GFits::image(std::string&)"
55 #define G_TABLE1 "GFits::table(int&)"
56 #define G_TABLE2 "GFits::table(std::string&)"
57 #define G_SET1 "GFits::set(int&, GFitsHDU&)"
58 #define G_SET2 "GFits::set(std::string&, GFitsHDU&)"
59 #define G_INSERT1 "GFits::insert(const int& extno, GFitsHDU& hdu)"
60 #define G_INSERT2 "GFits::insert(std::string&, GFitsHDU& hdu)"
61 #define G_REMOVE1 "GFits::remove(int&)"
62 #define G_REMOVE2 "GFits::remove(std::string&)"
63 #define G_OPEN "GFits::open(GFilename&, bool&)"
64 #define G_SAVE "GFits::save(bool&)"
65 #define G_SAVETO "GFits::saveto(GFilename&, bool&)"
66 #define G_PUBLISH "GFits::publish(std::string&, std::string&)"
67 #define G_FREE_MEM "GFits::free_members()"
68 #define G_NEW_IMAGE "GFits::new_image()"
69 
70 /* __ Macros _____________________________________________________________ */
71 
72 /* __ Coding definitions _________________________________________________ */
73 #define G_DELETE_EMPTY_FITS_FILES //!< Do not write empty FITS files
74 
75 /* __ Debug definitions __________________________________________________ */
76 //#define G_DEBUG //!< Debug methods
77 
78 
79 /*==========================================================================
80  = =
81  = Constructors/destructors =
82  = =
83  ==========================================================================*/
84 
85 /***********************************************************************//**
86  * @brief Void constructor
87  *
88  * Constructs empty FITS file.
89  ***************************************************************************/
91 {
92  // Initialise class members
93  init_members();
94 
95  // Return
96  return;
97 }
98 
99 
100 /***********************************************************************//**
101  * @brief FITS file constructor
102  *
103  * @param[in] filename FITS file name.
104  * @param[in] create Create file if it does not exist? (default: false)
105  *
106  * Construct an object by opening a FITS file. If the file does not exist it
107  * will be created if @p create is set to true.
108  ***************************************************************************/
109 GFits::GFits(const GFilename& filename, const bool& create)
110 {
111  // Initialise class members
112  init_members();
113 
114  // Open specified FITS file
115  open(filename, create);
116 
117  // Return
118  return;
119 }
120 
121 
122 /***********************************************************************//**
123  * @brief Copy constructor
124  *
125  * @param[in] fits FITS file.
126  ***************************************************************************/
127 GFits::GFits(const GFits& fits)
128 {
129  // Initialise class members
130  init_members();
131 
132  // Copy members
133  copy_members(fits);
134 
135  // Return
136  return;
137 }
138 
139 
140 /***********************************************************************//**
141  * @brief Destructor
142  ***************************************************************************/
144 {
145  // Free members
146  free_members();
147 
148  // Return
149  return;
150 }
151 
152 
153 /*==========================================================================
154  = =
155  = Operators =
156  = =
157  ==========================================================================*/
158 
159 /***********************************************************************//**
160  * @brief Assignment operator
161  *
162  * @param[in] fits FITS file.
163  * @return FITS file.
164  ***************************************************************************/
166 {
167  // Execute only if object is not identical
168  if (this != &fits) {
169 
170  // Free members
171  free_members();
172 
173  // Initialise private members for clean destruction
174  init_members();
175 
176  // Copy members
177  copy_members(fits);
178 
179  } // endif: object was not identical
180 
181  // Return this object
182  return *this;
183 }
184 
185 
186 /*==========================================================================
187  = =
188  = Public methods =
189  = =
190  ==========================================================================*/
191 
192 /***********************************************************************//**
193  * @brief Clear FITS file
194  ***************************************************************************/
195 void GFits::clear(void)
196 {
197  // Close file and free all members
198  free_members();
199 
200  // Initialise members
201  init_members();
202 
203  // Return
204  return;
205 }
206 
207 
208 /***********************************************************************//**
209  * @brief Clone FITS file
210  *
211  * @return Pointer to deep copy of FITS file.
212  ***************************************************************************/
213 GFits* GFits::clone(void) const
214 {
215  return new GFits(*this);
216 }
217 
218 
219 /***********************************************************************//**
220  * @brief Get pointer to HDU
221  *
222  * @param[in] extno Extension number [0,...,size()-1].
223  * @return Pointer to HDU.
224  *
225  * @exception GException::out_of_range
226  * Extension number is out of range.
227  * @exception GException::runtime_error
228  * No HDU found for specified extension number.
229  *
230  * Returns a pointer to the HDU with the specified extension number @p extno.
231  * An exception is thrown if the HDU is not valid (i.e. NULL).
232  ***************************************************************************/
233 GFitsHDU* GFits::at(const int& extno)
234 {
235  // Compile option: raise an exception if index is out of range
236  #if defined(G_RANGE_CHECK)
237  if (extno < 0 || extno >= size()) {
238  throw GException::out_of_range(G_AT1, "FITS extension number",
239  extno, size());
240  }
241  #endif
242 
243  // Get HDU pointer
244  GFitsHDU* ptr = m_hdu[extno];
245 
246  // Throw an error if HDU has not been found
247  if (ptr == NULL) {
248  std::string msg = "FITS HDU with extension number "+
249  gammalib::str(extno)+" not found in FITS file.";
250  throw GException::runtime_error(G_AT1, msg);
251  }
252 
253  // Return pointer
254  return ptr;
255 }
256 
257 
258 /***********************************************************************//**
259  * @brief Get pointer to HDU (const version)
260  *
261  * @param[in] extno Extension number [0,...,size()-1].
262  * @return Pointer to HDU.
263  *
264  * @exception GException::out_of_range
265  * Extension number is out of range.
266  * @exception GException::runtime_error
267  * No HDU found for specified extension number.
268  *
269  * Returns a pointer to the HDU with the specified extension number @p extno.
270  * An exception is thrown if the HDU is not valid (i.e. NULL).
271  ***************************************************************************/
272 const GFitsHDU* GFits::at(const int& extno) const
273 {
274  // Compile option: raise an exception if index is out of range
275  #if defined(G_RANGE_CHECK)
276  if (extno < 0 || extno >= size()) {
277  throw GException::out_of_range(G_AT1, "FITS extension number",
278  extno, size());
279  }
280  #endif
281 
282  // Get HDU pointer
283  GFitsHDU* ptr = m_hdu[extno];
284 
285  // Throw an error if HDU has not been found
286  if (ptr == NULL) {
287  std::string msg = "FITS HDU with extension number "+
288  gammalib::str(extno)+" not found in FITS file.";
289  throw GException::runtime_error(G_AT1, msg);
290  }
291 
292  // Return pointer
293  return ptr;
294 }
295 
296 
297 /***********************************************************************//**
298  * @brief Get pointer to HDU
299  *
300  * @param[in] extname Name of HDU extension.
301  * @return Pointer to HDU.
302  *
303  * @exception GException::invalid_argument
304  * No HDU with specified name has been found.
305  *
306  * Returns a pointer to the HDU with the specified @p extname. An exception
307  * is thrown if the HDU is not valid (i.e. NULL).
308  ***************************************************************************/
309 GFitsHDU* GFits::at(const std::string& extname)
310 {
311  // Get extenion number
312  int extno = this->extno(extname);
313  if (extno == -1) {
314  std::string msg = "FITS extension \""+extname+"\" not found in FITS "
315  "file. Please specify a valid extension name.";
317  }
318 
319  // Get HDU pointer
320  GFitsHDU* ptr = m_hdu[extno];
321 
322  // Return pointer
323  return ptr;
324 }
325 
326 
327 /***********************************************************************//**
328  * @brief Get pointer to HDU (const version)
329  *
330  * @param[in] extname Name of HDU extension.
331  * @return Pointer to HDU.
332  *
333  * @exception GException::invalid_argument
334  * No HDU with specified name has been found.
335  *
336  * Returns a pointer to the HDU with the specified @p extname. An exception
337  * is thrown if the HDU is not valid (i.e. NULL).
338  ***************************************************************************/
339 const GFitsHDU* GFits::at(const std::string& extname) const
340 {
341  // Get extenion number
342  int extno = this->extno(extname);
343  if (extno == -1) {
344  std::string msg = "FITS extension \""+extname+"\" not found in FITS "
345  "file. Please specify a valid extension name.";
347  }
348 
349  // Get HDU pointer
350  GFitsHDU* ptr = m_hdu[extno];
351 
352  // Return pointer
353  return ptr;
354 }
355 
356 
357 /***********************************************************************//**
358  * @brief Get pointer to image HDU
359  *
360  * @param[in] extno Extension number [0,...,size()-1].
361  * @return Pointer to FITS image.
362  *
363  * @exception GException::invalid_argument
364  * Requested HDU is not an image.
365  *
366  * Returns a pointer to the image HDU with extension number extno.
367  ***************************************************************************/
368 GFitsImage* GFits::image(const int& extno)
369 {
370  // Get HDU pointer
371  GFitsImage* ptr = dynamic_cast<GFitsImage*>(at(extno));
372 
373  // Throw an error if HDU is not an image
374  if (ptr == NULL) {
375  std::string msg = "FITS extension number "+gammalib::str(extno)+" is "
376  "not an image. Please specify the extension number "
377  "of an image.";
379  }
380 
381  // Return pointer
382  return ptr;
383 }
384 
385 
386 /***********************************************************************//**
387  * @brief Get pointer to image HDU (const version)
388  *
389  * @param[in] extno Extension number [0,...,size()-1].
390  * @return Pointer to FITS image.
391  *
392  * @exception GException::invalid_argument
393  * Requested HDU is not an image.
394  *
395  * Returns a pointer to the image HDU with extension number extno.
396  ***************************************************************************/
397 const GFitsImage* GFits::image(const int& extno) const
398 {
399  // Get HDU pointer
400  const GFitsImage* ptr = dynamic_cast<const GFitsImage*>(at(extno));
401 
402  // Throw an error if HDU is not an image
403  if (ptr == NULL) {
404  std::string msg = "FITS extension number "+gammalib::str(extno)+" is "
405  "not an image. Please specify the extension number "
406  "of an image.";
408  }
409 
410  // Return pointer
411  return ptr;
412 }
413 
414 
415 /***********************************************************************//**
416  * @brief Get pointer to image HDU
417  *
418  * @param[in] extname Name of HDU extension.
419  * @return Pointer to FITS image.
420  *
421  * @exception GException::invalid_argument
422  * Requested HDU is not an image.
423  *
424  * Returns a pointer to the image HDU with extension name extname.
425  ***************************************************************************/
426 GFitsImage* GFits::image(const std::string& extname)
427 {
428  // Get HDU pointer
429  GFitsImage* ptr = dynamic_cast<GFitsImage*>(at(extname));
430 
431  // Throw an error if HDU is not an image
432  if (ptr == NULL) {
433  std::string msg = "FITS extension \""+extname+"\" is not an image. "
434  "Please specify the name of an image extension.";
436  }
437 
438  // Return pointer
439  return ptr;
440 }
441 
442 
443 /***********************************************************************//**
444  * @brief Get pointer to image HDU (const version)
445  *
446  * @param[in] extname Name of HDU extension.
447  * @return Pointer to FITS image.
448  *
449  * @exception GException::invalid_argument
450  * Requested HDU is not an image.
451  *
452  * Returns a pointer to the image HDU with extension name extname.
453  ***************************************************************************/
454 const GFitsImage* GFits::image(const std::string& extname) const
455 {
456  // Get HDU pointer
457  const GFitsImage* ptr = dynamic_cast<const GFitsImage*>(at(extname));
458 
459  // Throw an error if HDU is not an image
460  if (ptr == NULL) {
461  std::string msg = "FITS extension \""+extname+"\" is not an image. "
462  "Please specify the name of an image extension.";
464  }
465 
466  // Return pointer
467  return ptr;
468 }
469 
470 
471 /***********************************************************************//**
472  * @brief Get pointer to table HDU
473  *
474  * @param[in] extno Extension number [0,...,size()-1].
475  * @return Pointer to FITS table.
476  *
477  * @exception GException::invalid_argument
478  * Requested HDU is not a table.
479  *
480  * Returns a pointer to the table HDU with extension number extno.
481  ***************************************************************************/
482 GFitsTable* GFits::table(const int& extno)
483 {
484  // Get HDU pointer
485  GFitsTable* ptr = dynamic_cast<GFitsTable*>(at(extno));
486 
487  // Throw an error if HDU is not an image
488  if (ptr == NULL) {
489  std::string msg = "FITS extension number "+gammalib::str(extno)+" is "
490  "not a table. Please specify the extension number "
491  "of a table.";
493  }
494 
495  // Return pointer
496  return ptr;
497 }
498 
499 
500 /***********************************************************************//**
501  * @brief Get pointer to table HDU (const version)
502  *
503  * @param[in] extno Extension number [0,...,size()-1].
504  * @return Pointer to FITS table.
505  *
506  * @exception GException::invalid_argument
507  * Requested HDU is not a table.
508  *
509  * Returns a pointer to the table HDU with extension number extno.
510  ***************************************************************************/
511 const GFitsTable* GFits::table(const int& extno) const
512 {
513  // Get HDU pointer
514  const GFitsTable* ptr = dynamic_cast<const GFitsTable*>(at(extno));
515 
516  // Throw an error if HDU is not an image
517  if (ptr == NULL) {
518  std::string msg = "FITS extension number "+gammalib::str(extno)+" is "
519  "not a table. Please specify the extension number "
520  "of a table.";
522  }
523 
524  // Return pointer
525  return ptr;
526 }
527 
528 
529 /***********************************************************************//**
530  * @brief Get pointer to table HDU
531  *
532  * @param[in] extname Name of HDU extension.
533  * @return Pointer to FITS table.
534  *
535  * @exception GException::invalid_argument
536  * Requested HDU is not a table.
537  *
538  * Returns a pointer to the table HDU with extension name extname.
539  ***************************************************************************/
540 GFitsTable* GFits::table(const std::string& extname)
541 {
542  // Get HDU pointer
543  GFitsTable* ptr = dynamic_cast<GFitsTable*>(at(extname));
544 
545  // Throw an error if HDU is not an image
546  if (ptr == NULL) {
547  std::string msg = "FITS extension \""+extname+"\" is not a table. "
548  "Please specify the name of a table extension.";
550  }
551 
552  // Return pointer
553  return ptr;
554 }
555 
556 
557 /***********************************************************************//**
558  * @brief Get pointer to table HDU
559  *
560  * @param[in] extname Name of HDU extension.
561  * @return Pointer to FITS table.
562  *
563  * @exception GException::invalid_argument
564  * Requested HDU is not a table.
565  *
566  * Returns a pointer to the table HDU with extension name extname.
567  ***************************************************************************/
568 const GFitsTable* GFits::table(const std::string& extname) const
569 {
570  // Get HDU pointer
571  const GFitsTable* ptr = dynamic_cast<const GFitsTable*>(at(extname));
572 
573  // Throw an error if HDU is not an image
574  if (ptr == NULL) {
575  std::string msg = "FITS extension \""+extname+"\" is not a table. "
576  "Please specify the name of a table extension.";
578  }
579 
580  // Return pointer
581  return ptr;
582 }
583 
584 
585 /***********************************************************************//**
586  * @brief Set HDU for the specified extension number
587  *
588  * @param[in] extno Extension number [0,...,size()-1].
589  * @param[in] hdu HDU.
590  * @return Pointer to cloned HDU.
591  *
592  * @exception GException::out_of_range
593  * Extension number out of range.
594  * @exception GException::invalid_argument
595  * Attempt to insert non-image HDU in first slot.
596  ***************************************************************************/
597 GFitsHDU* GFits::set(const int& extno, const GFitsHDU& hdu)
598 {
599  // Compile option: raise exception if index is out of range
600  #if defined(G_RANGE_CHECK)
601  if (extno < 0 || extno >= size()) {
602  throw GException::out_of_range(G_SET1, "FITS extension number",
603  extno, size());
604  }
605  #endif
606 
607  // Throw an exception if a non-image extension should be set in the
608  // first slot
609  if (extno == 0 && hdu.exttype() != GFitsHDU::HT_IMAGE) {
610  std::string msg = "Attempt to set a table extension as the first"
611  " extension of a FITS file.\n"
612  "The first extension of a FITS file must be an"
613  " image, hence use the next slot to set a"
614  " table.";
616  }
617 
618  // Free existing HDU only if it differs from current HDU. This prevents
619  // unintential deallocation of the argument
620  if ((m_hdu[extno] != NULL) && (m_hdu[extno] != &hdu)) {
621  delete m_hdu[extno];
622  }
623 
624  // Assign new HDU by cloning
625  m_hdu[extno] = hdu.clone();
626 
627  // If FITS file exists then connect cloned HDU to FITS file
628  if (m_fitsfile != NULL) {
629  __fitsfile fptr = *(FPTR(m_fitsfile));
630  fptr.HDUposition = extno;
631  m_hdu[extno]->connect(&fptr);
632  }
633  else {
634  m_hdu[extno]->extno(extno);
635  }
636 
637  // Return pointer to HDU
638  return (m_hdu[extno]);
639 }
640 
641 
642 /***********************************************************************//**
643  * @brief Set HDU for the specified extension name
644  *
645  * @param[in] extname Name of HDU extension.
646  * @param[in] hdu HDU.
647  * @return Pointer to cloned HDU.
648  *
649  * @exception GException::invalid_argument
650  * Extension name not found.
651  ***************************************************************************/
652 GFitsHDU* GFits::set(const std::string& extname, const GFitsHDU& hdu)
653 {
654  // Get extenion number
655  int extno = this->extno(extname);
656  if (extno == -1) {
657  std::string msg = "No extension with name \""+extname+"\" found in "
658  "FITS file. Please specify a valid extension name.";
660  }
661 
662  // Set HDU and return pointer
663  return (set(extno, hdu));
664 }
665 
666 
667 /***********************************************************************//**
668  * @brief Append HDU to FITS file
669  *
670  * @param[in] hdu HDU.
671  * @return Pointer to appended HDU.
672  *
673  * Append HDU to the next free position in a FITS file. In case that no HDU
674  * exists so far in the FITS file and if the HDU to append is NOT an image,
675  * an empty primary image will be inserted as first HDU in the FITS file.
676  * This guarantees the compatibility with the FITS standard.
677  ***************************************************************************/
679 {
680  // Debug header
681  #if defined(G_DEBUG)
682  std::cout << "GFits::append(";
683  switch (hdu.exttype()) {
684  case GFitsHDU::HT_IMAGE:
685  std::cout << "GFitsImage";
686  break;
688  std::cout << "GFitsAsciiTable";
689  break;
691  std::cout << "GFitsBinTable";
692  break;
693  default:
694  std::cout << "<unknown header type>";
695  break;
696  }
697  std::cout << ") (size=" << size() << "): entry" << std::endl;
698  #endif
699 
700  // Determine next free HDU number
701  int n_hdu = size();
702 
703  // Add primary image if required
704  if (n_hdu == 0 && hdu.exttype() != GFitsHDU::HT_IMAGE) {
705 
706  // Allocate primary image
707  GFitsHDU* primary = new_primary();
708 
709  // If FITS file exists then connect primary image to FITS file
710  if (m_fitsfile != NULL) {
711  __fitsfile fptr = *(FPTR(m_fitsfile));
712  fptr.HDUposition = n_hdu;
713  primary->connect(&fptr);
714  }
715 
716  // Debug information
717  #if defined(G_DEBUG)
718  std::cout << "GFits::append: append primary image to HDU.";
719  std::cout << std::endl;
720  #endif
721 
722  // Push back primary image
723  m_hdu.push_back(primary);
724 
725  // Increment HDU number
726  n_hdu++;
727  }
728 
729  // Clone FITS HDU
730  GFitsHDU* ptr = hdu.clone();
731 
732  // Append HDU if it's valid
733  if (ptr != NULL) {
734 
735  // If FITS file exists then connect cloned HDU to FITS file
736  if (m_fitsfile != NULL) {
737  __fitsfile fptr = *(FPTR(m_fitsfile));
738  fptr.HDUposition = n_hdu;
739  ptr->connect(&fptr);
740  }
741  else {
742  ptr->extno(n_hdu);
743  }
744 
745  // Debug information
746  #if defined(G_DEBUG)
747  std::cout << "GFits::append: append HDU (extno=" << n_hdu << ")";
748  std::cout << std::endl;
749  #endif
750 
751  // Push back HDU
752  m_hdu.push_back(ptr);
753 
754  // Debug trailer
755  #if defined(G_DEBUG)
756  std::cout << "GFits::append: exit" << std::endl;
757  #endif
758 
759  } // endif: HDU was valid
760 
761  // Return
762  return ptr;
763 }
764 
765 
766 /***********************************************************************//**
767  * @brief Set HDU for the specified extension number
768  *
769  * @param[in] extno Extension number [0,...,size()-1].
770  * @param[in] hdu HDU.
771  * @return Pointer to cloned HDU.
772  *
773  * @exception GException::out_of_range
774  * Extension number out of range.
775  * @exception GException::invalid_argument
776  * Attempt to insert non-image HDU in first slot.
777  ***************************************************************************/
778 GFitsHDU* GFits::insert(const int& extno, const GFitsHDU& hdu)
779 {
780  // Compile option: raise exception if extension number is out of range
781  #if defined(G_RANGE_CHECK)
782  if (extno < 0 || extno >= size()) {
783  throw GException::out_of_range(G_INSERT1, "FITS extension number",
784  extno, size());
785  }
786  #endif
787 
788  // Throw an exception if a non-image extension should be inserted in the
789  // first slot
790  if (extno == 0 && hdu.exttype() != GFitsHDU::HT_IMAGE) {
791  std::string msg = "Attempt to insert a table extension as the first "
792  "extension of a FITS file. The first extension of "
793  "a FITS file must be an image, hence use the next "
794  "slot to insert a table.";
796  }
797 
798  // Create deep copy of HDU
799  GFitsHDU* ptr = hdu.clone();
800 
801  // Inserts deep copy of HDU
802  m_hdu.insert(m_hdu.begin() + extno, ptr);
803 
804  // If FITS file exists then connect cloned HDU to FITS file
805  if (m_fitsfile != NULL) {
806  __fitsfile fptr = *(FPTR(m_fitsfile));
807  fptr.HDUposition = extno;
808  ptr->connect(&fptr);
809  }
810  else {
811  ptr->extno(extno);
812  }
813 
814  // Update extno for all subsequent HDUs
815  for (int i = extno + 1; i < size(); ++i) {
816  m_hdu[i]->extno(i);
817  }
818 
819  // Return pointer to HDU
820  return ptr;
821 }
822 
823 
824 /***********************************************************************//**
825  * @brief Insert HDU into FITS file
826  *
827  * @param[in] extname Extension name.
828  * @param[in] hdu Extension.
829  * @return Pointer to deep copy of extension @p hdu.
830  *
831  * @exception GException::invalid_argument
832  * Extension name not found.
833  *
834  * Inserts an extension @p hdu into the FITS file before the extension with
835  * the specified @p extname.
836  ***************************************************************************/
837 GFitsHDU* GFits::insert(const std::string& extname, const GFitsHDU& hdu)
838 {
839  // Get extenion number
840  int extno = this->extno(extname);
841  if (extno == -1) {
842  std::string msg = "No extension with name \""+extname+"\" found in "
843  "FITS file. Please specify a valid extension name.";
845  }
846 
847  // Insert HDU and return pointer
848  return (insert(extno, hdu));
849 }
850 
851 
852 /***********************************************************************//**
853  * @brief Remove HDU from FITS file
854  *
855  * @param[in] extno Extension number [0,...,size()-1].
856  *
857  * @exception GException::out_of_range
858  * Specified @p extno is out of range.
859  *
860  * @todo Handle HDU update in FITS file.
861  ***************************************************************************/
862 void GFits::remove(const int& extno)
863 {
864  // Compile option: raise an exception if index is out of range
865  #if defined(G_RANGE_CHECK)
866  if (extno < 0 || extno >= size()) {
867  throw GException::out_of_range(G_REMOVE1, "FITS extension number",
868  extno, size());
869  }
870  #endif
871 
872  // Throw an exception if removal would lead to a non-image extension
873  // in the first slot
874  if (extno == 0 && size() > 1 && m_hdu[1]->exttype() != GFitsHDU::HT_IMAGE) {
875  std::string msg = "Attempt to remove primary image from a FITS file "
876  "with a table extension as second extension. The "
877  "removal of the primary image would result in "
878  "having a table as the first extension of the FITS "
879  "file, which is not a valid FITS file.";
881  }
882 
883  // Delete HDU
884  if (m_hdu[extno] != NULL) delete m_hdu[extno];
885 
886  // Erase HDU from FITS file
887  m_hdu.erase(m_hdu.begin() + extno);
888 
889  // Update extno for all subsequent HDUs
890  for (int i = extno; i < size(); ++i) {
891  m_hdu[i]->extno(i);
892  }
893 
894  // Return
895  return;
896 }
897 
898 
899 /***********************************************************************//**
900  * @brief Remove HDU from FITS file
901  *
902  * @param[in] extname Name of HDU extension.
903  *
904  * @exception GException::invalid_argument
905  * Specified @p extname not found in FITS file.
906  ***************************************************************************/
907 void GFits::remove(const std::string& extname)
908 {
909  // Get extenion number
910  int extno = this->extno(extname);
911  if (extno == -1) {
912  std::string msg = "No extension with name \""+extname+"\" found in "
913  "FITS file. Please specify a valid extension name.";
915  }
916 
917  // Remove extension
918  remove(extno);
919 
920  // Return
921  return;
922 }
923 
924 
925 /***********************************************************************//**
926  * @brief Append FITS file
927  *
928  * @param[in] fits FITS file.
929  *
930  * Append all extension of @p fits file to FITS file.
931  *
932  * @todo Handle HDU update in FITS file.
933  ***************************************************************************/
934 void GFits::extend(const GFits& fits)
935 {
936  // Do nothing if FITS file is empty
937  if (!fits.is_empty()) {
938 
939  // Get size. Note that we extract the size first to avoid an
940  // endless loop that arises when a container is appended to
941  // itself.
942  int num = fits.size();
943 
944  // Reserve enough space
945  reserve(size() + num);
946 
947  // Loop over all HDUs and append pointers to deep copies
948  for (int i = 0; i < num; ++i) {
949 
950  // Clone HDU
951  GFitsHDU* ptr = m_hdu[i]->clone();
952 
953  // Push back HDU on stack
954  m_hdu.push_back(ptr);
955 
956  // Retrieve extno
957  int extno = m_hdu.size()-1;
958 
959  // If FITS file exists then connect cloned HDU to FITS file
960  if (m_fitsfile != NULL) {
961  __fitsfile fptr = *(FPTR(m_fitsfile));
962  fptr.HDUposition = extno;
963  ptr->connect(&fptr);
964  }
965  else {
966  ptr->extno(extno);
967  }
968 
969  } // endfor: looped over all HDUs
970 
971  } // endif: FITS file was not empty
972 
973  // Return
974  return;
975 }
976 
977 
978 /***********************************************************************//**
979  * @brief Return extension number for specified extension name
980  *
981  * @param[in] extname Name of HDU extension.
982  * @return Extension number (-1 of not found).
983  *
984  * Returns the extension number for a specified extension name @p extname. If
985  * the extension name if "PRIMARY" an extension number of 0 is returned.
986  * If the extension is not found, -1 is returned.
987  *
988  * @p extname is case in-sensitive.
989  ***************************************************************************/
990 int GFits::extno(const std::string& extname) const
991 {
992  // Initialise extension number
993  int extno = -1;
994 
995  // Get extname in upper case
996  std::string uextname = gammalib::toupper(extname);
997 
998  // Return primary HDU if requested ...
999  if (uextname == "PRIMARY") {
1000  if (size() > 0) {
1001  extno = 0;
1002  }
1003  }
1004 
1005  // ... otherwise search for specified extension
1006  else {
1007  for (int i = 0; i < size(); ++i) {
1008  if (gammalib::toupper(m_hdu[i]->extname()) == uextname) {
1009  extno = i;
1010  break;
1011  }
1012  }
1013  }
1014 
1015  // Return extno
1016  return extno;
1017 }
1018 
1019 
1020 /***********************************************************************//**
1021  * @brief Open or (optionally) create FITS file
1022  *
1023  * @param[in] filename Name of FITS file to be opened.
1024  * @param[in] create Create FITS file if it does not exist (default: false)
1025  *
1026  * @exception GException::invalid_argument
1027  * Class instance contains already an opened FITS file.
1028  * Close file before opening a new one using GFits::close().
1029  * @exception GException::fits_error
1030  * Unable to determine number of HDUs in the FITS file.
1031  * @exception GException::runtime_error
1032  * Unknown HDU type encountered.
1033  *
1034  * This method opens all HDUs that are found in the specified FITS file.
1035  * If the file does not exist, and if create=true, a new FITS file is created.
1036  * For each HDU, a GFitsHDU object is associated to the GFits object.
1037  * The HDUs can then be accessed using the hdu(const std::string&) or
1038  * hdu(int extno) method.
1039  * Any environment variable present in the filename will be expanded.
1040  ***************************************************************************/
1041 void GFits::open(const GFilename& filename, const bool& create)
1042 {
1043  // Don't allow opening if a file is already open
1044  if (m_fitsfile != NULL) {
1045  std::string msg;
1046  if (m_filename == filename) {
1047  msg = "FITS file \""+m_filename+"\" has already been "
1048  "opened, cannot open it again.";
1049  }
1050  else {
1051  msg = "A FITS file \""+m_filename+"\" has already "
1052  "been opened, cannot open another FITS file \""+
1053  filename+"\" before closing the existing one.";
1054  }
1056  }
1057 
1058  // Remove any HDUs
1059  m_hdu.clear();
1060 
1061  // Initialise FITS file as readwrite and non created
1062  m_readwrite = true;
1063  m_created = false;
1064 
1065  // Try opening FITS file with readwrite access
1066  int status = 0;
1067  status = __ffopen(FHANDLE(m_fitsfile), std::string(filename).c_str(), 1, &status);
1068 
1069  // If failed then try opening as readonly.
1070  //
1071  // Possible error codes are:
1072  // 104: FILE_NOT_OPENED (could not open the named file)
1073  // 107: END_OF_FILE (tried to move past end of file, occurs if extname is invalid)
1074  // 112: READONLY_FILE (cannot write to readonly file)
1075  if (status == 104 || status == 107 || status == 112) {
1076  status = 0;
1077  status = __ffopen(FHANDLE(m_fitsfile), std::string(filename).c_str(), 0, &status);
1078  m_readwrite = false;
1079  }
1080 
1081  // If failed and if we are allowed to create a new FITS file then create
1082  // FITS file now. We pass the URL as filename since we want to create a
1083  // fresh file without any extension. We also remove any file before in
1084  // case that file exists but is not a FITS file.
1085  //
1086  // Possible error codes are:
1087  // 104: FILE_NOT_OPENED (could not open the named file)
1088  // 107: END_OF_FILE (tried to move past end of file, occurs if extname is invalid)
1089  if (create && (status == 104 || status == 107)) {
1090  filename.remove();
1091  status = 0;
1092  status = __ffinit(FHANDLE(m_fitsfile), filename.url().c_str(), &status);
1093  m_readwrite = true;
1094  m_created = true;
1095  }
1096 
1097  // Throw special exception if status=202 (keyword not found). This error
1098  // may occur if the file is opened with an expression
1099  if (status == 202) {
1100  std::string msg = "Keyword not found when opening FITS file \""+
1101  filename+"\". Please specify a valid keyword.";
1102  throw GException::fits_error(G_OPEN, status, msg);
1103  }
1104 
1105  // Throw any other error
1106  else if (status != 0) {
1107  std::string msg = "Error occured when opening FITS file \""+
1108  filename+"\".";
1109  throw GException::fits_error(G_OPEN, status, msg);
1110 
1111  }
1112 
1113  // Store FITS file name as GFilename object
1114  m_filename = filename;
1115 
1116  // Determine number of HDUs
1117  int num_hdu = 0;
1118  status = __ffthdu(FPTR(m_fitsfile), &num_hdu, &status);
1119  if (status != 0) {
1120  throw GException::fits_error(G_OPEN, status);
1121  }
1122 
1123  // Open and append all HDUs
1124  for (int i = 0; i < num_hdu; ++i) {
1125 
1126  // Move to HDU and get HDU type
1127  int type = gammalib::fits_move_to_hdu(G_OPEN, m_fitsfile, i+1);
1128 
1129  // Perform type dependent HDU allocation
1130  GFitsHDU* hdu = NULL;
1131  switch (type) {
1132  case GFitsHDU::HT_IMAGE:
1133  hdu = new_image();
1134  break;
1136  hdu = new GFitsAsciiTable;
1137  break;
1139  hdu = new GFitsBinTable;
1140  break;
1141  default:
1142  std::string msg = "Unknown HDU type \""+gammalib::str(type)+"\" "
1143  "encountered.";
1144  throw GException::runtime_error(G_OPEN, msg);
1145  break;
1146  }
1147 
1148  // Open HDU
1149  hdu->open(FPTR(m_fitsfile), i);
1150 
1151  // Append HDU
1152  m_hdu.push_back(hdu);
1153 
1154  } // endfor: looped over all HDUs
1155 
1156  // Return
1157  return;
1158 }
1159 
1160 
1161 /***********************************************************************//**
1162  * @brief Saves FITS file
1163  *
1164  * @param[in] clobber Overwrite existing FITS file? (default: false).
1165  *
1166  * @exception GException::invalid_value
1167  * Attemting to overwrite an existing file without having specified
1168  * clobber=true.
1169  * FITS file needs to be opened before saving.
1170  *
1171  * Saves all HDUs of an open FITS file to disk. After saving, the FITS file
1172  * remains open. Invoke the close() method if explicit closing is needed.
1173  * Note that de-allocation of the GFits object also closes the FITS file.
1174  *
1175  * In the special case that no first HDU exists an empty primary image is
1176  * created.
1177  ***************************************************************************/
1178 void GFits::save(const bool& clobber)
1179 {
1180  // Debug header
1181  #if defined(G_DEBUG)
1182  std::cout << "GFits::save (size=" << size() << "): entry" << std::endl;
1183  #endif
1184 
1185  // Initialise cfitsio status
1186  int status = 0;
1187 
1188  // If we attempt to save an existing file without overwriting permission
1189  // then throw an error
1190  if (!m_created && !clobber) {
1191  std::string msg = "Attempted to overwrite FITS file \""+m_filename+
1192  "\". Please set clobber flag to true.";
1193  throw GException::invalid_value(G_SAVE, msg);
1194  }
1195 
1196  // If no FITS file has been opened then throw an error
1197  if (m_fitsfile == NULL) {
1198  std::string msg = "Cannot save FITS file that was not open. Please "
1199  "open FITS file before saving or use saveto() "
1200  "method.";
1201  throw GException::invalid_value(G_SAVE, msg);
1202  }
1203 
1204  // We have to make sure that all data are online before we can save them.
1205  // This needs to be done before the data are deleted from disk. We do
1206  // this by cloning all HDUs, as cloning will force all data to get online.
1207  for (int i = 0; i < size(); ++i) {
1208  GFitsHDU* hdu = m_hdu[i]->clone();
1209  delete m_hdu[i];
1210  m_hdu[i] = hdu;
1211  }
1212 
1213  // Determine number of HDUs in current FITS file on disk
1214  int num_hdu = 0;
1215  status = __ffthdu(FPTR(m_fitsfile), &num_hdu, &status);
1216  if (status != 0) {
1217  throw GException::fits_error(G_SAVE, status);
1218  }
1219 
1220  // Debug
1221  #if defined(G_DEBUG)
1222  std::cout << "GFits::save: currently " << num_hdu << " in file";
1223  std::cout << std::endl;
1224  #endif
1225 
1226  // Delete all HDUs (except of the primary HDU) since we will write them
1227  // all freshly
1228  for (int i = num_hdu-1; i >= 0; --i) {
1229  status = __ffdhdu(FPTR(m_fitsfile), NULL, &status);
1230  if (status != 0) {
1231  throw GException::fits_error(G_SAVE, status);
1232  }
1233  }
1234 
1235  // Debug
1236  #if defined(G_DEBUG)
1237  std::cout << "GFits::save: deleted HDUs" << std::endl;
1238  #endif
1239 
1240  // If no HDUs exist in the FITS object then write an empty primary image
1241  if (size() == 0) {
1242  status = __ffmahd(FPTR(m_fitsfile), 1, NULL, &status);
1243  status = __ffcrim(FPTR(m_fitsfile), 8, 0, NULL, &status);
1244  if (status != 0) {
1245  throw GException::fits_error(G_SAVE, status);
1246  }
1247  }
1248 
1249  // ... otherwise save all HDUs
1250  else {
1251  for (int i = 0; i < size(); ++i) {
1252  m_hdu[i]->extno(i);
1253  m_hdu[i]->save();
1254  }
1255  }
1256 
1257  // Debug
1258  #if defined(G_DEBUG)
1259  std::cout << "GFits::save: saved HDUs" << std::endl;
1260  #endif
1261 
1262  // Flush file to disk
1263  status = __ffflus(FPTR(m_fitsfile), &status);
1264  if (status != 0) {
1265  throw GException::fits_error(G_SAVE, status);
1266  }
1267 
1268  // Signal that file needs not to be created anymore
1269  m_created = false;
1270 
1271  // Debug trailer
1272  #if defined(G_DEBUG)
1273  std::cout << "GFits::save: exit" << std::endl;
1274  #endif
1275 
1276  // Return
1277  return;
1278 }
1279 
1280 
1281 /***********************************************************************//**
1282  * @brief Saves to specified FITS file
1283  *
1284  * @param[in] filename Filename.
1285  * @param[in] clobber Overwrite existing FITS file?
1286  *
1287  * @exception GException::invalid_value
1288  * Specified file exists already. Overwriting requires
1289  * clobber=true.
1290  *
1291  * Saves the FITS object into a specific file.
1292  ***************************************************************************/
1293 void GFits::saveto(const GFilename& filename, const bool& clobber)
1294 {
1295  // Debug header
1296  #if defined(G_DEBUG)
1297  std::cout << "GFits::saveto(\"" << filename << "\", " << clobber << ")"
1298  << " (size=" << size() << "): entry" << std::endl;
1299  #endif
1300 
1301  // If overwriting has been specified then remove the file. Otherwise,
1302  // throw an exception if the file exists.
1303  if (clobber) {
1304  filename.remove();
1305  }
1306  else if (filename.exists()) {
1307  std::string msg = "Attempted to overwrite FITS file \""+filename.url()+
1308  "\". Please set clobber flag to true.";
1309  throw GException::invalid_value(G_SAVETO, msg);
1310  }
1311 
1312  // Create or open FITS file
1313  GFits fits(filename, true);
1314 
1315  // Append all headers
1316  for (int i = 0; i < size(); ++i) {
1317  fits.append(*m_hdu[i]);
1318  }
1319 
1320  // Save FITS file
1321  fits.save();
1322 
1323  // Close FITS file
1324  fits.close();
1325 
1326  // Debug trailer
1327  #if defined(G_DEBUG)
1328  std::cout << "GFits::saveto: exit" << std::endl;
1329  #endif
1330 
1331  // Return
1332  return;
1333 }
1334 
1335 
1336 /***********************************************************************//**
1337  * @brief Close FITS file
1338  *
1339  * Closing detaches a FITS file from the GFits object and returns a clean
1340  * empty object.
1341  ***************************************************************************/
1342 void GFits::close(void)
1343 {
1344  // Close file and free all members
1345  free_members();
1346 
1347  // Initialise members
1348  init_members();
1349 
1350  // Return
1351  return;
1352 }
1353 
1354 
1355 /***********************************************************************//**
1356  * @brief Publish FITS HDU on VO Hub
1357  *
1358  * @param[in] extno Extension number.
1359  * @param[in] name Optional name of published ressource.
1360  *
1361  * Publishes FITS HDU with specified @p extno.
1362  ***************************************************************************/
1363 void GFits::publish(const int& extno, const std::string& name) const
1364 {
1365  // If HDU is an image then publish it as a FITS image
1366  if (this->at(extno)->exttype() == GFitsHDU::HT_IMAGE) {
1367 
1368  // Get HDU as an image
1369  const GFitsImage* image = this->image(extno);
1370 
1371  // Save extension name
1372  std::string extname = image->extname();
1373 
1374  // Optionally set extension name
1375  if (!name.empty()) {
1376  const_cast<GFitsImage*>(image)->extname(name);
1377  }
1378 
1379  // Create VO Client
1380  GVOClient client;
1381 
1382  // Publish image using VO client
1383  client.publish(*image);
1384 
1385  // Restore extension name
1386  const_cast<GFitsImage*>(image)->extname(extname);
1387 
1388  } // endif: HDU was an image
1389 
1390  // ... otherwise publish the HDU as a VO table
1391  else {
1392 
1393  // Get HDU as a table
1394  const GFitsTable* table = this->table(extno);
1395 
1396  // Save extension name
1397  std::string extname = table->extname();
1398 
1399  // Optionally set extension name
1400  if (!name.empty()) {
1401  const_cast<GFitsTable*>(table)->extname(name);
1402  }
1403 
1404  // Create VO Client
1405  GVOClient client;
1406 
1407  // Generate VO table FITS table
1408  GVOTable votable(*table);
1409 
1410  // Publish VO table using VO client
1411  client.publish(votable);
1412 
1413  // Restore extension name
1414  const_cast<GFitsTable*>(table)->extname(extname);
1415 
1416  } // endelse: published HDU as VO table
1417 
1418  // Return
1419  return;
1420 }
1421 
1422 
1423 /***********************************************************************//**
1424  * @brief Publish FITS HDU on VO Hub
1425  *
1426  * @param[in] extname Extension name.
1427  * @param[in] name Optional name of published ressource.
1428  *
1429  * Publishes FITS HDU with specified @p extname.
1430  ***************************************************************************/
1431 void GFits::publish(const std::string& extname, const std::string& name) const
1432 {
1433  // Get extenion number
1434  int extno = this->extno(extname);
1435  if (extno == -1) {
1436  std::string msg = "Extension \""+extname+"\" not found in FITS file. "
1437  "Please specify a valid extension name.";
1439  }
1440 
1441  // Publish FITS HDU
1442  publish(extno, name);
1443 
1444  // Return
1445  return;
1446 }
1447 
1448 
1449 /***********************************************************************//**
1450  * @brief Print FITS information
1451  *
1452  * @param[in] chatter Chattiness.
1453  * @return String containing FITS information.
1454  ***************************************************************************/
1455 std::string GFits::print(const GChatter& chatter) const
1456 {
1457  // Initialise result string
1458  std::string result;
1459 
1460  // Continue only if chatter is not silent
1461  if (chatter != SILENT) {
1462 
1463  // Append header
1464  result.append("=== GFits ===");
1465 
1466  // Append file information
1467  result.append("\n"+gammalib::parformat("Filename"));
1468  result.append(m_filename);
1469  result.append("\n"+gammalib::parformat("History"));
1470  if (m_created) {
1471  result.append("new file");
1472  }
1473  else {
1474  result.append("existing file");
1475  }
1476  result.append("\n"+gammalib::parformat("Mode"));
1477  if (m_readwrite) {
1478  result.append("read/write");
1479  }
1480  else {
1481  result.append("read only");
1482  }
1483  result.append("\n"+gammalib::parformat("Number of HDUs"));
1484  result.append(gammalib::str(size()));
1485 
1486  // NORMAL: Append HDUs
1487  if (chatter >= NORMAL) {
1488  for (int i = 0; i < size(); ++i) {
1489  result.append("\n");
1490  result.append(m_hdu[i]->print(gammalib::reduce(chatter)));
1491  }
1492  }
1493 
1494  } // endif: chatter was not silent
1495 
1496  // Return result
1497  return result;
1498 }
1499 
1500 
1501 /*==========================================================================
1502  = =
1503  = Private methods =
1504  = =
1505  ==========================================================================*/
1506 
1507 /***********************************************************************//**
1508  * @brief Initialise class members
1509  ***************************************************************************/
1511 {
1512  // Initialise GFits members
1513  m_hdu.clear();
1514  m_filename.clear();
1515  m_fitsfile = NULL;
1516  m_readwrite = true;
1517  m_created = true;
1518 
1519  // Return
1520  return;
1521 }
1522 
1523 
1524 /***********************************************************************//**
1525  * @brief Copy class members
1526  *
1527  * @param fits Object to be copied
1528  *
1529  * The method does not copy the FITS file name and FITS file pointer.
1530  * This prevents that several copies of the FITS file pointer exist in
1531  * different instances of GFits, which would lead to confusion since one
1532  * instance could close the file while for another it still would be
1533  * opened. The rule ONE INSTANCE - ONE FILE applies.
1534  ***************************************************************************/
1535 void GFits::copy_members(const GFits& fits)
1536 {
1537  // Reset FITS file attributes
1538  m_filename.clear();
1539  m_fitsfile = NULL;
1540  m_readwrite = true;
1541  m_created = true;
1542 
1543  // Clone HDUs
1544  m_hdu.clear();
1545  for (int i = 0; i < fits.m_hdu.size(); ++i) {
1546  m_hdu.push_back((fits.m_hdu[i]->clone()));
1547  }
1548 
1549  // Return
1550  return;
1551 }
1552 
1553 
1554 /***********************************************************************//**
1555  * @brief Delete class members
1556  *
1557  * Closes the FITS file if it had been opened and deallocate all HDUs.
1558  *
1559  * If the G_DELETE_EMPTY_FITS_FILES option is defined, files without HDUs
1560  * or corrupted files will be deleted. This prevents leaving corrupted files
1561  * on disk (yet, corrupted files may be generated by another application,
1562  * thus this is not 100% safe; better make the code solid against reading
1563  * corrupted FITS files).
1564  ***************************************************************************/
1566 {
1567  // If FITS file has been opened then close it now
1568  if (m_fitsfile != NULL) {
1569 
1570  // If file has been created but not yet saved then delete the file
1571  // now. We do not worry about the status in this case.
1572  if (m_created) {
1573  int status = 0;
1574  __ffdelt(FPTR(m_fitsfile), &status);
1575  }
1576 
1577  // Compile option: If there are no HDUs then delete the file (don't
1578  // worry about error)
1579  #if defined(G_DELETE_EMPTY_FITS_FILES)
1580  else if (size() == 0) {
1581  int status = 0;
1582  __ffdelt(FPTR(m_fitsfile), &status);
1583  }
1584  #endif
1585 
1586  // ... otherwise close the file
1587  else {
1588  int status = 0;
1589  status = __ffclos(FPTR(m_fitsfile), &status);
1590  if (status == 252) {
1591  int new_status = 0;
1592  __ffdelt(FPTR(m_fitsfile), &new_status);
1593  if (new_status != 0) {
1594  throw GException::fits_error(G_FREE_MEM, new_status);
1595  }
1596  }
1597  else if (status != 0) {
1598  throw GException::fits_error(G_FREE_MEM, status);
1599  }
1600 
1601  } // endelse: there was an open FITS file
1602 
1603  // Set FITS file pointer to NULL
1604  m_fitsfile = NULL;
1605 
1606  } // endif: FITS file had been opened
1607 
1608  // Free HDUs
1609  for (int i = 0; i < m_hdu.size(); ++i) {
1610  if (m_hdu[i] != NULL) {
1611  delete m_hdu[i];
1612  m_hdu[i] = NULL;
1613  }
1614  }
1615 
1616  // Clear HDUs
1617  m_hdu.clear();
1618 
1619  // Return
1620  return;
1621 }
1622 
1623 
1624 /***********************************************************************//**
1625  * @brief Allocate new FITS image and return memory pointer
1626  *
1627  * @exception runtime_error
1628  * Invalid number of bits per pixel.
1629  *
1630  * Depending on the number of bits per pixel, a FITS image is allocated
1631  * and the pointer is returned. The following FITS image classes are
1632  * handled:
1633  * GFitsImageByte (bitpix=8)
1634  * GFitsImageShort (bitpix=16)
1635  * GFitsImageLong (bitpix=32)
1636  * GFitsImageLongLong (bitpix=64)
1637  * GFitsImageFloat (bitpix=-32)
1638  * GFitsImageDouble (bitpix=-64)
1639  * The information about the number of bits per pixels is extracted from
1640  * the actual HDU.
1641  ***************************************************************************/
1643 {
1644  // Initialise return value
1645  GFitsImage* image = NULL;
1646 
1647  // Get number of bits per pixel
1648  int status = 0;
1649  int bitpix = -64;
1650  status = __ffgiet(FPTR(m_fitsfile), &bitpix, &status);
1651  if (status != 0) {
1652  throw GException::fits_error(G_NEW_IMAGE, status);
1653  }
1654 
1655  // Check for unsigned image
1656  char keyname[10];
1657  long long bzero;
1658  long long bscale;
1659  std::sprintf(keyname, "BZERO");
1660  if (__ffgky(FPTR(m_fitsfile), __TLONGLONG, keyname, &bzero, NULL, &status) != 0) {
1661  bzero = 0;
1662  }
1663  std::sprintf(keyname, "BSCALE");
1664  if (__ffgky(FPTR(m_fitsfile), __TLONGLONG, keyname, &bscale, NULL, &status) != 0) {
1665  bscale = 0;
1666  }
1667  if (bitpix == 8 && bzero == -128 && bscale == 1) {
1668  bitpix = 10;
1669  }
1670  else if (bitpix == 16 && bzero == 32768u && bscale == 1) {
1671  bitpix = 20;
1672  }
1673  else if (bitpix == 32 && bzero == 2147483648u && bscale == 1) {
1674  bitpix = 40;
1675  }
1676 
1677  // Allocate bitpix dependent image
1678  switch (bitpix) {
1679  case 8:
1680  image = new GFitsImageByte;
1681  break;
1682  case 10:
1683  image = new GFitsImageSByte;
1684  break;
1685  case 16:
1686  image = new GFitsImageShort;
1687  break;
1688  case 20:
1689  image = new GFitsImageUShort;
1690  break;
1691  case 32:
1692  image = new GFitsImageLong;
1693  break;
1694  case 40:
1695  image = new GFitsImageULong;
1696  break;
1697  case 64:
1698  image = new GFitsImageLongLong;
1699  break;
1700  case -32:
1701  image = new GFitsImageFloat;
1702  break;
1703  case -64:
1704  image = new GFitsImageDouble;
1705  break;
1706  default:
1707  std::string msg = "Invalid number of bits per pixel (bitpix="+
1708  gammalib::str(bitpix)+").";
1710  break;
1711  }
1712 
1713  // Return image pointer
1714  return image;
1715 }
1716 
1717 
1718 /***********************************************************************//**
1719  * @brief Return minimal primary HDU
1720  *
1721  * Creates a primary HDU in memory and open it using the GFitsHDU::open()
1722  * method.
1723  ***************************************************************************/
1725 {
1726  // Allocate an empty image
1728 
1729  // Create primary image in memory
1730  int status = 0;
1731  __fitsfile* fptr = NULL;
1732  status = __ffinit(&fptr, "mem://", &status);
1733  status = __ffcrim(fptr, 8, 0, NULL, &status);
1734 
1735  // Open HDU
1736  image->open(fptr,0);
1737 
1738  // Close FITS file in memory
1739  status = __ffclos(fptr, &status);
1740  if (status == 252) {
1741  status = 0;
1742  status = __ffdelt(fptr, &status);
1743  }
1744 
1745  // Initialise FITS file pointer
1746  FPTR(image->m_fitsfile)->HDUposition = 0;
1747  FPTR(image->m_fitsfile)->Fptr = NULL;
1748 
1749  // Return image
1750  return image;
1751 }
1752 
1753 
1754 /*==========================================================================
1755  = =
1756  = FITS utility functions =
1757  = =
1758  ==========================================================================*/
1759 
1760 /***********************************************************************//**
1761  * @brief Move to FITS extension
1762  *
1763  * @param[in] caller Name of caller.
1764  * @param[in] vptr FITS file void pointer.
1765  * @param[in] hdunum HDU number (optional)
1766  *
1767  * @exception GException::fits_error
1768  * cfitsio error occured.
1769  *
1770  * If @p hdunum is >0, moves the FITS file void pointer to the HDU
1771  * specified by @p hdunum. Otherwise, the FITS file void pointer is moved
1772  * to the HDU specified by the @p HDUposition attribute of the void pointer.
1773  ***************************************************************************/
1774 int gammalib::fits_move_to_hdu(const std::string& caller, void* vptr,
1775  const int& hdunum)
1776 {
1777  // Initialise status and HDU type
1778  int status = 0;
1779  int type = 0;
1780 
1781  // Set HDU position
1782  int position = (hdunum > 0) ? hdunum : (FPTR(vptr)->HDUposition)+1;
1783 
1784  // Move to HDU
1785  status = __ffmahd(FPTR(vptr), position, &type, &status);
1786 
1787  // Throw exception in case of an error
1788  if (status != 0) {
1789  std::string msg = "Unable to move FITS file pointer to extension "
1790  "number "+gammalib::str(position-1)+".";
1791  throw GException::fits_error(caller, status, msg);
1792  }
1793 
1794  // Return HDU type
1795  return type;
1796 }
Long long integer FITS image class definition.
std::vector< GFitsHDU * > m_hdu
Pointers to HDUs.
Definition: GFits.hpp:143
GFitsImage * new_image(void)
Allocate new FITS image and return memory pointer.
Definition: GFits.cpp:1642
GFitsHDU * at(const int &extno)
Get pointer to HDU.
Definition: GFits.cpp:233
VO client class.
Definition: GVOClient.hpp:59
Double precision FITS image class definition.
Abstract FITS image base class.
Definition: GFitsImage.hpp:43
#define __ffdhdu(A, B, C)
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition: GFits.cpp:482
GFits(void)
Void constructor.
Definition: GFits.cpp:90
#define G_INSERT1
Definition: GFits.cpp:59
#define G_OPEN
Definition: GFits.cpp:63
#define FPTR(A)
void clear(void)
Clear FITS file.
Definition: GFits.cpp:195
Signed Byte FITS image class definition.
void publish(const GFitsHDU &hdu)
Publish FITS HDU.
Definition: GVOClient.cpp:372
void open(void *vptr, int hdunum)
Open HDU.
Definition: GFitsHDU.cpp:285
Single precision FITS image class definition.
Abstract FITS extension base class.
Definition: GFitsHDU.hpp:51
Signed Byte FITS image class.
void connect(void *fptr)
Connect HDU to FITS file.
Definition: GFitsHDU.cpp:192
#define __ffflus(A, B)
#define G_INSERT2
Definition: GFits.cpp:60
const GFilename & filename(void) const
Return FITS filename.
Definition: GFits.hpp:313
Unsigned long image FITS class.
virtual HDUType exttype(void) const =0
GFilename m_filename
FITS file name.
Definition: GFits.hpp:144
#define __ffinit(A, B, C)
#define G_PUBLISH
Definition: GFits.cpp:66
void * m_fitsfile
FITS file pointer pointing on actual HDU.
Definition: GFitsHDU.hpp:133
#define FHANDLE(A)
Unsigned short FITS image class definition.
Gammalib tools definition.
#define __ffcrim(A, B, C, D, E)
#define __ffthdu(A, B, C)
FITS file class.
Definition: GFits.hpp:63
FITS file class interface definition.
Long integer FITS image class definition.
#define __ffclos(A, B)
const int & extno(void) const
Return extension number.
Definition: GFitsHDU.hpp:176
GFits & operator=(const GFits &fits)
Assignment operator.
Definition: GFits.cpp:165
Long long integer FITS image class.
VO client class interface definition.
void saveto(const GFilename &filename, const bool &clobber=false)
Saves to specified FITS file.
Definition: GFits.cpp:1293
Unsigned long image FITS class definition.
GFitsImage * new_primary(void)
Return minimal primary HDU.
Definition: GFits.cpp:1724
void save(const bool &clobber=false)
Saves FITS file.
Definition: GFits.cpp:1178
#define __ffmahd(A, B, C, D)
void remove(const int &extno)
Remove HDU from FITS file.
Definition: GFits.cpp:862
Double precision FITS image class.
GFits * clone(void) const
Clone FITS file.
Definition: GFits.cpp:213
#define G_TABLE1
Definition: GFits.cpp:55
CFITSIO interface header.
Filename class.
Definition: GFilename.hpp:62
void open(const GFilename &filename, const bool &create=false)
Open or (optionally) create FITS file.
Definition: GFits.cpp:1041
void publish(const int &extno, const std::string &name="") const
Publish FITS HDU on VO Hub.
Definition: GFits.cpp:1363
void free_members(void)
Delete class members.
Definition: GFits.cpp:1565
#define __ffopen(A, B, C, D)
bool exists(void) const
Checks whether file exists.
Definition: GFilename.cpp:223
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
virtual GFitsHDU * clone(void) const =0
Clones object.
#define G_NEW_IMAGE
Definition: GFits.cpp:68
#define __ffgky(A, B, C, D, E, F)
GChatter
Definition: GTypemaps.hpp:33
#define G_SAVE
Definition: GFits.cpp:64
#define __ffdelt(A, B)
Short integer FITS image class.
FITS ASCII table class.
#define G_REMOVE1
Definition: GFits.cpp:61
Single precision FITS image class.
Unsigned short FITS image class.
void remove(void) const
Remove file from disk.
Definition: GFilename.cpp:345
const std::string & extname(void) const
Return extension name.
Definition: GFitsHDU.hpp:162
Short integer FITS image class definition.
#define G_REMOVE2
Definition: GFits.cpp:62
FITS Byte image class definition.
#define G_AT2
Definition: GFits.cpp:52
Long integer FITS image class.
GFitsImage * image(const int &extno)
Get pointer to image HDU.
Definition: GFits.cpp:368
#define G_IMAGE1
Definition: GFits.cpp:53
void reserve(const int &num)
Reserves space for HDUs in FITS file.
Definition: GFits.hpp:265
std::string url(void) const
Return Uniform Resource Locator (URL)
Definition: GFilename.hpp:189
void init_members(void)
Initialise class members.
Definition: GFits.cpp:1510
#define __TLONGLONG
void clear(void)
Clear file name.
Definition: GFilename.cpp:188
bool is_empty(void) const
Signals if there are no HDUs in FITS file.
Definition: GFits.hpp:251
#define G_SAVETO
Definition: GFits.cpp:65
#define G_TABLE2
Definition: GFits.cpp:56
FITS binary table class.
int size(void) const
Return number of HDUs in FITS file.
Definition: GFits.hpp:237
GChatter reduce(const GChatter &chatter)
Reduce chattiness by one level.
Definition: GTypemaps.hpp:65
FITS Byte image class.
GFitsHDU * insert(const int &extno, const GFitsHDU &hdu)
Set HDU for the specified extension number.
Definition: GFits.cpp:778
virtual ~GFits(void)
Destructor.
Definition: GFits.cpp:143
bool m_created
FITS file has been created (true/false)
Definition: GFits.hpp:147
Exception handler interface definition.
bool m_readwrite
FITS file is readwrite (true/false)
Definition: GFits.hpp:146
GFitsHDU * append(const GFitsHDU &hdu)
Append HDU to FITS file.
Definition: GFits.cpp:678
std::string toupper(const std::string &s)
Convert string to upper case.
Definition: GTools.cpp:941
VO table class definition.
FITS binary table class definition.
void copy_members(const GFits &fits)
Copy class members.
Definition: GFits.cpp:1535
#define G_AT1
Definition: GFits.cpp:51
int fits_move_to_hdu(const std::string &caller, void *vptr, const int &hdunum=0)
Move to FITS extension.
Definition: GFits.cpp:1774
#define G_FREE_MEM
Definition: GFits.cpp:67
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
#define G_SET2
Definition: GFits.cpp:58
VOTable class.
Definition: GVOTable.hpp:55
GFitsHDU * set(const int &extno, const GFitsHDU &hdu)
Set HDU for the specified extension number.
Definition: GFits.cpp:597
FITS ASCII table class definition.
std::string print(const GChatter &chatter=NORMAL) const
Print FITS information.
Definition: GFits.cpp:1455
#define G_SET1
Definition: GFits.cpp:57
void close(void)
Close FITS file.
Definition: GFits.cpp:1342
int extno(const std::string &extname) const
Return extension number for specified extension name.
Definition: GFits.cpp:990
void * m_fitsfile
FITS file pointer.
Definition: GFits.hpp:145
#define G_IMAGE2
Definition: GFits.cpp:54
#define __ffgiet(A, B, C)
void extend(const GFits &fits)
Append FITS file.
Definition: GFits.cpp:934
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489