Magnitudes.hxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// SVN: $Id: Magnitudes.hxx 5292 2008-11-21 16:38:11Z RossCollins $
//-----------------------------------------------------------------------------
#ifndef MAGNITUDES_H
#define MAGNITUDES_H

#include <float.h>

#include "Constants.h"
#include "DataMethod.hxx"
#include "FitsFile.h"
#include "StringOps.h"
#include "IllumTable.hxx"
//-----------------------------------------------------------------------------
template<typename DataType>
class CalibratedMagnitude : public DataMethod<DataType>
{
public:
  CalibratedMagnitude(const FitsFile& aCatFile, IllumTable& anIllum)
    : mCatFile(aCatFile)
    { 
      DataMethod<DataType>::methodName = "CalibratedMagnitude"; 
      mIllum = &anIllum;
    }

  virtual ~CalibratedMagnitude() { }

  void setSource(IntMap& tbls, IntMap& cols, StringMap units)
  {
    mExtnumColno = cols[0];
    mFluxColno = cols[1];
    mXColno = cols[2];
    mYColno = cols[3];
    mExtnumTblno = tbls[0];
    mFluxTblno = tbls[1];
    mXTblno = tbls[2];
    mYTblno = tbls[3];
    if (mFluxColno < 0)
    {
      throw DataMethodException("CalibratedMagnitude method: cannot source");
    }
  }

  void setTarget(IntMap& cols, StringMap units) { mMagColno = cols[0]; }

