RaDecByAst.hxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// CVS: $Revision: 4695 $, $Date: 2008-04-09 16:18:42 +0100 (Wed, 09 Apr 2008) $, $Author: EckhardSutorius $
//-----------------------------------------------------------------------------
#ifndef RADECBYAST_HXX
#define RADECBYAST_HXX

// We will make use of AST for operations involving WCS
extern "C"
{
  #include "ast.h"
}

#include <cmath>

#include "Converter.h"
#include "FitsFile.h"
#include "Logger.h"
#include "DataMethod.hxx"
#include "StringOps.h"
//-----------------------------------------------------------------------------
template<typename DataType>
class RaDecByAst : public DataMethod<DataType>
{
public:
  RaDecByAst(const FitsFile& aFile)
    : mFile(aFile)
    { DataMethod<DataType>::methodName = "RaDecByAst"; }

  virtual ~RaDecByAst() { }

  void setSource(IntMap& tbls, IntMap& cols, StringMap units)
  {
    extnumColno = cols[0];
    xColno = cols[1];
    yColno = cols[2];
    extnumTblno = tbls[0];
    xTblno = tbls[1];
    yTblno = tbls[2];
  }

  void setTarget(IntMap& cols, StringMap units)
  {
    raColno  = cols[0];
    decColno = cols[1];

    // AST will calculate angles in radians, so we need the conversion
    // factors for the required units in the schema
    Converter conv;
    raFact  = conv.getFactor("RADIANS", units[0]);
    decFact = conv.getFactor("RADIANS", units[1]);
  }

  void doit2(TableData<DataType>& inpData, TableData<DataType>& data, int row1, int row2)
  {
    double pi = M_PI;
    string filename = mFile.getFileName();
    double crval1,crval2;

    // Set up AST
    astBegin;

    // Create a new AST fits channel
    AstFitsChan *fitschan = astFitsChan( NULL, NULL, "" );
    AstFrameSet *wcsinfo0 = 0;
    AstFrameSet *wcsinfo1 = 0;
    AstFrameSet *wcsinfo2 = 0;
    AstFrameSet *wcsinfo3 = 0;

    for (int hduNo = 2; hduNo <= mFile.getNumHdus(); ++hduNo)
    {
      int ext = hduNo - 2;

      // Move to extension on the FITS file
      mFile.movabsHdu(hduNo);

      // 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::cCRVAL1,crval1);
      mFile.readKey(FitsFile::cCRVAL2,crval2);
      if (crval1 == 0 && crval2 == 0) {
          log.addError("Could not read WCS info for " + filename);
      } else {
          if (ext == 0) {
              wcsinfo0 = (AstFrameSet*)astRead(fitschan);
          } else if (ext == 1) {
              wcsinfo1 = (AstFrameSet*)astRead(fitschan);
          } else if (ext == 2) {
              wcsinfo2 = (AstFrameSet*)astRead(fitschan);
          } else if (ext == 3) {
              wcsinfo3 = (AstFrameSet*)astRead(fitschan);
          }
      }
    }

    // Use these arrays for AST input
    double xpix[1], ypix[1], xsky[1], ysky[1];

    for (int rowNo = row1; rowNo < row2; ++rowNo)
    {
      // Get the FITS extension number for the data in this row. Note
      // that the extension number is assumed to be in the fitsio
      // convention
      Numeric value;
      data.value(extnumColno, rowNo, value);

      int extNum = (int)value;

      if (crval1 != 0 || crval2 != 0)
      {
        // Input X and Y
        //double xpix, ypix;
        inpData.value(xColno, rowNo, xpix[0]);
        inpData.value(yColno, rowNo, ypix[0]);

        // Get the sky coordinates
        // Note these will be in radians. Angles can also be negative.
        if (extNum == 2)
        {
          astTran2(wcsinfo0, 1, xpix, ypix, 1, xsky, ysky);
        }
        else if (extNum == 3)
        {
          astTran2(wcsinfo1, 1, xpix, ypix, 1, xsky, ysky);
        }
        else if (extNum == 4)
        {
          astTran2(wcsinfo2, 1, xpix, ypix, 1, xsky, ysky);
        }
        else if (extNum == 5)
        {
          astTran2(wcsinfo3, 1, xpix, ypix, 1, xsky, ysky);
        }

        // Make sure all RAs are in the range 0 to 2*pi
        if (xsky[0] < 0) xsky[0] += 2 * pi;

        Numeric ra = raFact * xsky[0];
        Numeric dec = decFact * ysky[0];

        data.assign(raColno, rowNo, ra);
        data.assign(decColno, rowNo, dec);
      }
    }
    // Finshed with AST
    astEnd;
  }

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

private:
  FitsFile mFile;
  int extnumColno, raColno, decColno, xColno, yColno;
  int extnumTblno, xTblno, yTblno;
  double decFact, raFact;
  Logger log;
};
//-----------------------------------------------------------------------------
#endif
//-----------------------------------------------------------------------------
// Change log:
//
//  7-Jun-2005, ETWS: Original version.
//  7-Apr-2008, ETWS: Upgraded to use new detection table layout.

---