HlRadCirc.hxx

Go to the documentation of this file.

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

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

#include <algorithm>
#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 HlRadCirc : public DataMethod<DataType>
{
public:
    HlRadCirc(const FitsFile& aCatFile, const FitsFile& aFile)
      : mCatFile(aCatFile),mFile(aFile)
   {
     DataMethod<DataType>::methodName = "HlRadCirc";
   }

  virtual ~HlRadCirc() { }

  //---------------------------------------------------------------------------
  void setSource(IntMap& tbls, IntMap& cols, StringMap units)
  {
    afColNo[0]=0;
    afeColNo[0]=0;

    for (int i=0;i<NOAPERS-1;++i)
    {
      afColNo[i+1]  = cols[2*i];
      afeColNo[i+1] = cols[2*i+1];
    }
    pfColNo    = cols[26];
    extColNo   = cols[27];

    //xpsColNo   = cols[30]; // Both x and y?
    //ypsColNo   = cols[31];


  }
  //---------------------------------------------------------------------------
  void setTarget(IntMap& cols, StringMap units)
  {
    hlCircColNo   = cols[0];
    hlCErrColNo   = cols[1];

  }
  //---------------------------------------------------------------------------
  void doit(TableData<DataType>& data, int row1, int row2)
  {
    Logger log;

    double apertures[NOAPERS]={0.,0.5,sqrt(0.5),1.,sqrt(2.),2.,2.*sqrt(2.),4.,5.,6.,7.,8.,10.,12.}; // in arcsec -- use coreRadius to convert.
    double rCore[mCatFile.getNumHdus()];
    double crval1;
    double crval2;
    int npixx,npixy;

    double centre[2];
    double sizeFact[2];
    double pixelScale[mCatFile.getNumHdus()];


    // 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("RCORE", rCore[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" );

      // Calc xpixSize,ypixSize

      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);
      pixelScale[extNo]=sqrt(xPixSize*yPixSize);

    }

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

      // Input aperture fluxes and errors
      double aperFlux[NOAPERS];
      double aperFluxErr[NOAPERS];
      //double xPixSize,yPixSize;

      double petroFlux;
      double totalFlux;
      // Use petrosian
      data.value(pfColNo, rowNo, petroFlux);

      double hlCirc=DefaultValues::realdef;
      double hlCErr=DefaultValues::realdef;
      // Get hlCirc

      // Put distort correction in magnitudes in separate
      // if no distortion correction in mags, none in radii

      // Get 5 nearest apertures
      // Fit quadratic using Least Squares

      if (petroFlux>0.)
      {
        int bestAper=0;
        bool isAperSet=false;
        totalFlux=petroFlux/0.9; // Over estimate for disks, under for ells

        for (int i=1;i<NOAPERS;++i)
        {
          data.value(afColNo[i], rowNo, aperFlux[i]);
          data.value(afeColNo[i], rowNo, aperFluxErr[i]);
          // Find radius where aperFlux~0.5*kronFlux
          if (aperFlux[i]>0.5*totalFlux && !isAperSet)
          {
            bestAper=i;
            isAperSet=true;

          }

        }
        aperFlux[0]=0.;
        aperFluxErr[0]=aperFluxErr[1];
        if (!isAperSet) bestAper=NOAPERS-1;

        // Find five best apertures
        int minAper = max(bestAper-int(floor(NDATAPNTS/2.)),0);
        int maxAper = min(minAper+NDATAPNTS,NOAPERS-1);
        if ((maxAper-minAper)<NDATAPNTS) minAper=maxAper-NDATAPNTS;

        // Now Get data, set up matrix for LSQ via SVN
        double matA[NDATAPNTS][NPOLY];
        double vecB[NDATAPNTS];
        double matW[NPOLY];
        double matV[NPOLY][NPOLY];
        double work[NPOLY];
        double xfit[NPOLY];
        double matC[NPOLY][NPOLY];
        int jstat;
        double minRadius=rCore[extNo]*apertures[minAper];
        double maxRadius=rCore[extNo]*apertures[minAper+NDATAPNTS-1];


        for (int j=0;j<NDATAPNTS;++j)
        {

            double radius;
            radius=rCore[extNo]*apertures[minAper+j];
            if (j==0) minRadius=radius;
            if (j==NDATAPNTS-1) maxRadius=radius;
            double frac=0.;
            double fracErr=1.;

            if (aperFlux[minAper+j]>=0.) {
                frac=aperFlux[minAper+j]/totalFlux;
                if (aperFluxErr[minAper+j]>0.) {
                    fracErr=aperFluxErr[minAper+j]/totalFlux;
                }
            }
            double weight=1./fracErr;
            matA[j][0]=weight;
            vecB[j]=frac*weight;
            for (int k=0;k<NPOLY-1;++k) {
                matA[j][k+1]=pow(radius,k+1)*weight;
            }
        }

         // Use SVD. Ax=0
        //
        cslaSvd(NDATAPNTS,NPOLY,NDATAPNTS,NPOLY, (double *) matA,matW, (double *) matV, work, &jstat);
         // @todo: set smallest SVD components to 0.
        // ratio should be no larger than 10^10 if dp.
        double maxSvd=*max_element(matW,matW+NPOLY);
        for (int j=0;j<NPOLY;++j)
        {
            if (matW[j]<=minSvdRatio*maxSvd) matW[j]=0;

        }

        cslaSvdsol(NDATAPNTS,NPOLY,NDATAPNTS,NPOLY,vecB,(double *) matA,matW,(double *) matV, work,xfit);
        cslaSvdcov(NPOLY,NPOLY,NPOLY,matW,(double *) matV, work, (double *) matC);
        //
        // Now get polynomial and find roots
        // Numerical recipes: zroots
        // Make sure roots are sensible
        double roots[2];
        // Find when polynomial fits data
        xfit[0]-=0.5;
        // Use Newton Raphson? --> general
        quadroots(xfit,roots); // replace by newraph
        // Find most likely root. Should be within range of data
        // Deblended objs and ones where there is a sudden turn down just
        // after half-light radius have larger hlradius than petroAper.
        // Do a few more checks
        bool hasMatch=false;

        for (int i=0;i<NPOLY-1;++i) {
          if (roots[i]>=minRadius && roots[i]<=maxRadius && !hasMatch) {
            hlCirc=roots[i];
            hasMatch=true;
          }
        }
        double sumCov=0.;
        double dfda[NPOLY];
        double dfdr=0;
        dfda[0]=1.;

        for (int i=0;i<NPOLY-1;++i)
        {
          dfda[i+1]=pow(hlCirc,i+1);
          dfdr+=pow(hlCirc,i)*(i+1)*xfit[i+1];
        }
        for (int i=0;i<NPOLY;++i) {
          for (int j=0;j<NPOLY;++j) {
              sumCov+=dfda[i]*dfda[j]*matC[i][j];
          }
        }
        double ferr=sqrt(sumCov);
        // Add pixel scale error in quarature
        if (hlCirc<0.) {
          hlCErr=DefaultValues::realdef;
        }  else if (abs(dfdr)>0.) {
          hlCirc=pixelScale[extNo]*hlCirc;
          hlCErr=pixelScale[extNo]*sqrt((ferr/dfdr)*(ferr/dfdr)+0.5*0.5);
        } else {
          hlCirc=pixelScale[extNo]*hlCirc;
          hlCErr=pixelScale[extNo]*0.5;
        }


      }

      data.assign(hlCircColNo, rowNo, hlCirc);
      data.assign(hlCErrColNo, rowNo, hlCErr);

    }

  }

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

  void quadroots(double a[3], double roots[2]) {
    // Simple method to find quadratic roots
    double sq=a[1]*a[1]-4.*a[0]*a[2];
    if (sq>=0) {
        roots[0]=(-a[1]+sqrt(sq))/(2.*a[2]);
        roots[1]=(-a[1]-sqrt(sq))/(2.*a[2]);
    } else {
        roots[0]=DefaultValues::realdef;
        roots[1]=DefaultValues::realdef;
    }
  }
  //---------------------------------------------------------------------------
  double newrph(double poly[], double hlrE, double mnr, double mxr) {
    const double finalErr=1.0e-4;
    const int maxIter=5;
    double error=1.;
    int it=0;
    double hlr=hlrE;
    double dr=0;
    while (error>finalErr and it<maxIter) {
        double func=0;
        double dfdr=0;
        for (int i=0;i<NPOLY;++i) {
            func=poly[i]*pow(hlr,i);
            if (i>0) {
                dfdr+=i*poly[i]*pow(hlr,(i-1));
            }
        }
        dr=func/dfdr;
        error=sqrt(func);
        hlr-=dr;

        it++;
    }
    // @todo: think about this
    if (hlr<mnr) hlr=mnr;
    if (hlr>mxr) hlr=mxr;
    return hlr;
  }
private:
  const static unsigned short int NOAPERS=14;
  const static unsigned short int NDATAPNTS=5;
  const static unsigned short int NPOLY=3; // Quadratic
  const static double minSvdRatio=1.0e-10;
  FitsFile mCatFile;
  FitsFile mFile;

  int hlCircColNo;
  int hlCErrColNo;
  int afColNo[NOAPERS];
  int afeColNo[NOAPERS];
  int pfColNo;
  int extColNo;


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