  void doit2(TableData<DataType>& inpData, TableData<DataType>& data, int row1, int row2)
  {
    double pi = M_PI;
    double flux2mag = 2.5/M_LN10;
    double amss[4], amse[4];
    double zpmag[4], zperr[4];
    double extinct[4];
    double exptime[4];
    std::string ctype1[4],ctype2[4];
    double crpix1[4],crpix2[4],cd11[4],cd12[4],cd21[4],cd22[4];
    double pv21[4],pv23[4];

    // get frame type
    std::string frameType = DataOps::getFrameType(mCatFile.getFileName());
    
    // Get observation data from file
    for (int hduNo = 2; hduNo <= mCatFile.getNumHdus(); ++hduNo)
    {
      int extNo = hduNo - 2;
      // Move to extension on the FITS file
      mCatFile.movabsHdu(hduNo);

      // Read the airmass values
      mCatFile.readKey(FitsFile::cAMSTART, amss[extNo]);
      mCatFile.readKey(FitsFile::cAMEND, amse[extNo]);

      // Read the zero point
      mCatFile.readKey(FitsFile::cMAGZPT, zpmag[extNo]);
      mCatFile.readKey(FitsFile::cMAGZRR, zperr[extNo]);

      // Read the extinction
      mCatFile.readKey(FitsFile::cEXTINCT, extinct[extNo]);

      // Read the exposure time
      mCatFile.readKey(FitsFile::cEXP_TIME, exptime[extNo]);
      
      // Read in the astrometry values
      mCatFile.readKey(FitsFile::cCTYPE1 ,ctype1[extNo]);
      mCatFile.readKey(FitsFile::cCTYPE2 ,ctype2[extNo]);
      mCatFile.readKey(FitsFile::cCRPIX1, crpix1[extNo]);
      mCatFile.readKey(FitsFile::cCRPIX2, crpix2[extNo]);
      mCatFile.readKey(FitsFile::cCD11, cd11[extNo]);
      mCatFile.readKey(FitsFile::cCD21, cd21[extNo]);
      mCatFile.readKey(FitsFile::cCD12, cd12[extNo]);
      mCatFile.readKey(FitsFile::cCD22, cd22[extNo]);
      if (ctype1[extNo].find("ZPN")==5 && ctype2[extNo].find("ZPN")==5) 
      {
          mCatFile.readKey(FitsFile::cPV21, pv21[extNo]);
          mCatFile.readKey(FitsFile::cPV23, pv23[extNo]);
      }
      else
      {
          pv21[extNo]=1.;
          pv23[extNo]=0.;
      }
    }
    for (int rowNo = row1; rowNo < row2; ++rowNo)
    {
      Numeric value;
      data.value(mExtnumColno, rowNo, value);
      int extNum = (int)value - 2;
      
      // Read the "flux"
      double flux;
      inpData.value(mFluxColno, rowNo, flux);

      // Read the x-position 
      double xpos;
      inpData.value(mXColno, rowNo, xpos);

      // Read the y-position 
      double ypos;
      inpData.value(mYColno, rowNo, ypos);

      // Calculate xi and eta
      double xi = (cd11[extNum]*pi/180.)*(xpos-crpix1[extNum]) + 
                  (cd12[extNum]*pi/180.)*(ypos-crpix2[extNum]);
      double eta = (cd21[extNum]*pi/180.)*(xpos-crpix1[extNum]) + 
                   (cd22[extNum]*pi/180.)*(ypos-crpix2[extNum]);
      double rad = sqrt(xi*xi+eta*eta);

      // Read Aperture correction (written into AperMag1-7)
      double aperCor = 0.0;
      if (data.isDefined(mMagColno, rowNo))
      {
        data.value(mMagColno, rowNo, aperCor);
      }
      // Instrumental magnitude corresponding to this flux
      // Need to put in Luptitudes - very rough for now
      //if (flux > FLT_EPSILON) {
      if (inpData.isDefined(mFluxColno, rowNo) && flux > FLT_EPSILON)
      {
        double instMag = -flux2mag * log(flux);

        // Correction calculations
        // Extinction correction
        double extCorr = -extinct[extNum] 
                         * ((amss[extNum]+amse[extNum])/2.0 - 1.);

        // Exposure time correction
        double texpCorr = flux2mag * log(exptime[extNum]);

        // Distortion correction
        double distortCorr;
        if (frameType.find(PreSuffixNames::mosaSuffix) != string::npos)
        {
          // don't apply to mosaics
          distortCorr = 0.0;
        }
        else
        {
          if (ctype1[extNum].find("ZPN")!=std::string::npos && 
              ctype2[extNum].find("ZPN")!=std::string::npos)
          {
            distortCorr = 1.0 + 3.0*pv23[extNum]*rad*rad/pv21[extNum];
            distortCorr *= 1.0 + pv23[extNum]*rad*rad/pv21[extNum];
          } else if (ctype1[extNum].find("TAN")!=std::string::npos && 
                     ctype2[extNum].find("TAN")!=std::string::npos) {
            if (rad==0.)
            {
              distortCorr = 1.0;
            } else {
              distortCorr = 1.0/(1.0+tan(rad)*tan(rad));
              distortCorr *= atan(rad)/rad;
            }
          } else {
            distortCorr = 1.0;
          }
            
          distortCorr = flux2mag * log(distortCorr);
        }
        
        // Illumination correction
        double illumCorr;
        if (frameType.find(PreSuffixNames::deepSuffix) != string::npos)
        {
          // don't apply to deep images
          illumCorr = 0.0;
        }
        else
        {
          Intervall xiIntervall, etaIntervall;
          mIllum->setXiIntervall(xi);
          mIllum->setEtaIntervall(eta);
          xiIntervall = mIllum->getXiIntervall();
          etaIntervall = mIllum->getEtaIntervall();
          double f1=0.0, f2=0.0, f3=0.0, f4=0.0, d1=0.0, d2=0.0;
          d1 = xiIntervall.val2 - xiIntervall.val1;
          d2 = etaIntervall.val2 - etaIntervall.val1;
          if (d1 == 0.0 and d2 != 0.0)
          {
            f1 = mIllum->getMag(xiIntervall.x1, etaIntervall.x1);
            f2 = mIllum->getMag(xiIntervall.x1, etaIntervall.x2) - f1;
            f2 *= (eta - etaIntervall.val1)/d2;
            illumCorr = f1 + f2;
          }
          else if (d2 == 0.0 and d1 != 0.0)
          {
            f1 = mIllum->getMag(xiIntervall.x1, etaIntervall.x1);
            f2 = mIllum->getMag(xiIntervall.x2, etaIntervall.x1) - f1;
            f2 *= (xi - xiIntervall.val1)/d1;
            illumCorr = f1 + f2;
          }
          else if (d1 == 0.0 and d2 == 0.0)
          {
            illumCorr = mIllum->getMag(xiIntervall.x1, etaIntervall.x1);
          }
          else
          {
            f1 = mIllum->getMag(xiIntervall.x1, etaIntervall.x1);
            f1 *= (xiIntervall.val2 - xi);
            f2 = mIllum->getMag(xiIntervall.x2, etaIntervall.x1);
            f2 *= (xi - xiIntervall.val1);
            illumCorr = (f1 + f2) * (etaIntervall.val2 - eta);
            f3 = mIllum->getMag(xiIntervall.x1, etaIntervall.x2);
            f3 *= (xiIntervall.val2 - xi);
            f4 = mIllum->getMag(xiIntervall.x2, etaIntervall.x2);
            f4 *= (xi - xiIntervall.val1);
            illumCorr += (f3 + f4) * (eta - etaIntervall.val1);
            illumCorr /= (d1 * d2);
          }
        }
        
        // Calibrated magnitude
        double calMag = instMag + zpmag[extNum] 
          + extCorr + texpCorr - aperCor 
          + distortCorr + illumCorr;
        
        data.assign(mMagColno, rowNo, calMag);
      }
      else
      {
        data.assign(mMagColno, rowNo, DefaultValues::realdef);
      }
    }
  }

