00001
00002
00003
00004
00005
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
00083 std::string frameType = DataOps::getFrameType(mCatFile.getFileName());
00084
00085 if (StringOps::ToLower(options.getArchive()) == "vsa")
00086 {
00087 mCatFile.movabsHdu(1);
00088
00089 mCatFile.readKey(FitsFile::cAMSTART, pamss);
00090 mCatFile.readKey(FitsFile::cAMEND, pamse);
00091 mCatFile.readKey(FitsFile::cEXP_TIME, pexptime);
00092 }
00093
00094
00095 for (int hduNo = 2; hduNo <= numHDUs; ++hduNo)
00096 {
00097 int extNo = hduNo - 2;
00098
00099 mCatFile.movabsHdu(hduNo);
00100
00101
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
00116 mCatFile.readKey(FitsFile::cMAGZPT, zpmag[extNo]);
00117 mCatFile.readKey(FitsFile::cMAGZRR, zperr[extNo]);
00118
00119
00120 mCatFile.readKey(FitsFile::cEXTINCT, extinct[extNo]);
00121
00122
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
00150 double flux;
00151 inpData.value(mFluxColno, rowNo, flux);
00152
00153
00154 double xpos;
00155 inpData.value(mXColno, rowNo, xpos);
00156
00157
00158 double ypos;
00159 inpData.value(mYColno, rowNo, ypos);
00160
00161
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
00169 double aperCor = 0.0;
00170 if (data.isDefined(mMagColno, rowNo))
00171 {
00172 data.value(mMagColno, rowNo, aperCor);
00173 }
00174
00175
00176
00177 if (inpData.isDefined(mFluxColno, rowNo) && flux > FLT_EPSILON)
00178 {
00179 double instMag = -flux2mag * log(flux);
00180
00181
00182
00183 double extCorr = -extinct[extNum]
00184 * ((amss[extNum]+amse[extNum])/2.0 - 1.);
00185
00186
00187 double texpCorr = flux2mag * log(exptime[extNum]);
00188
00189
00190 double distortCorr;
00191 if (frameType.find(PreSuffixNames::mosaSuffix) != string::npos)
00192 {
00193
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
00220 double illumCorr;
00221 if (frameType.find(PreSuffixNames::deepSuffix) != string::npos)
00222 {
00223
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
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
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
00333 double flux, eflux;
00334 inpData.value(mFluxColno, rowNo, flux);
00335 inpData.value(mFluxErrColno, rowNo, eflux);
00336
00337
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
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
00396
00397 if (data.isDefined(mag1Colno,rowNo) && data.isDefined(mag2Colno,rowNo))
00398 {
00399
00400 Numeric mag1, mag2;
00401 data.value(mag1Colno, rowNo, mag1);
00402 data.value(mag2Colno, rowNo, mag2);
00403
00404
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
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
00456 if (data.isDefined(err1Colno,rowNo) &&
00457 data.isDefined(err2Colno,rowNo))
00458 {
00459
00460 Numeric err1, err2;
00461 data.value(err1Colno, rowNo, err1);
00462 data.value(err2Colno, rowNo, err2);
00463
00464
00465
00466
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
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
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505