[160] | 1 | #include <iostream> |
---|
| 2 | #include <string> |
---|
| 3 | #include <wcs.h> |
---|
| 4 | #include <wcshdr.h> |
---|
| 5 | #include <fitshdr.h> |
---|
| 6 | #include <wcsfix.h> |
---|
| 7 | #include <wcsunits.h> |
---|
| 8 | #define WCSLIB_GETWCSTAB // define this so that we don't try and redefine |
---|
| 9 | // wtbarr (this is a problem when using gcc v.4+ |
---|
| 10 | #include <fitsio.h> |
---|
| 11 | #include <math.h> |
---|
| 12 | #include <duchamp.hh> |
---|
| 13 | #include <Cubes/cubes.hh> |
---|
| 14 | |
---|
| 15 | int FitsHeader::defineWCS(string fname, Param &par) |
---|
| 16 | { |
---|
| 17 | /** |
---|
| 18 | * FitsHeader::defineWCS(wcsprm *wcs, char *hdr, Param &par) |
---|
| 19 | * A function that reads the WCS header information from the |
---|
| 20 | * FITS file given by fptr. |
---|
| 21 | * It will also sort out the spectral axis, and covert to the correct |
---|
| 22 | * velocity type, or frequency type if need be. |
---|
[161] | 23 | * It calls FitsHeader::getBUNIT so that the Integrated Flux units can |
---|
| 24 | * be calculated by FitsHeader::fixUnits. |
---|
[160] | 25 | */ |
---|
| 26 | |
---|
| 27 | fitsfile *fptr; |
---|
| 28 | int bitpix,numAxes; |
---|
| 29 | int noComments = 1; //so that fits_hdr2str will ignore COMMENT, HISTORY etc |
---|
| 30 | int nExc = 0,nkeys; |
---|
| 31 | char *hdr; |
---|
[161] | 32 | |
---|
| 33 | // Open the FITS file |
---|
[160] | 34 | int status = 0; |
---|
| 35 | if( fits_open_file(&fptr,fname.c_str(),READONLY,&status) ){ |
---|
| 36 | fits_report_error(stderr, status); |
---|
| 37 | return FAILURE; |
---|
| 38 | } |
---|
[161] | 39 | |
---|
| 40 | // Get the dimensions of the FITS file -- number of axes and size of each. |
---|
[160] | 41 | status = 0; |
---|
| 42 | if(fits_get_img_dim(fptr, &numAxes, &status)){ |
---|
| 43 | fits_report_error(stderr, status); |
---|
| 44 | return FAILURE; |
---|
| 45 | } |
---|
| 46 | long *dimAxes = new long[numAxes]; |
---|
| 47 | for(int i=0;i<numAxes;i++) dimAxes[i]=1; |
---|
| 48 | status = 0; |
---|
| 49 | if(fits_get_img_size(fptr, numAxes, dimAxes, &status)){ |
---|
| 50 | fits_report_error(stderr, status); |
---|
| 51 | return FAILURE; |
---|
| 52 | } |
---|
[161] | 53 | |
---|
| 54 | // Read in the entire PHU of the FITS file to a string. |
---|
| 55 | // This will be read by the wcslib functions to extract the WCS. |
---|
[160] | 56 | status = 0; |
---|
| 57 | fits_hdr2str(fptr, noComments, NULL, nExc, &hdr, &nkeys, &status); |
---|
| 58 | if( status ){ |
---|
| 59 | duchampWarning("defineWCS","Error whilst reading FITS header to string: "); |
---|
| 60 | fits_report_error(stderr, status); |
---|
| 61 | return FAILURE; |
---|
| 62 | } |
---|
[161] | 63 | |
---|
| 64 | // Close the FITS file -- not needed any more in this function. |
---|
[160] | 65 | status = 0; |
---|
| 66 | fits_close_file(fptr, &status); |
---|
| 67 | if (status){ |
---|
[168] | 68 | duchampWarning("defineWCS","Error closing file: "); |
---|
[160] | 69 | fits_report_error(stderr, status); |
---|
| 70 | } |
---|
| 71 | |
---|
| 72 | wcsprm *wcs = new wcsprm; |
---|
[161] | 73 | |
---|
| 74 | // Initialise the wcsprm structure |
---|
[163] | 75 | int flag; |
---|
[160] | 76 | if(flag = wcsini(true,numAxes,wcs)){ |
---|
| 77 | std::stringstream errmsg; |
---|
| 78 | errmsg << "wcsini failed! Code=" << flag |
---|
| 79 | << ": " << wcs_errmsg[flag] << std::endl; |
---|
| 80 | duchampError("defineWCS",errmsg.str()); |
---|
| 81 | return FAILURE; |
---|
| 82 | } |
---|
| 83 | wcs->flag=-1; |
---|
[161] | 84 | |
---|
| 85 | int relax=1; // for wcspih -- admit all recognised informal WCS extensions |
---|
| 86 | int ctrl=2; // for wcspih -- report each rejected card and its reason for |
---|
| 87 | // rejection |
---|
[163] | 88 | int nwcs, nreject; |
---|
[161] | 89 | // Parse the FITS header to fill in the wcsprm structure |
---|
[160] | 90 | if(flag = wcspih(hdr, nkeys, relax, ctrl, &nreject, &nwcs, &wcs)) { |
---|
[161] | 91 | // if here, something went wrong -- report what. |
---|
[160] | 92 | std::stringstream errmsg; |
---|
| 93 | errmsg << "wcspih failed! Code="<<flag<<": "<<wcs_errmsg[flag]<<std::endl; |
---|
| 94 | duchampWarning("defineWCS",errmsg.str()); |
---|
| 95 | } |
---|
| 96 | else{ |
---|
| 97 | int stat[6]; |
---|
| 98 | int axes[3]={dimAxes[wcs->lng],dimAxes[wcs->lat],dimAxes[wcs->spec]}; |
---|
[161] | 99 | |
---|
| 100 | // Applies all necessary corrections to the wcsprm structure |
---|
| 101 | // (missing cards, non-standard units or spectral types, ...) |
---|
[160] | 102 | if(flag=wcsfix(1,axes,wcs,stat)) { |
---|
| 103 | std::stringstream errmsg; |
---|
| 104 | errmsg << "wcsfix failed:\n"; |
---|
| 105 | for(int i=0; i<NWCSFIX; i++) |
---|
| 106 | if (stat[i] > 0) |
---|
| 107 | errmsg <<" flag="<<flag<<": "<< wcsfix_errmsg[stat[i]]<<std::endl; |
---|
| 108 | duchampWarning("defineWCS", errmsg.str() ); |
---|
| 109 | } |
---|
[161] | 110 | |
---|
| 111 | // Set up the wcsprm struct. Report if something goes wrong. |
---|
[160] | 112 | if(flag=wcsset(wcs)){ |
---|
| 113 | std::stringstream errmsg; |
---|
| 114 | errmsg<<"wcsset failed! Code="<<flag <<": "<<wcs_errmsg[flag]<<std::endl; |
---|
| 115 | duchampWarning("defineWCS",errmsg.str()); |
---|
| 116 | } |
---|
| 117 | |
---|
| 118 | if(wcs->naxis>2){ // if there is a spectral axis |
---|
| 119 | |
---|
| 120 | int index = wcs->spec; |
---|
| 121 | string desiredType,specType = wcs->ctype[index]; |
---|
| 122 | if(wcs->restfrq != 0){ |
---|
| 123 | // Set the spectral axis to a standard specification: VELO-F2V |
---|
| 124 | desiredType = duchampVelocityType; |
---|
| 125 | if(wcs->restwav == 0) wcs->restwav = 299792458.0 / wcs->restfrq; |
---|
[183] | 126 | this->spectralDescription = duchampSpectralDescription[VELOCITY]; |
---|
[160] | 127 | } |
---|
| 128 | else{ |
---|
| 129 | // No rest frequency defined, so put spectral dimension in frequency. |
---|
| 130 | // Set the spectral axis to a standard specification: FREQ |
---|
| 131 | duchampWarning("defineWCS", |
---|
| 132 | "No rest frequency defined. Using frequency units in spectral axis.\n"); |
---|
| 133 | desiredType = duchampFrequencyType; |
---|
| 134 | par.setSpectralUnits("MHz"); |
---|
| 135 | if(strcmp(wcs->cunit[index],"")==0){ |
---|
| 136 | duchampWarning("defineWCS", |
---|
| 137 | "No frequency unit given. Assuming frequency axis is in Hz.\n"); |
---|
| 138 | strcpy(wcs->cunit[index],"Hz"); |
---|
| 139 | } |
---|
[183] | 140 | this->spectralDescription = duchampSpectralDescription[FREQUENCY]; |
---|
[160] | 141 | } |
---|
[161] | 142 | |
---|
| 143 | // Check to see if the spectral type (eg VELO-F2V) matches that wanted |
---|
| 144 | // from duchamp.hh. Only first four characters checked. |
---|
[160] | 145 | if(strncmp(specType.c_str(),desiredType.c_str(),4)!=0){ |
---|
| 146 | index = -1; |
---|
[161] | 147 | // If not a match, translate the spectral axis to the desired type |
---|
[160] | 148 | if( flag = wcssptr(wcs, &index, (char *)desiredType.c_str())){ |
---|
| 149 | std::stringstream errmsg; |
---|
| 150 | errmsg<<"wcssptr failed! Code="<<flag <<": " |
---|
| 151 | <<wcs_errmsg[flag]<<std::endl; |
---|
| 152 | duchampWarning("defineWCS",errmsg.str()); |
---|
| 153 | } |
---|
| 154 | } |
---|
| 155 | |
---|
| 156 | } // end of if(numAxes>2) |
---|
| 157 | |
---|
[161] | 158 | // Save the wcs to the FitsHeader class that is running this function |
---|
[160] | 159 | this->setWCS(wcs); |
---|
| 160 | this->setNWCS(nwcs); |
---|
| 161 | |
---|
| 162 | wcsfree(wcs); |
---|
| 163 | |
---|
| 164 | } |
---|
| 165 | |
---|
[161] | 166 | // Get the brightness unit, so that we can set the units for the |
---|
| 167 | // integrated flux when we go to fixUnits. |
---|
| 168 | |
---|
[160] | 169 | this->getBUNIT(fname); |
---|
| 170 | |
---|
| 171 | this->fixUnits(par); |
---|
| 172 | |
---|
| 173 | return SUCCESS; |
---|
| 174 | |
---|
| 175 | } |
---|