  // No simple data model available
  void doit(TableData<DataType>& data, int row1, int row2) { }

private:
  FitsFile mCatFile;
  IllumTable* mIllum;
  int mExtnumColno;
  int mFluxColno;
  int mXColno;
  int mYColno;
  int mMagColno;
  int mExtnumTblno, mFluxTblno, mXTblno, mYTblno;
};
//-----------------------------------------------------------------------------
template<typename DataType>
class MagnitudeError : public DataMethod<DataType>
{
public:
  MagnitudeError() { DataMethod<DataType>::methodName = "MagnitudeError"; }

  virtual ~MagnitudeError() { }

  void setSource(IntMap& tbls, IntMap& cols, StringMap units)
  {
    mFluxColno = cols[0];
    mFluxErrColno = cols[1];
    mFluxTblno = tbls[0];
    mFluxErrTblno = tbls[1];
    if (mFluxColno < 0 || mFluxErrColno < 0)
    {
      throw DataMethodException("MagnitudeError method : cannot source");
    }
  }

  void setTarget(IntMap& cols, StringMap units) { mMagErrColno = cols[0]; }

  void doit2(TableData<DataType>& inpData, TableData<DataType>& data, int row1, int row2)
  {
    double mag2flux = M_LN10/2.5;
    for (int rowNo = row1; rowNo < row2; ++rowNo)
    {
      // Read the "flux" and error in the flux
      double flux, eflux;
      inpData.value(mFluxColno, rowNo, flux);
      inpData.value(mFluxErrColno, rowNo, eflux);

      // Error in magnitude
      if (inpData.isDefined(mFluxColno, rowNo) &&
          inpData.isDefined(mFluxErrColno, rowNo) &&
          flux > FLT_EPSILON)
      {
        double merr = eflux / (mag2flux * flux);
        data.assign(mMagErrColno, rowNo, merr);
      }
      else
      {
        data.assign(mMagErrColno, rowNo, DefaultValues::realdef);
      }
    }
  }

  // No simple data model available
   void doit(TableData<DataType>& data, int row1, int row2) { }

private:
  int mFluxColno;
  int mFluxErrColno;
  int mMagErrColno;
  int mFluxTblno, mFluxErrTblno;
};
//-----------------------------------------------------------------------------
template<typename DataType>
class ColourIndex : public DataMethod<DataType>
{
public:
  ColourIndex() { DataMethod<DataType>::methodName = "ColourIndex"; }

  virtual ~ColourIndex() { }

  void setSource(IntMap& tbls, IntMap& cols, StringMap units)
  {
    mag1Colno = cols[0];
    mag2Colno = cols[1];
    if (mag1Colno < 0 || mag2Colno < 0)
    {
      throw DataMethodException("ColourIndex method : cannot source; "
                                + StringOps::NumToString(mag1Colno) + ";"
                                + StringOps::NumToString(mag2Colno));
    }
  }

