// ----------------------------------------------------------------------- // headerIO.cc: Read in various PHU keywords from a FITS file. // ----------------------------------------------------------------------- // Copyright (C) 2006, Matthew Whiting, ATNF // // This program is free software; you can redistribute it and/or modify it // under the terms of the GNU General Public License as published by the // Free Software Foundation; either version 2 of the License, or (at your // option) any later version. // // Duchamp is distributed in the hope that it will be useful, but WITHOUT // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or // FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License // for more details. // // You should have received a copy of the GNU General Public License // along with Duchamp; if not, write to the Free Software Foundation, // Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA // // Correspondence concerning Duchamp may be directed to: // Internet email: Matthew.Whiting [at] atnf.csiro.au // Postal address: Dr. Matthew Whiting // Australia Telescope National Facility, CSIRO // PO Box 76 // Epping NSW 1710 // AUSTRALIA // ----------------------------------------------------------------------- #include #include #include #include #define WCSLIB_GETWCSTAB // define this so that we don't try and redefine // wtbarr (this is a problem when using gcc v.4+ #include #include #include #include #include namespace duchamp { OUTCOME FitsHeader::readHeaderInfo(std::string fname, Param &par) { OUTCOME returnValue = SUCCESS; // Open the FITS file fitsfile *fptr; int status = 0; if( fits_open_file(&fptr,fname.c_str(),READONLY,&status) ){ DUCHAMPWARN("Cube Reader","Error opening file "<readHeaderInfo(fptr,par); // Close the FITS file. status = 0; fits_close_file(fptr, &status); if (status){ DUCHAMPWARN("Cube Reader","Error closing file: "); fits_report_error(stderr, status); } } return returnValue; } OUTCOME FitsHeader::readHeaderInfo(Param &par) { OUTCOME returnValue = SUCCESS; if(this->fptr == 0) { returnValue = FAILURE; DUCHAMPERROR("Cube Reader", "FITS file not opened."); } else { returnValue = this->readHeaderInfo(this->fptr, par); } return returnValue; } OUTCOME FitsHeader::readHeaderInfo(fitsfile *fptr, Param &par) // OUTCOME FitsHeader::readHeaderInfo(std::string fname, Param &par) { /// A simple front-end function to the three header access /// functions defined below. OUTCOME returnValue = SUCCESS; // Get the brightness unit, so that we can set the units for the // integrated flux when we go to fixSpectralUnits. OUTCOME bunitResult = this->readBUNIT(fptr); // if(this->readBLANKinfo(fname, par)==FAILURE) OUTCOME blankResult = this->readBLANKinfo(fptr, par); OUTCOME beamResult = this->readBeamInfo(fptr, par); if(this->wcs->spec>=0) this->fixSpectralUnits(par.getSpectralUnits()); if(bunitResult == FAILURE || blankResult==FAILURE || beamResult==FAILURE ){ DUCHAMPWARN("Cube Reader","Header could not be read completely"); returnValue = FAILURE; } return returnValue; } ////////////////////////////////////////////////// OUTCOME FitsHeader::readBUNIT(fitsfile *fptr) { /// Read the BUNIT header keyword, to store the units of brightness (flux). /// \param fname The name of the FITS file. char *comment = new char[FLEN_COMMENT]; char *unit = new char[FLEN_VALUE]; OUTCOME returnStatus=SUCCESS; int status = 0; // Read the BUNIT keyword, and translate to standard unit format if needs be std::string header("BUNIT"); fits_read_key(fptr, TSTRING, (char *)header.c_str(), unit, comment, &status); if (status){ DUCHAMPWARN("Cube Reader","Error reading BUNIT keyword: "); fits_report_error(stderr, status); returnStatus = FAILURE; } else{ wcsutrn(0,unit); this->fluxUnits = unit; this->originalFluxUnits = unit; returnStatus = SUCCESS; } delete [] comment; delete [] unit; return returnStatus; } ////////////////////////////////////////////////// OUTCOME FitsHeader::readBLANKinfo(fitsfile *fptr, Param &par) { /// Reading in the Blank pixel value keywords, which is only done /// if requested via the flagBlankPix parameter. /// /// If the BLANK keyword is in the header, use that and store the relevant /// values. Also copy them into the parameter set. /// If not, use the default value (either the default from param.cc or /// from the param file) and assume simple values for the keywords: ///
  • The scale keyword is the same as the blank value, ///
  • The blank keyword (which is an int) is 1 and ///
  • The bzero (offset) is 0. ///
/// \param fname The name of the FITS file. /// \param par The Param set: to know the flagBlankPix value and to /// store the keywords. OUTCOME returnStatus = SUCCESS; int status = 0; char *comment = new char[FLEN_COMMENT]; int blank; float bscale, bzero; // Read the BLANK keyword. // If this isn't present, make sure flagTrim is false (if it is // currently true, let the user know you're doing this) and set // flagBlankPix to false so that the isBlank functions all return // false // If it is, read the other two necessary keywords, and then set // the values accordingly. std::string header("BLANK"); if(fits_read_key(fptr, TINT, (char *)header.c_str(), &blank, comment, &status)){ par.setFlagBlankPix(false); if(par.getFlagTrim()){ par.setFlagTrim(false); if(status == KEY_NO_EXIST){ DUCHAMPWARN("Cube Reader", "There is no BLANK keyword present. Not doing any trimming."); } else{ DUCHAMPWARN("Cube Reader", "Error reading BLANK keyword, so not doing any trimming."); fits_report_error(stderr, status); } returnStatus = FAILURE; //only return a failure if we care whether the BLANK keyword is present, which is only if the Trim flag is set. } } else{ status = 0; header="BZERO"; fits_read_key(fptr, TFLOAT, (char *)header.c_str(), &bzero, comment, &status); status = 0; header="BSCALE"; fits_read_key(fptr, TFLOAT, (char *)header.c_str(), &bscale, NULL, &status); this->blankKeyword = blank; this->bscaleKeyword = bscale; this->bzeroKeyword = bzero; par.setFlagBlankPix(true); par.setBlankKeyword(blank); par.setBscaleKeyword(bscale); par.setBzeroKeyword(bzero); par.setBlankPixVal( blank*bscale + bzero ); } delete [] comment; return returnStatus; } ////////////////////////////////////////////////// OUTCOME FitsHeader::readBeamInfo(fitsfile *fptr, Param &par) { /// Reading in the beam parameters from the header. /// Use these, plus the basic WCS parameters to calculate the size of /// the beam in pixels. Copy the beam size into the parameter set. /// If information not present in FITS header, use the parameter /// set to define the beam size. /// \param fname The name of the FITS file. /// \param par The Param set. this->itsBeam.readFromFITS(fptr, par, this->getAvPixScale()); return SUCCESS; } }