VSAMagnitudes.hxx

Go to the documentation of this file.
00001 
00002 
00003 
00004 //-----------------------------------------------------------------------------
00005 // SVN: $Id: VSAMagnitudes.hxx 7300 2010-08-05 13:29:29Z EckhardSutorius $
00006 //-----------------------------------------------------------------------------
00007 #ifndef VSAMAGNITUDES_H
00008 #define VSAMAGNITUDES_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 CalibratedVSAMagnitude : public DataMethod<DataType>
00026 {
00027 public:
00032   CalibratedVSAMagnitude(const FitsFile& aCatFile, IllumTable& anIllum,
00033                          const std::string &magName)
00034     : mCatFile(aCatFile)
00035     { 
00036       DataMethod<DataType>::methodName = "CalibratedVSAMagnitude"; 
00037       mIllum = &anIllum;
00038       mMagName = magName;
00039     }
00040 
00042   virtual ~CalibratedVSAMagnitude() { }
00043 
00045   void setSource(IntMap& tbls, IntMap& cols, StringMap units)
00046   {
00047     mExtnumColno = cols[0];
00048     mFluxColno = cols[1];
00049     mXColno = cols[2];
00050     mYColno = cols[3];
00051     mIsoFluxColno = cols[4];
00052     mAperFlux3Colno = cols[5];
00053     mAperFlux4Colno = cols[6];
00054     mAperFlux5Colno = cols[7];
00055     mSkyColno = cols[8];
00056     mPHeightColno = cols[9];
00057 
00058     mExtnumTblno = tbls[0];
00059     mFluxTblno = tbls[1];
00060     mXTblno = tbls[2];
00061     mYTblno = tbls[3];
00062     mIsoFluxTblno = tbls[4];
00063     mAperFlux3Tblno = tbls[5];
00064     mAperFlux4Tblno = tbls[6];
00065     mAperFlux5Tblno = tbls[7];
00066     mSkyTblno = tbls[8];
00067     mPHeightTblno = tbls[9];
00068 
00069     if (mFluxColno < 0)
00070     {
00071       throw DataMethodException("CalibratedVSAMagnitude method: cannot source");
00072     }
00073   }
00074 
00076   void setTarget(IntMap& cols, StringMap units) { 
00077     mMagColno = cols[0]; 
00078     mDeltaMagColno = cols[1];
00079     mIllumCorrColno = cols[2];
00080     mDistortCorrColno = cols[3];
00081     mSaturatCorrColno = cols[4];
00082     mMagNoCorrColno = cols[5];
00083   }
00084 
00089   void doit2(TableData<DataType>& inpData, TableData<DataType>& data, int row1, int row2)
00090   {
00091     Options options;
00092     double pi = M_PI;
00093     double flux2mag = 2.5/M_LN10;
00094 
00095     int numHDUs = mCatFile.getNumHdus();
00096     double pamss, pamse, amss[numHDUs], amse[numHDUs];
00097     double zpmag[numHDUs], zperr[numHDUs];
00098     double extinct[numHDUs];
00099     double pexptime, exptime[numHDUs];
00100     std::string ctype1[numHDUs],ctype2[numHDUs];
00101     double crpix1[numHDUs],crpix2[numHDUs];
00102     double cd11[numHDUs],cd12[numHDUs],cd21[numHDUs],cd22[numHDUs];
00103     double pv21[numHDUs],pv23[numHDUs];
00104     double apcor3[numHDUs],apcor4[numHDUs],apcor5[numHDUs],saturate[numHDUs];
00105 
00106     // get frame type
00107     std::string frameType = DataOps::getFrameType(mCatFile.getFileName());
00108 
00109     if (StringOps::ToLower(options.getArchive()) == "vsa")
00110     {
00111       mCatFile.movabsHdu(1);
00112       // Read the airmass values from primary header for VSA
00113       mCatFile.readKey(FitsFile::cAMSTART, pamss);
00114       mCatFile.readKey(FitsFile::cAMEND, pamse);
00115       mCatFile.readKey(FitsFile::cEXP_TIME, pexptime);
00116     }
00117 
00118     // Get observation data from file
00119     for (int hduNo = 2; hduNo <= numHDUs; ++hduNo)
00120     {
00121       int extNo = hduNo - 2;
00122       // Move to extension on the FITS file
00123       mCatFile.movabsHdu(hduNo);
00124 
00125       // Set/Read the airmass values and the exposure time
00126       if (StringOps::ToLower(options.getArchive()) == "vsa")
00127       {
00128         amss[extNo] = pamss;
00129         amse[extNo] = pamse;
00130         exptime[extNo] = pexptime;
00131       }
00132       else
00133       {
00134         mCatFile.readKey(FitsFile::cAMSTART, amss[extNo]);
00135         mCatFile.readKey(FitsFile::cAMEND, amse[extNo]);
00136         mCatFile.readKey(FitsFile::cEXP_TIME, exptime[extNo]);
00137       }
00138 
00139       // Read the zero point
00140       mCatFile.readKey(FitsFile::cMAGZPT, zpmag[extNo]);
00141       mCatFile.readKey(FitsFile::cMAGZRR, zperr[extNo]);
00142 
00143       // Read the extinction
00144       mCatFile.readKey(FitsFile::cEXTINCT, extinct[extNo]);
00145 
00146       // Read the aperture corrections and saturate
00147       mCatFile.readKey("APCOR3", apcor3[extNo]);
00148       mCatFile.readKey("APCOR4", apcor4[extNo]);
00149       mCatFile.readKey("APCOR5", apcor5[extNo]);
00150       mCatFile.readKey("SATURATE", saturate[extNo]);
00151 
00152       // Read in the astrometry values
00153       mCatFile.readKey(FitsFile::cCTYPE1 ,ctype1[extNo]);
00154       mCatFile.readKey(FitsFile::cCTYPE2 ,ctype2[extNo]);
00155       mCatFile.readKey(FitsFile::cCRPIX1, crpix1[extNo]);
00156       mCatFile.readKey(FitsFile::cCRPIX2, crpix2[extNo]);
00157       mCatFile.readKey(FitsFile::cCD11, cd11[extNo]);
00158       mCatFile.readKey(FitsFile::cCD21, cd21[extNo]);
00159       mCatFile.readKey(FitsFile::cCD12, cd12[extNo]);
00160       mCatFile.readKey(FitsFile::cCD22, cd22[extNo]);
00161       if (ctype1[extNo].find("ZPN") != std::string::npos 
00162           && ctype2[extNo].find("ZPN") != std::string::npos) 
00163       {
00164         mCatFile.readKey(FitsFile::cPV21, pv21[extNo]);
00165         mCatFile.readKey(FitsFile::cPV23, pv23[extNo]);
00166       }
00167       else
00168       {
00169         pv21[extNo] = 1.;
00170         pv23[extNo] = 0.;
00171       }
00172     }
00173 
00174     for (int rowNo = row1; rowNo < row2; ++rowNo)
00175     {
00176       Numeric value;
00177       data.value(mExtnumColno, rowNo, value);
00178       int extNum = (int)value - 2;
00179 
00180       // Read the "fluxes"
00181       double flux, isoFlux, aperFlux3, aperFlux4, aperFlux5;
00182       inpData.value(mFluxColno, rowNo, flux);
00183       inpData.value(mIsoFluxColno, rowNo, isoFlux);
00184       inpData.value(mAperFlux3Colno, rowNo, aperFlux3);
00185       inpData.value(mAperFlux4Colno, rowNo, aperFlux4);
00186       inpData.value(mAperFlux5Colno, rowNo, aperFlux5);
00187 
00188       // Read the peak height
00189       double pHeight;
00190       inpData.value(mPHeightColno, rowNo, pHeight);
00191 
00192       // Read the sky
00193       double sky;
00194       inpData.value(mSkyColno, rowNo, sky);
00195 
00196       // Read the x-position 
00197       double xpos;
00198       inpData.value(mXColno, rowNo, xpos);
00199 
00200       // Read the y-position 
00201       double ypos;
00202       inpData.value(mYColno, rowNo, ypos);
00203 
00204       // Calculate xi and eta
00205       double xi = (cd11[extNum]*pi/180.)*(xpos-crpix1[extNum]) + 
00206                   (cd12[extNum]*pi/180.)*(ypos-crpix2[extNum]);
00207       double eta = (cd21[extNum]*pi/180.)*(xpos-crpix1[extNum]) + 
00208                    (cd22[extNum]*pi/180.)*(ypos-crpix2[extNum]);
00209       double rad = sqrt(xi*xi+eta*eta);
00210 
00211       // Read Aperture correction (written into AperMag1-7)
00212       double aperCor = 0.0;
00213       if (data.isDefined(mMagColno, rowNo))
00214       {
00215         data.value(mMagColno, rowNo, aperCor);
00216       }
00217 
00218       // Correction calculations
00219       // calculate only if deltaMag is not set yet
00220       double deltaMag = DefaultValues::realdef;
00221       double illumCorr = DefaultValues::realdef;
00222       double distortCorr = DefaultValues::realdef;
00223       double saturatCorr = DefaultValues::realdef;   
00224       if (data.isDefined(mDeltaMagColno, rowNo))
00225       {
00226         data.value(mDeltaMagColno, rowNo, deltaMag);
00227       }
00228       if (data.isDefined(mIllumCorrColno, rowNo))
00229       {
00230         data.value(mIllumCorrColno, rowNo, illumCorr);
00231       }
00232       if (data.isDefined(mDistortCorrColno, rowNo))
00233       {
00234         data.value(mDistortCorrColno, rowNo, distortCorr);
00235       }
00236       if (data.isDefined(mSaturatCorrColno, rowNo))
00237       {
00238         data.value(mSaturatCorrColno, rowNo, saturatCorr);
00239       }
00240 
00241       // Instrumental magnitude corresponding to this flux
00242       // Need to put in Luptitudes - very rough for now
00243       if (inpData.isDefined(mFluxColno, rowNo) && flux > FLT_EPSILON)
00244       {
00245         double instMag = -flux2mag * log(flux);
00246 
00247         if (mMagName == "aperMag3" && saturatCorr < DefaultValues::realdef/2.0)
00248         {
00249           // Saturation correction          
00250           double p1, p2, p3, dat2, dat3;
00251           double top = 0.0;
00252           double bot = 0.0;
00253           double profMag = 0.0;
00254           double diffMag = 0.0;
00255           double estMag = 0.0;
00256           double corrMag = 0.0;
00257           double apcor = 0.0;
00258 
00259           // profile estimate
00260           p1 = pow(10, -0.4*max(0.0, apcor3[extNum]));
00261           p2 = pow(10, -0.4*max(0.0, apcor4[extNum])) 
00262             - pow(10, -0.4*max(0.0, apcor3[extNum]));
00263           p3 = pow(10, -0.4*max(0.0, apcor5[extNum])) 
00264             - pow(10, -0.4*max(0.0, apcor4[extNum]));
00265           dat2 = max(0.0, aperFlux4 - aperFlux3);
00266           dat3 = max(0.0, aperFlux5 - aperFlux4);
00267           top = dat2 * p2 + dat3 * p3;
00268           bot = p2 * p2 + p3 * p3;
00269           profMag = 2.5 * log10(max(1.0, top / bot));
00270 
00271           // alternative based on isoFlux
00272           diffMag = 2.5 * log10(max(1.0, (isoFlux - aperFlux3) / (1.0 - p1)));
00273 
00274           // calculate saturation correction
00275           if (sky + pHeight > 0.9 * saturate[extNum])
00276           {
00277             // usual estimate
00278             apcor = pow(10, 0.4 * apcor3[extNum]);
00279             estMag = 2.5 * log10(max(1.0, aperFlux3 * apcor));
00280             if(abs(profMag - diffMag) < 0.1)
00281             {
00282               corrMag = -max(estMag - 0.025, max(profMag, diffMag));
00283             }
00284             else
00285             {
00286               corrMag = -max(estMag - 0.025, min(profMag, diffMag));
00287             }
00288             saturatCorr = min(0.0, corrMag + estMag);
00289           }
00290           else
00291           {
00292             saturatCorr = 0.0;
00293           }
00294         }
00295 
00296         if (distortCorr < DefaultValues::realdef/2.0)
00297         {
00298           // Distortion correction
00299           if (frameType.find(PreSuffixNames::mosaSuffix) != string::npos
00300               || mCatFile.isTile())
00301           {
00302             // don't apply to mosaics or tiles
00303             distortCorr = 0.0;
00304           }
00305           else
00306           {
00307             if (ctype1[extNum].find("ZPN") != std::string::npos && 
00308                 ctype2[extNum].find("ZPN") != std::string::npos)
00309             {
00310               distortCorr = 1.0 + 3.0*pv23[extNum]*rad*rad/pv21[extNum];
00311               distortCorr *= 1.0 + pv23[extNum]*rad*rad/pv21[extNum];
00312             } else if (ctype1[extNum].find("TAN") != std::string::npos && 
00313                        ctype2[extNum].find("TAN") != std::string::npos) {
00314               if (rad == 0.)
00315               {
00316                 distortCorr = 1.0;
00317               } else {
00318                 distortCorr = 1.0/(1.0+tan(rad)*tan(rad));
00319                 distortCorr *= atan(rad)/rad;
00320               }
00321             } else {
00322               distortCorr = 1.0;
00323             }
00324             distortCorr = flux2mag * log(distortCorr);
00325           }
00326         }
00327 
00328         if (illumCorr < DefaultValues::realdef/2.0)
00329         {
00330           // Illumination correction
00331           if (frameType.find(PreSuffixNames::deepSuffix) != string::npos)
00332           {
00333             // don't apply to deep images
00334             illumCorr = 0.0;
00335           }
00336           else
00337           {
00338             Intervall xiIntervall, etaIntervall;
00339             mIllum->setXiIntervall(xi);
00340             mIllum->setEtaIntervall(eta);
00341             xiIntervall = mIllum->getXiIntervall();
00342             etaIntervall = mIllum->getEtaIntervall();
00343             double f1=0.0, f2=0.0, f3=0.0, f4=0.0, d1=0.0, d2=0.0;
00344             d1 = xiIntervall.val2 - xiIntervall.val1;
00345             d2 = etaIntervall.val2 - etaIntervall.val1;
00346             if (d1 == 0.0 and d2 != 0.0)
00347             {
00348               f1 = mIllum->getMag(xiIntervall.x1, etaIntervall.x1);
00349               f2 = mIllum->getMag(xiIntervall.x1, etaIntervall.x2) - f1;
00350               f2 *= (eta - etaIntervall.val1)/d2;
00351               illumCorr = f1 + f2;
00352             }
00353             else if (d2 == 0.0 and d1 != 0.0)
00354             {
00355               f1 = mIllum->getMag(xiIntervall.x1, etaIntervall.x1);
00356               f2 = mIllum->getMag(xiIntervall.x2, etaIntervall.x1) - f1;
00357               f2 *= (xi - xiIntervall.val1)/d1;
00358               illumCorr = f1 + f2;
00359             }
00360             else if (d1 == 0.0 and d2 == 0.0)
00361             {
00362               illumCorr = mIllum->getMag(xiIntervall.x1, etaIntervall.x1);
00363             }
00364             else
00365             {
00366               f1 = mIllum->getMag(xiIntervall.x1, etaIntervall.x1);
00367               f1 *= (xiIntervall.val2 - xi);
00368               f2 = mIllum->getMag(xiIntervall.x2, etaIntervall.x1);
00369               f2 *= (xi - xiIntervall.val1);
00370               illumCorr = (f1 + f2) * (etaIntervall.val2 - eta);
00371               f3 = mIllum->getMag(xiIntervall.x1, etaIntervall.x2);
00372               f3 *= (xiIntervall.val2 - xi);
00373               f4 = mIllum->getMag(xiIntervall.x2, etaIntervall.x2);
00374               f4 *= (xi - xiIntervall.val1);
00375               illumCorr += (f3 + f4) * (eta - etaIntervall.val1);
00376               illumCorr /= (d1 * d2);
00377             }
00378           }
00379         }
00380 
00381         if (deltaMag < DefaultValues::realdef/2.0)
00382         {
00383           // Extinction correction
00384           double extCorr = -extinct[extNum] 
00385             * ((amss[extNum]+amse[extNum])/2.0 - 1.);
00386 
00387           // Exposure time correction
00388           double texpCorr = flux2mag * log(exptime[extNum]);
00389           
00390           deltaMag = extCorr + texpCorr + distortCorr + illumCorr;
00391         }
00392 
00393         // Calibrated magnitude
00394         double calMag = instMag + zpmag[extNum] - aperCor + deltaMag;
00395         if (mMagName == "aperMag3")
00396         {
00397           calMag += saturatCorr;
00398         }
00399 
00400         data.assign(mMagColno, rowNo, calMag);
00401         data.assign(mDeltaMagColno, rowNo, deltaMag);
00402         data.assign(mIllumCorrColno, rowNo, illumCorr);
00403         data.assign(mDistortCorrColno, rowNo, distortCorr);
00404         data.assign(mSaturatCorrColno, rowNo, saturatCorr);
00405         if (mMagNoCorrColno != mMagColno)
00406         {
00407           data.assign(mMagNoCorrColno, rowNo, calMag += aperCor);
00408         }
00409       }
00410       else
00411       {
00412         data.assign(mMagColno, rowNo, DefaultValues::realdef);
00413         data.assign(mMagNoCorrColno, rowNo, DefaultValues::realdef);
00414         if (deltaMag < DefaultValues::realdef/2.0)
00415         {
00416           data.assign(mDeltaMagColno, rowNo, 0.0);
00417           data.assign(mIllumCorrColno, rowNo, 0.0);
00418           data.assign(mDistortCorrColno, rowNo, 0.0);
00419           data.assign(mSaturatCorrColno, rowNo, 0.0);
00420         }
00421       }
00422     }
00423   }
00424 
00425   // No simple data model available
00426   void doit(TableData<DataType>& data, int row1, int row2) { }
00427 
00428 private:
00429   FitsFile mCatFile;
00430   IllumTable* mIllum;
00431   std::string mMagName;
00432   int mExtnumColno;
00433   int mFluxColno;
00434   int mXColno, mYColno;
00435   int mIsoFluxColno, mAperFlux3Colno, mAperFlux4Colno, mAperFlux5Colno;
00436   int mSkyColno, mPHeightColno;
00437   int mMagColno, mMagNoCorrColno;
00438   int mDeltaMagColno, mIllumCorrColno, mDistortCorrColno, mSaturatCorrColno;
00439   int mExtnumTblno, mFluxTblno, mXTblno, mYTblno;
00440   int mIsoFluxTblno, mAperFlux3Tblno, mAperFlux4Tblno, mAperFlux5Tblno;
00441   int mSkyTblno, mPHeightTblno;
00442 };
00443 //-----------------------------------------------------------------------------
00444 #endif
00445 //-----------------------------------------------------------------------------
00446 // Change log:
00447 //
00448 //  9-Dec-2009, ETWS: Original version based on Magnitudes.hxx.