RaDecByAst.hxx

Go to the documentation of this file.
00001 
00002 
00003 
00004 //-----------------------------------------------------------------------------
00005 // CVS: $Revision: 6830 $, $Date: 2010-04-16 15:33:45 +0100 (Fri, 16 Apr 2010) $, $Author: EckhardSutorius $
00006 //-----------------------------------------------------------------------------
00007 #ifndef RADECBYAST_HXX
00008 #define RADECBYAST_HXX
00009 
00010 // We will make use of AST for operations involving WCS
00011 extern "C"
00012 {
00013   #include "ast.h"
00014 }
00015 
00016 #include <cmath>
00017 
00018 #include "Converter.h"
00019 #include "FitsFile.h"
00020 #include "Logger.h"
00021 #include "DataMethod.hxx"
00022 #include "StringOps.h"
00023 //-----------------------------------------------------------------------------
00026 template<typename DataType>
00027 class RaDecByAst : public DataMethod<DataType>
00028 {
00029 public:
00031   RaDecByAst(const FitsFile& aFile)
00032     : mFile(aFile)
00033     { DataMethod<DataType>::methodName = "RaDecByAst"; }
00034 
00036   virtual ~RaDecByAst() { }
00037 
00040   void setSource(IntMap& tbls, IntMap& cols, StringMap units)
00041   {
00042     extnumColno = cols[0];
00043     xColno = cols[1];
00044     yColno = cols[2];
00045     extnumTblno = tbls[0];
00046     xTblno = tbls[1];
00047     yTblno = tbls[2];
00048   }
00049 
00052   void setTarget(IntMap& cols, StringMap units)
00053   {
00054     raColno  = cols[0];
00055     decColno = cols[1];
00056 
00057     // AST will calculate angles in radians, so we need the conversion
00058     // factors for the required units in the schema
00059     Converter conv;
00060     raFact  = conv.getFactor("RADIANS", units[0]);
00061     decFact = conv.getFactor("RADIANS", units[1]);
00062   }
00063 
00065   void doit2(TableData<DataType>& inpData, TableData<DataType>& data, int row1, int row2)
00066   {
00067     double pi = M_PI;
00068     string filename = mFile.getFileName();
00069     int numHDUs = mFile.getNumHdus();
00070     double crval1[numHDUs], crval2[numHDUs];
00071 
00072     // Set up AST
00073     astBegin;
00074 
00075     // Create a new AST fits channel
00076     AstFitsChan *fitschan = astFitsChan(NULL, NULL, "%.*s", 0, "NULL");
00077     AstFrameSet *wcsinfo0 = 0;
00078     AstFrameSet *wcsinfo1 = 0;
00079     AstFrameSet *wcsinfo2 = 0;
00080     AstFrameSet *wcsinfo3 = 0;
00081     AstFrameSet *wcsinfo4 = 0;
00082     AstFrameSet *wcsinfo5 = 0;
00083     AstFrameSet *wcsinfo6 = 0;
00084     AstFrameSet *wcsinfo7 = 0;
00085     AstFrameSet *wcsinfo8 = 0;
00086     AstFrameSet *wcsinfo9 = 0;
00087     AstFrameSet *wcsinfo10 = 0;
00088     AstFrameSet *wcsinfo11 = 0;
00089     AstFrameSet *wcsinfo12 = 0;
00090     AstFrameSet *wcsinfo13 = 0;
00091     AstFrameSet *wcsinfo14 = 0;
00092     AstFrameSet *wcsinfo15 = 0;
00093 
00094     for (int hduNo = 2; hduNo <= numHDUs; ++hduNo)
00095     {
00096       int ext = hduNo - 2;
00097 
00098       // Move to extension on the FITS file
00099       mFile.movabsHdu(hduNo);
00100 
00101       // Put all cards in the header of this extension into the AST
00102       // fits channel
00103       int ncards = mFile.getNumCards();
00104 
00105       for (int n = 0; n < ncards; ++n) {
00106         char card[80];
00107         mFile.readCard(n, card);
00108         string cardStr = card;
00109         if (cardStr.find("TCTYP3") == 0) {
00110           cardStr.replace(0, 6, "CTYPE1");
00111         } else if (cardStr.find("TCTYP5") == 0) {
00112           cardStr.replace(0, 6, "CTYPE2");
00113         } else if (cardStr.find("TCRVL3") == 0) {
00114           cardStr.replace(0, 6, "CRVAL1");
00115         } else if (cardStr.find("TCRVL5") == 0) {
00116           cardStr.replace(0, 6, "CRVAL2");
00117         } else if (cardStr.find("TCRPX3") == 0) {
00118           cardStr.replace(0, 6, "CRPIX1");
00119         } else if (cardStr.find("TCRPX5") == 0) {
00120           cardStr.replace(0, 6, "CRPIX2");
00121         } else if (cardStr.find("TC3_3") == 0) {
00122           cardStr.replace(0, 5, "CD1_1");
00123         } else if (cardStr.find("TC3_5") == 0) {
00124           cardStr.replace(0, 5, "CD1_2");
00125         } else if (cardStr.find("TC5_3") == 0) {
00126           cardStr.replace(0, 5, "CD2_1");
00127         } else if (cardStr.find("TC5_5") == 0) {
00128           cardStr.replace(0, 5, "CD2_2");
00129         } else if (cardStr.find("TV5_1") == 0) {
00130           cardStr.replace(0, 5, "PV2_1");
00131         } else if (cardStr.find("TV5_2") == 0) {
00132           cardStr.replace(0, 5, "PV2_2");
00133         } else if (cardStr.find("TV5_3") == 0) {
00134           cardStr.replace(0, 5, "PV2_3");
00135         } else if (cardStr.find("TV5_4") == 0) {
00136           cardStr.replace(0, 5, "PV2_4");
00137         } else if (cardStr.find("TV5_5") == 0) {
00138           cardStr.replace(0, 5, "PV2_5");
00139         }
00140         astPutFits(fitschan, cardStr.c_str(), 0);
00141       }
00142       // Make AST wcsinfo object
00143       astClear( fitschan, "Card" );
00144 
00145       mFile.readKey(FitsFile::cCRVAL1, crval1[ext]);
00146       mFile.readKey(FitsFile::cCRVAL2, crval2[ext]);
00147 
00148       if (crval1[ext] == 0 && crval2[ext] == 0) {
00149           log.addError("Could not read WCS info in RDBA for " + filename);
00150       } else {
00151         if (ext == 0) {
00152           wcsinfo0 = (AstFrameSet*)astRead(fitschan);
00153         } else if (ext == 1) {
00154           wcsinfo1 = (AstFrameSet*)astRead(fitschan);
00155         } else if (ext == 2) {
00156           wcsinfo2 = (AstFrameSet*)astRead(fitschan);
00157         } else if (ext == 3) {
00158           wcsinfo3 = (AstFrameSet*)astRead(fitschan);
00159         } else if (ext == 4) {
00160           wcsinfo4 = (AstFrameSet*)astRead(fitschan);
00161         } else if (ext == 5) {
00162           wcsinfo5 = (AstFrameSet*)astRead(fitschan);
00163         } else if (ext == 6) {
00164           wcsinfo6 = (AstFrameSet*)astRead(fitschan);
00165         } else if (ext == 7) {
00166           wcsinfo7 = (AstFrameSet*)astRead(fitschan);
00167         } else if (ext == 8) {
00168           wcsinfo8 = (AstFrameSet*)astRead(fitschan);
00169         } else if (ext == 9) {
00170           wcsinfo9 = (AstFrameSet*)astRead(fitschan);
00171         } else if (ext == 10) {
00172           wcsinfo10 = (AstFrameSet*)astRead(fitschan);
00173         } else if (ext == 11) {
00174           wcsinfo11 = (AstFrameSet*)astRead(fitschan);
00175         } else if (ext == 12) {
00176           wcsinfo12 = (AstFrameSet*)astRead(fitschan);
00177         } else if (ext == 13) {
00178           wcsinfo13 = (AstFrameSet*)astRead(fitschan);
00179         } else if (ext == 14) {
00180           wcsinfo14 = (AstFrameSet*)astRead(fitschan);
00181         } else if (ext == 15) {
00182           wcsinfo15 = (AstFrameSet*)astRead(fitschan);
00183         }
00184       }
00185     }
00186 
00187     // Use these arrays for AST input
00188     double xpix[1], ypix[1], xsky[1], ysky[1];
00189 
00190     for (int rowNo = row1; rowNo < row2; ++rowNo)
00191     {
00192       // Get the FITS extension number for the data in this row. Note
00193       // that the extension number is assumed to be in the fitsio
00194       // convention
00195       Numeric value;
00196       data.value(extnumColno, rowNo, value);
00197       int extNum = (int)value;
00198 
00199       if (crval1[extNum] == 0 && crval2[extNum] == 0) {
00200         log.addError("No WCS info for " + filename
00201                      + "[" + StringOps::NumToString(extNum) + "]");
00202       } else {
00203         // Input X and Y
00204         //double xpix, ypix;
00205         inpData.value(xColno, rowNo, xpix[0]);
00206         inpData.value(yColno, rowNo, ypix[0]);
00207 
00208         // Get the sky coordinates
00209         // Note these will be in radians. Angles can also be negative.
00210         if (extNum == 2) {
00211           astTran2(wcsinfo0, 1, xpix, ypix, 1, xsky, ysky);
00212         } else if (extNum == 3) {
00213           astTran2(wcsinfo1, 1, xpix, ypix, 1, xsky, ysky);
00214         } else if (extNum == 4) {
00215           astTran2(wcsinfo2, 1, xpix, ypix, 1, xsky, ysky);
00216         } else if (extNum == 5) {
00217           astTran2(wcsinfo3, 1, xpix, ypix, 1, xsky, ysky);
00218         } else if (extNum == 6) {
00219           astTran2(wcsinfo4, 1, xpix, ypix, 1, xsky, ysky);
00220         } else if (extNum == 7) {
00221           astTran2(wcsinfo5, 1, xpix, ypix, 1, xsky, ysky);
00222         } else if (extNum == 8) {
00223           astTran2(wcsinfo6, 1, xpix, ypix, 1, xsky, ysky);
00224         } else if (extNum == 9) {
00225           astTran2(wcsinfo7, 1, xpix, ypix, 1, xsky, ysky);
00226         } else if (extNum == 10) {
00227           astTran2(wcsinfo8, 1, xpix, ypix, 1, xsky, ysky);
00228         } else if (extNum == 11) {
00229           astTran2(wcsinfo9, 1, xpix, ypix, 1, xsky, ysky);
00230         } else if (extNum == 12) {
00231           astTran2(wcsinfo10, 1, xpix, ypix, 1, xsky, ysky);
00232         } else if (extNum == 13) {
00233           astTran2(wcsinfo11, 1, xpix, ypix, 1, xsky, ysky);
00234         } else if (extNum == 14) {
00235           astTran2(wcsinfo12, 1, xpix, ypix, 1, xsky, ysky);
00236         } else if (extNum == 15) {
00237           astTran2(wcsinfo13, 1, xpix, ypix, 1, xsky, ysky);
00238         } else if (extNum == 16) {
00239           astTran2(wcsinfo14, 1, xpix, ypix, 1, xsky, ysky);
00240         } else if (extNum == 17) {
00241           astTran2(wcsinfo15, 1, xpix, ypix, 1, xsky, ysky);
00242         }
00243         // Make sure all RAs are in the range 0 to 2*pi
00244         if (xsky[0] < 0) xsky[0] += 2 * pi;
00245 
00246         Numeric ra = raFact * xsky[0];
00247         Numeric dec = decFact * ysky[0];
00248 
00249         data.assign(raColno, rowNo, ra);
00250         data.assign(decColno, rowNo, dec);
00251       }
00252     }
00253     // Finshed with AST
00254     astEnd;
00255   }
00256 
00257   // No simple data model available
00258   void doit(TableData<DataType>& data, int row1, int row2) { }
00259 
00260 private:
00261   FitsFile mFile;
00262   int extnumColno, raColno, decColno, xColno, yColno;
00263   int extnumTblno, xTblno, yTblno;
00264   double decFact, raFact;
00265   Logger log;
00266 };
00267 //-----------------------------------------------------------------------------
00268 #endif
00269 //-----------------------------------------------------------------------------
00270 // Change log:
00271 //
00272 //  7-Jun-2005, ETWS: Original version.
00273 //  7-Apr-2008, ETWS: Upgraded to use new detection table layout.