00001
00002
00003
00004
00005
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
00107 std::string frameType = DataOps::getFrameType(mCatFile.getFileName());
00108
00109 if (StringOps::ToLower(options.getArchive()) == "vsa")
00110 {
00111 mCatFile.movabsHdu(1);
00112
00113 mCatFile.readKey(FitsFile::cAMSTART, pamss);
00114 mCatFile.readKey(FitsFile::cAMEND, pamse);
00115 mCatFile.readKey(FitsFile::cEXP_TIME, pexptime);
00116 }
00117
00118
00119 for (int hduNo = 2; hduNo <= numHDUs; ++hduNo)
00120 {
00121 int extNo = hduNo - 2;
00122
00123 mCatFile.movabsHdu(hduNo);
00124
00125
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
00140 mCatFile.readKey(FitsFile::cMAGZPT, zpmag[extNo]);
00141 mCatFile.readKey(FitsFile::cMAGZRR, zperr[extNo]);
00142
00143
00144 mCatFile.readKey(FitsFile::cEXTINCT, extinct[extNo]);
00145
00146
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
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
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
00189 double pHeight;
00190 inpData.value(mPHeightColno, rowNo, pHeight);
00191
00192
00193 double sky;
00194 inpData.value(mSkyColno, rowNo, sky);
00195
00196
00197 double xpos;
00198 inpData.value(mXColno, rowNo, xpos);
00199
00200
00201 double ypos;
00202 inpData.value(mYColno, rowNo, ypos);
00203
00204
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
00212 double aperCor = 0.0;
00213 if (data.isDefined(mMagColno, rowNo))
00214 {
00215 data.value(mMagColno, rowNo, aperCor);
00216 }
00217
00218
00219
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
00242
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
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
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
00272 diffMag = 2.5 * log10(max(1.0, (isoFlux - aperFlux3) / (1.0 - p1)));
00273
00274
00275 if (sky + pHeight > 0.9 * saturate[extNum])
00276 {
00277
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
00299 if (frameType.find(PreSuffixNames::mosaSuffix) != string::npos
00300 || mCatFile.isTile())
00301 {
00302
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
00331 if (frameType.find(PreSuffixNames::deepSuffix) != string::npos)
00332 {
00333
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
00384 double extCorr = -extinct[extNum]
00385 * ((amss[extNum]+amse[extNum])/2.0 - 1.);
00386
00387
00388 double texpCorr = flux2mag * log(exptime[extNum]);
00389
00390 deltaMag = extCorr + texpCorr + distortCorr + illumCorr;
00391 }
00392
00393
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
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
00447
00448