Magnitudes.hxx

Go to the documentation of this file.
00001 
00002 
00003 
00004 //-----------------------------------------------------------------------------
00005 // SVN: $Id: Magnitudes.hxx 7080 2010-06-10 15:34:33Z EckhardSutorius $
00006 //-----------------------------------------------------------------------------
00007 #ifndef MAGNITUDES_H
00008 #define MAGNITUDES_H
00009 
00010 #include <float.h>
00011 
00012 #include "Constants.h"
00013 #include "DataMethod.hxx"
00014 #include "FitsFile.h"
00015 #include "StringOps.h"
00016 #include "IllumTable.hxx"
00017 #include "Options.h"
00018 //-----------------------------------------------------------------------------
00024 template<typename DataType>
00025 class CalibratedMagnitude : public DataMethod<DataType>
00026 {
00027 public:
00032   CalibratedMagnitude(const FitsFile& aCatFile, IllumTable& anIllum)
00033     : mCatFile(aCatFile)
00034     { 
00035       DataMethod<DataType>::methodName = "CalibratedMagnitude"; 
00036       mIllum = &anIllum;
00037     }
00038 
00040   virtual ~CalibratedMagnitude() { }
00041 
00043   void setSource(IntMap& tbls, IntMap& cols, StringMap units)
00044   {
00045     mExtnumColno = cols[0];
00046     mFluxColno = cols[1];
00047     mXColno = cols[2];
00048     mYColno = cols[3];
00049     mExtnumTblno = tbls[0];
00050     mFluxTblno = tbls[1];
00051     mXTblno = tbls[2];
00052     mYTblno = tbls[3];
00053     if (mFluxColno < 0)
00054     {
00055       throw DataMethodException("CalibratedMagnitude method: cannot source");
00056     }
00057   }
00058 
00060   void setTarget(IntMap& cols, StringMap units) { mMagColno = cols[0]; }
00061 
00066   void doit2(TableData<DataType>& inpData, TableData<DataType>& data, int row1, int row2)
00067   {
00068     Options options;
00069     double pi = M_PI;
00070     double flux2mag = 2.5/M_LN10;
00071 
00072     int numHDUs = mCatFile.getNumHdus();
00073     double pamss, pamse, amss[numHDUs], amse[numHDUs];
00074     double zpmag[numHDUs], zperr[numHDUs];
00075     double extinct[numHDUs];
00076     double pexptime, exptime[numHDUs];
00077     std::string ctype1[numHDUs],ctype2[numHDUs];
00078     double crpix1[numHDUs],crpix2[numHDUs];
00079     double cd11[numHDUs],cd12[numHDUs],cd21[numHDUs],cd22[numHDUs];
00080     double pv21[numHDUs],pv23[numHDUs];
00081 
00082     // get frame type
00083     std::string frameType = DataOps::getFrameType(mCatFile.getFileName());
00084 
00085     if (StringOps::ToLower(options.getArchive()) == "vsa")
00086     {
00087       mCatFile.movabsHdu(1);
00088       // Read the airmass values from primary header for VSA
00089       mCatFile.readKey(FitsFile::cAMSTART, pamss);
00090       mCatFile.readKey(FitsFile::cAMEND, pamse);
00091       mCatFile.readKey(FitsFile::cEXP_TIME, pexptime);
00092     }
00093 
00094     // Get observation data from file
00095     for (int hduNo = 2; hduNo <= numHDUs; ++hduNo)
00096     {
00097       int extNo = hduNo - 2;
00098       // Move to extension on the FITS file
00099       mCatFile.movabsHdu(hduNo);
00100 
00101       // Set/Read the airmass values and the exposure time
00102       if (StringOps::ToLower(options.getArchive()) == "vsa")
00103       {
00104         amss[extNo] = pamss;
00105         amse[extNo] = pamse;
00106         exptime[extNo] = pexptime;
00107       }
00108       else
00109       {
00110         mCatFile.readKey(FitsFile::cAMSTART, amss[extNo]);
00111         mCatFile.readKey(FitsFile::cAMEND, amse[extNo]);
00112         mCatFile.readKey(FitsFile::cEXP_TIME, exptime[extNo]);
00113       }
00114 
00115       // Read the zero point
00116       mCatFile.readKey(FitsFile::cMAGZPT, zpmag[extNo]);
00117       mCatFile.readKey(FitsFile::cMAGZRR, zperr[extNo]);
00118 
00119       // Read the extinction
00120       mCatFile.readKey(FitsFile::cEXTINCT, extinct[extNo]);
00121 
00122       // Read in the astrometry values
00123       mCatFile.readKey(FitsFile::cCTYPE1 ,ctype1[extNo]);
00124       mCatFile.readKey(FitsFile::cCTYPE2 ,ctype2[extNo]);
00125       mCatFile.readKey(FitsFile::cCRPIX1, crpix1[extNo]);
00126       mCatFile.readKey(FitsFile::cCRPIX2, crpix2[extNo]);
00127       mCatFile.readKey(FitsFile::cCD11, cd11[extNo]);
00128       mCatFile.readKey(FitsFile::cCD21, cd21[extNo]);
00129       mCatFile.readKey(FitsFile::cCD12, cd12[extNo]);
00130       mCatFile.readKey(FitsFile::cCD22, cd22[extNo]);
00131       if (ctype1[extNo].find("ZPN") != std::string::npos 
00132           && ctype2[extNo].find("ZPN") != std::string::npos) 
00133       {
00134           mCatFile.readKey(FitsFile::cPV21, pv21[extNo]);
00135           mCatFile.readKey(FitsFile::cPV23, pv23[extNo]);
00136       }
00137       else
00138       {
00139           pv21[extNo] = 1.;
00140           pv23[extNo] = 0.;
00141       }
00142     }
00143     for (int rowNo = row1; rowNo < row2; ++rowNo)
00144     {
00145       Numeric value;
00146       data.value(mExtnumColno, rowNo, value);
00147       int extNum = (int)value - 2;
00148       
00149       // Read the "flux"
00150       double flux;
00151       inpData.value(mFluxColno, rowNo, flux);
00152 
00153       // Read the x-position 
00154       double xpos;
00155       inpData.value(mXColno, rowNo, xpos);
00156 
00157       // Read the y-position 
00158       double ypos;
00159       inpData.value(mYColno, rowNo, ypos);
00160 
00161       // Calculate xi and eta
00162       double xi = (cd11[extNum]*pi/180.)*(xpos-crpix1[extNum]) + 
00163                   (cd12[extNum]*pi/180.)*(ypos-crpix2[extNum]);
00164       double eta = (cd21[extNum]*pi/180.)*(xpos-crpix1[extNum]) + 
00165                    (cd22[extNum]*pi/180.)*(ypos-crpix2[extNum]);
00166       double rad = sqrt(xi*xi+eta*eta);
00167 
00168       // Read Aperture correction (written into AperMag1-7)
00169       double aperCor = 0.0;
00170       if (data.isDefined(mMagColno, rowNo))
00171       {
00172         data.value(mMagColno, rowNo, aperCor);
00173       }
00174       // Instrumental magnitude corresponding to this flux
00175       // Need to put in Luptitudes - very rough for now
00176       //if (flux > FLT_EPSILON) {
00177       if (inpData.isDefined(mFluxColno, rowNo) && flux > FLT_EPSILON)
00178       {
00179         double instMag = -flux2mag * log(flux);
00180 
00181         // Correction calculations
00182         // Extinction correction
00183         double extCorr = -extinct[extNum] 
00184                          * ((amss[extNum]+amse[extNum])/2.0 - 1.);
00185 
00186         // Exposure time correction
00187         double texpCorr = flux2mag * log(exptime[extNum]);
00188 
00189         // Distortion correction
00190         double distortCorr;
00191         if (frameType.find(PreSuffixNames::mosaSuffix) != string::npos)
00192         {
00193           // don't apply to mosaics
00194           distortCorr = 0.0;
00195         }
00196         else
00197         {
00198           if (ctype1[extNum].find("ZPN") != std::string::npos && 
00199               ctype2[extNum].find("ZPN") != std::string::npos)
00200           {
00201             distortCorr = 1.0 + 3.0*pv23[extNum]*rad*rad/pv21[extNum];
00202             distortCorr *= 1.0 + pv23[extNum]*rad*rad/pv21[extNum];
00203           } else if (ctype1[extNum].find("TAN") != std::string::npos && 
00204                      ctype2[extNum].find("TAN") != std::string::npos) {
00205             if (rad == 0.)
00206             {
00207               distortCorr = 1.0;
00208             } else {
00209               distortCorr = 1.0/(1.0+tan(rad)*tan(rad));
00210               distortCorr *= atan(rad)/rad;
00211             }
00212           } else {
00213             distortCorr = 1.0;
00214           }
00215             
00216           distortCorr = flux2mag * log(distortCorr);
00217         }
00218         
00219         // Illumination correction
00220         double illumCorr;
00221         if (frameType.find(PreSuffixNames::deepSuffix) != string::npos)
00222         {
00223           // don't apply to deep images
00224           illumCorr = 0.0;
00225         }
00226         else
00227         {
00228           Intervall xiIntervall, etaIntervall;
00229           mIllum->setXiIntervall(xi);
00230           mIllum->setEtaIntervall(eta);
00231           xiIntervall = mIllum->getXiIntervall();
00232           etaIntervall = mIllum->getEtaIntervall();
00233           double f1=0.0, f2=0.0, f3=0.0, f4=0.0, d1=0.0, d2=0.0;
00234           d1 = xiIntervall.val2 - xiIntervall.val1;
00235           d2 = etaIntervall.val2 - etaIntervall.val1;
00236           if (d1 == 0.0 and d2 != 0.0)
00237           {
00238             f1 = mIllum->getMag(xiIntervall.x1, etaIntervall.x1);
00239             f2 = mIllum->getMag(xiIntervall.x1, etaIntervall.x2) - f1;
00240             f2 *= (eta - etaIntervall.val1)/d2;
00241             illumCorr = f1 + f2;
00242           }
00243           else if (d2 == 0.0 and d1 != 0.0)
00244           {
00245             f1 = mIllum->getMag(xiIntervall.x1, etaIntervall.x1);
00246             f2 = mIllum->getMag(xiIntervall.x2, etaIntervall.x1) - f1;
00247             f2 *= (xi - xiIntervall.val1)/d1;
00248             illumCorr = f1 + f2;
00249           }
00250           else if (d1 == 0.0 and d2 == 0.0)
00251           {
00252             illumCorr = mIllum->getMag(xiIntervall.x1, etaIntervall.x1);
00253           }
00254           else
00255           {
00256             f1 = mIllum->getMag(xiIntervall.x1, etaIntervall.x1);
00257             f1 *= (xiIntervall.val2 - xi);
00258             f2 = mIllum->getMag(xiIntervall.x2, etaIntervall.x1);
00259             f2 *= (xi - xiIntervall.val1);
00260             illumCorr = (f1 + f2) * (etaIntervall.val2 - eta);
00261             f3 = mIllum->getMag(xiIntervall.x1, etaIntervall.x2);
00262             f3 *= (xiIntervall.val2 - xi);
00263             f4 = mIllum->getMag(xiIntervall.x2, etaIntervall.x2);
00264             f4 *= (xi - xiIntervall.val1);
00265             illumCorr += (f3 + f4) * (eta - etaIntervall.val1);
00266             illumCorr /= (d1 * d2);
00267           }
00268         }
00269         
00270         // Calibrated magnitude
00271         double calMag = instMag + zpmag[extNum] 
00272           + extCorr + texpCorr - aperCor 
00273           + distortCorr + illumCorr;
00274         
00275         data.assign(mMagColno, rowNo, calMag);
00276       }
00277       else
00278       {
00279         data.assign(mMagColno, rowNo, DefaultValues::realdef);
00280       }
00281     }
00282   }
00283 
00284   // No simple data model available
00285   void doit(TableData<DataType>& data, int row1, int row2) { }
00286 
00287 private:
00288   FitsFile mCatFile;
00289   IllumTable* mIllum;
00290   int mExtnumColno;
00291   int mFluxColno;
00292   int mXColno;
00293   int mYColno;
00294   int mMagColno;
00295   int mExtnumTblno, mFluxTblno, mXTblno, mYTblno;
00296 };
00297 //-----------------------------------------------------------------------------
00300 template<typename DataType>
00301 class MagnitudeError : public DataMethod<DataType>
00302 {
00303 public:
00305   MagnitudeError() { DataMethod<DataType>::methodName = "MagnitudeError"; }
00306 
00308   virtual ~MagnitudeError() { }
00309 
00311   void setSource(IntMap& tbls, IntMap& cols, StringMap units)
00312   {
00313     mFluxColno = cols[0];
00314     mFluxErrColno = cols[1];
00315     mFluxTblno = tbls[0];
00316     mFluxErrTblno = tbls[1];
00317     if (mFluxColno < 0 || mFluxErrColno < 0)
00318     {
00319       throw DataMethodException("MagnitudeError method : cannot source");
00320     }
00321   }
00322 
00324   void setTarget(IntMap& cols, StringMap units) { mMagErrColno = cols[0]; }
00325 
00327   void doit2(TableData<DataType>& inpData, TableData<DataType>& data, int row1, int row2)
00328   {
00329     double mag2flux = M_LN10/2.5;
00330     for (int rowNo = row1; rowNo < row2; ++rowNo)
00331     {
00332       // Read the "flux" and error in the flux
00333       double flux, eflux;
00334       inpData.value(mFluxColno, rowNo, flux);
00335       inpData.value(mFluxErrColno, rowNo, eflux);
00336 
00337       // Error in magnitude
00338       if (inpData.isDefined(mFluxColno, rowNo) &&
00339           inpData.isDefined(mFluxErrColno, rowNo) &&
00340           flux > FLT_EPSILON)
00341       {
00342         double merr = eflux / (mag2flux * flux);
00343         data.assign(mMagErrColno, rowNo, merr);
00344       }
00345       else
00346       {
00347         data.assign(mMagErrColno, rowNo, DefaultValues::realdef);
00348       }
00349     }
00350   }
00351 
00352   // No simple data model available
00353    void doit(TableData<DataType>& data, int row1, int row2) { }
00354 
00355 private:
00356   int mFluxColno;
00357   int mFluxErrColno;
00358   int mMagErrColno;
00359   int mFluxTblno, mFluxErrTblno;
00360 };
00361 //-----------------------------------------------------------------------------
00364 template<typename DataType>
00365 class ColourIndex : public DataMethod<DataType>
00366 {
00367 public:
00369   ColourIndex() { DataMethod<DataType>::methodName = "ColourIndex"; }
00370 
00372   virtual ~ColourIndex() { }
00373 
00376   void setSource(IntMap& tbls, IntMap& cols, StringMap units)
00377   {
00378     mag1Colno = cols[0];
00379     mag2Colno = cols[1];
00380     if (mag1Colno < 0 || mag2Colno < 0)
00381     {
00382       throw DataMethodException("ColourIndex method : cannot source; "
00383                                 + StringOps::NumToString(mag1Colno) + ";"
00384                                 + StringOps::NumToString(mag2Colno));
00385     }
00386   }
00387 
00389   void setTarget(IntMap& cols, StringMap units) { ciColno = cols[0]; }
00390 
00391   void doit(TableData<DataType>& data, int row1, int row2)
00392   {
00393     for (int rowNo = row1; rowNo < row2; ++rowNo)
00394     {
00395       // Do the calculation if both magnitudes are defined - otherwise
00396       // leave them and a default will be assigned later.
00397       if (data.isDefined(mag1Colno,rowNo) && data.isDefined(mag2Colno,rowNo))
00398       {
00399         // Read the two magnitudes
00400         Numeric mag1, mag2;
00401         data.value(mag1Colno, rowNo, mag1);
00402         data.value(mag2Colno, rowNo, mag2);
00403 
00404         // Assign the colour index only if both magnitudes are not default:
00405         if (mag1 > DefaultValues::realdef/2.0 &&
00406             mag2 > DefaultValues::realdef/2.0)
00407         {
00408           data.assign(ciColno, rowNo, mag1 - mag2);
00409         }
00410       }
00411     }
00412   }
00413 
00414   // No data from raw data included
00415   void doit2(TableData<DataType>& inpData, TableData<DataType>& data, int row1, int row2) { }
00416 
00417 private:
00418   int mag1Colno;
00419   int mag2Colno;
00420   int ciColno;
00421 };
00422 //-----------------------------------------------------------------------------
00425 template<typename DataType>
00426 class ColourIndexError : public DataMethod<DataType>
00427 {
00428 public:
00430   ColourIndexError() { DataMethod<DataType>::methodName = "ColourIndexError"; }
00431 
00433   virtual ~ColourIndexError() { }
00434 
00437   void setSource(IntMap& tbls, IntMap& cols, StringMap units)
00438   {
00439     err1Colno = cols[0];
00440     err2Colno = cols[1];
00441     if (err1Colno < 0 || err2Colno < 0)
00442     {
00443       throw DataMethodException("ColourIndexError method : cannot source");
00444     }
00445   }
00446 
00448   void setTarget(IntMap& cols, StringMap units) { ciErrColno = cols[0]; }
00449 
00451   void doit(TableData<DataType>& data, int row1, int row2)
00452   {
00453     for (int rowNo = row1; rowNo < row2; ++rowNo)
00454     {
00455       // Do the calculation if both magnitude errors are defined
00456       if (data.isDefined(err1Colno,rowNo) &&
00457           data.isDefined(err2Colno,rowNo))
00458       {
00459         // Read the two magnitudes
00460         Numeric err1, err2;
00461         data.value(err1Colno, rowNo, err1);
00462         data.value(err2Colno, rowNo, err2);
00463 
00464         // Assign the colour index error which is simply the quadratic
00465         // sum of the errors in the individual errors if both errors are not
00466         // default:
00467         if (err1 > DefaultValues::realdef/2.0 &&
00468             err2 > DefaultValues::realdef/2.0)
00469         {
00470           Numeric err = sqrt(err1 * err1 + err2 * err2);
00471           data.assign(ciErrColno, rowNo, err);
00472         }
00473       }
00474     }
00475   }
00476 
00477   // No data from raw data included
00478   void doit2(TableData<DataType>& inpData, TableData<DataType>& data, int row1, int row2) { }
00479 
00480 private:
00481   int err1Colno;
00482   int err2Colno;
00483   int ciErrColno;
00484 };
00485 //-----------------------------------------------------------------------------
00486 #endif
00487 //-----------------------------------------------------------------------------
00488 // Change log:
00489 //
00490 // 27-May-2004,  IAB: Original version.
00491 //  6-Sep-2004,  NCH: Added this header text; fixed bugs in colour index
00492 //                    member functions where only first colour checked for
00493 //                    defined values.
00494 // 10-Mar-2005, ETWS: Included default value from Constants for magnitudes
00495 // 28-Apr-2005, ETWS: Included extinction, zeropoint, and exposure time
00496 //                    corrections
00497 //  6-Jun-2005,  NCH: Made ColourIndex and ColourIndexError default if
00498 //                    either of the source quantities used in the calculation
00499 //                    are defaults.
00500 // 22-Aug-2005, ETWS: Included aperture correction in Magnitude calculation.
00501 //  7-Dec-2005, ETWS: Included move of Fits Keys into class FitsFile
00502 //  7-Apr-2008, ETWS: Upgraded to use new detection table layout.
00503 // 17-Apr-2008, ETWS: Included illumination correction.
00504 //  5-Jun-2008, ETWS: Fixed illumination correction for cases where linear 
00505 //                    interpolation is needed.
Generated on Mon Oct 4 10:38:34 2010 for WfcamSrc by  doxygen 1.6.3