  void setTarget(IntMap& cols, StringMap units) { ciColno = cols[0]; }

  void doit(TableData<DataType>& data, int row1, int row2)
  {
    for (int rowNo = row1; rowNo < row2; ++rowNo)
    {
      // Do the calculation if both magnitudes are defined - otherwise
      // leave them and a default will be assigned later.
      if (data.isDefined(mag1Colno,rowNo) && data.isDefined(mag2Colno,rowNo))
      {
        // Read the two magnitudes
        Numeric mag1, mag2;
        data.value(mag1Colno, rowNo, mag1);
        data.value(mag2Colno, rowNo, mag2);

        // Assign the colour index only if both magnitudes are not default:
        if (mag1 > DefaultValues::realdef/2.0 &&
            mag2 > DefaultValues::realdef/2.0)
        {
          data.assign(ciColno, rowNo, mag1 - mag2);
        }
      }
    }
  }

  // No data from raw data included
  void doit2(TableData<DataType>& inpData, TableData<DataType>& data, int row1, int row2) { }

private:
  int mag1Colno;
  int mag2Colno;
  int ciColno;
};
//-----------------------------------------------------------------------------
template<typename DataType>
class ColourIndexError : public DataMethod<DataType>
{
public:
  ColourIndexError() { DataMethod<DataType>::methodName = "ColourIndexError"; }

  virtual ~ColourIndexError() { }

  void setSource(IntMap& tbls, IntMap& cols, StringMap units)
  {
    err1Colno = cols[0];
    err2Colno = cols[1];
    if (err1Colno < 0 || err2Colno < 0)
    {
      throw DataMethodException("ColourIndexError method : cannot source");
    }
  }

  void setTarget(IntMap& cols, StringMap units) { ciErrColno = cols[0]; }

  void doit(TableData<DataType>& data, int row1, int row2)
  {
    for (int rowNo = row1; rowNo < row2; ++rowNo)
    {
      // Do the calculation if both magnitude errors are defined
      if (data.isDefined(err1Colno,rowNo) &&
          data.isDefined(err2Colno,rowNo))
      {
        // Read the two magnitudes
        Numeric err1, err2;
        data.value(err1Colno, rowNo, err1);
        data.value(err2Colno, rowNo, err2);

        // Assign the colour index error which is simply the quadratic
        // sum of the errors in the individual errors if both errors are not
        // default:
        if (err1 > DefaultValues::realdef/2.0 &&
            err2 > DefaultValues::realdef/2.0)
        {
          Numeric err = sqrt(err1 * err1 + err2 * err2);
          data.assign(ciErrColno, rowNo, err);
        }
      }
    }
  }

  // No data from raw data included
  void doit2(TableData<DataType>& inpData, TableData<DataType>& data, int row1, int row2) { }

private:
  int err1Colno;
  int err2Colno;
  int ciErrColno;
};
//-----------------------------------------------------------------------------
#endif
//-----------------------------------------------------------------------------
// Change log:
//
// 27-May-2004,  IAB: Original version.
//  6-Sep-2004,  NCH: Added this header text; fixed bugs in colour index
//                    member functions where only first colour checked for
//                    defined values.
// 10-Mar-2005, ETWS: Included default value from Constants for magnitudes
// 28-Apr-2005, ETWS: Included extinction, zeropoint, and exposure time
//                    corrections
//  6-Jun-2005,  NCH: Made ColourIndex and ColourIndexError default if
//                    either of the source quantities used in the calculation
//                    are defaults.
// 22-Aug-2005, ETWS: Included aperture correction in Magnitude calculation.
//  7-Dec-2005, ETWS: Included move of Fits Keys into class FitsFile
//  7-Apr-2008, ETWS: Upgraded to use new detection table layout.
// 17-Apr-2008, ETWS: Included illumination correction.
//  5-Jun-2008, ETWS: Fixed illumination correction for cases where linear 
//                    interpolation is needed.

---