HlRadii.hxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// SVN: $Id: HlRadii.hxx 6558 2010-01-19 18:38:52Z RossCollins $
//-----------------------------------------------------------------------------
#ifndef HLRADII_H
#define HLRADII_H

#include <float.h>

#include "Constants.h"
#include "DataMethod.hxx"
#include "FitsFile.h"
#include "StringOps.h"
#include "Logger.h"
#include "DataOps.h"
#include "Converter.h"

//-----------------------------------------------------------------------------
template<typename DataType>
class HlRadii : public DataMethod<DataType>
{
public:
  HlRadii(const FitsFile& aCatFile, const FitsFile& aFile)
   : mCatFile(aCatFile),mFile(aFile)
   {
     DataMethod<DataType>::methodName = "HlRadii";
   }

  virtual ~HlRadii() { }

  //---------------------------------------------------------------------------
  void setSource(IntMap& tbls, IntMap& cols, StringMap units)
  {
    hlcColNo    = cols[0];
    ellColNo   = cols[1];
    extColNo   = cols[2];
    //xColNo= cols[3];
    //yColNo= cols[4];
    //xpsColNo   = cols[5]; // Both x and y?
    //ypsColNo   = cols[6];


  }
  //---------------------------------------------------------------------------
  void setTarget(IntMap& cols, StringMap units)
  {
    hlGeoColNo    = cols[0];
    hlSMnColNo    = cols[1];
    hlSMjColNo    = cols[2];
    hlCorSMnColNo = cols[3];
    hlCorSMjColNo = cols[4];
  }
  //---------------------------------------------------------------------------
  void doit(TableData<DataType>& data, int row1, int row2)
  {

    Logger log;

    // Read in seeing here

    double seeing[mCatFile.getNumHdus()];
    double crval1;
    double crval2;
    int npixx,npixy;

    double centre[2];
    double sizeFact[2];


    // Set up AST
    astBegin;
    AstFitsChan *fitschan = astFitsChan(NULL, NULL, "%.*s", 0, "NULL");


    Converter conv;
    sizeFact[0]  = conv.getFactor("RADIANS", "Arcseconds");
    sizeFact[1]  = conv.getFactor("RADIANS", "Arcseconds");


    for (int hduNo = 2; hduNo <= mCatFile.getNumHdus(); ++hduNo)
    {
      int extNo = hduNo - 2;

      AstFrameSet *wcsinfo = 0;


      // Move to extension on the FITS file
      mCatFile.movabsHdu(hduNo);
      mFile.movabsHdu(hduNo);
      mCatFile.readKey("SEEING", seeing[extNo]);

      // Put all cards in the header of this extension into the AST
      // fits channel
      int ncards = mFile.getNumCards();
      for (int n = 0; n < ncards; ++n) {
          char card[80];
          mFile.readCard(n, card);
          astPutFits(fitschan, card, 0);
      }

      // Make AST wcsinfo object
      astClear( fitschan, "Card" );

      mFile.readKey(FitsFile::CRVAL1, crval1);
      mFile.readKey(FitsFile::CRVAL2, crval2);
      if (crval1 == 0 && crval2 == 0) {
          log.addError("Could not read WCS info for " + mCatFile.getFileName());
      } else {
          wcsinfo = (AstFrameSet*)astRead(fitschan);
      }
      mCatFile.readKey("NXOUT", npixx);
      mCatFile.readKey("NYOUT", npixy);
      centre[0] = (double)npixx / 2;
      centre[1] = (double)npixy / 2;
      Numeric xPixSize,yPixSize;

      DataOps::calcPixSize(wcsinfo,centre,sizeFact,xPixSize,yPixSize);
      seeing[extNo]*=sqrt(xPixSize*yPixSize);

    }


    for (int rowNo = row1; rowNo < row2; ++rowNo)
    {
      double extNum;
      data.value(extColNo, rowNo, extNum);
      int extNo=(int)(extNum-2);

      mCatFile.movabsHdu((int)extNum);

      double hlCirc;
      double ell;

      data.value(hlcColNo, rowNo, hlCirc);


      // if no distortion correction in mags, none in radii

      // Now Get hlSMj, hlSmn,hlGeo
      data.value(ellColNo, rowNo, ell);

      double hlSMj=DefaultValues::realdef;
      double hlSMn=DefaultValues::realdef;
      double hlGeo=DefaultValues::realdef;
      double hlCorSMj=DefaultValues::realdef;
      double hlCorSMn=DefaultValues::realdef;
      if (hlCirc>0. && ell>0.)
      {
        hlSMj=hlCirc*circToEll(ell);
        hlSMn=(1-ell)*hlSMj;
        hlGeo=sqrt(hlSMj*hlSMn);

        double smjSq=hlSMj*hlSMj-seeFact*seeing[extNo]*seeing[extNo];
        if (smjSq>0) {
            hlCorSMj=sqrt(smjSq);
        } else {
            hlCorSMj=0.;
        }
        double smnSq=hlSMn*hlSMn-seeFact*seeing[extNo]*seeing[extNo];
        if (smnSq>0) {
            hlCorSMn=sqrt(smnSq);
        } else {
            hlCorSMn=0.;
        }

      }

      // Assign to correct columns
      data.assign(hlGeoColNo, rowNo, hlGeo);
      data.assign(hlSMnColNo, rowNo, hlSMn);
      data.assign(hlSMjColNo, rowNo, hlSMj);
      data.assign(hlCorSMnColNo, rowNo, hlCorSMn);
      data.assign(hlCorSMjColNo, rowNo, hlCorSMj);

    }
  }

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

  //------------------------------------------------------
  double circToEll(double ell) {
    // Convert from elliptical half-light radius to circular
    // Form well fit by moffat profile.
    const double mofPars[3]={1.8243302,0.30914569,0.24304360};
    double axisRatio=1.-ell;
    axisRatio /= mofPars[1];
    return mofPars[0]/pow((1+axisRatio*axisRatio),mofPars[2]);

  }
private:
  const static double seeFact=0.50; // Calculate from data.
  FitsFile mCatFile;
  FitsFile mFile;
  int hlGeoColNo;
  int hlSMnColNo;
  int hlSMjColNo;
  int hlCorSMnColNo;
  int hlCorSMjColNo;
  int hlcColNo;
  int extColNo;
  int ellColNo;

  //

};
//-----------------------------------------------------------------------------
#endif
//-----------------------------------------------------------------------------
// Change log:
//
// 12-Nov-2009,  NJC: Original version.