source: Sophya/trunk/SophyaExt/FitsIOServer/fabtcolread.cc@ 1657

/* Lecteur de colonne de table Fits (binaire ou ASCII) avec buffer */
#include "machdefs.h"
#include <stdlib.h>
#include <stdio.h>
#include "pexceptions.h"
#include "fabtcolread.h"
//! Class for reading a column in a FITS ASCII or BINARY table

/*!
  \class SOPHYA::FitsABTColRead
  \ingroup FitsIOServer
  Class for reading a column in a FITS ASCII or BINARY table
  \verbatim
  Exemple:
  FitsABTColRead fbt("myfits.fits","BoloMuv_28",0,1000,1,3);
  fbt.SetDebug(3);
  fbt.Print(3);
  for(long i=0;i<fbt.GetNbLine();i++) {
    double x = fbt.Read(i);
    if(i%lpmod==0) cout<<i<<": "<<x<<endl;
  }
  \endverbatim
*/

//////////////////////////////////////////////////////////////
/*!
  Constructor.
  \verbatim
  fname : FITS file name to be read
  collabel : label of the column to be read
  ihdu : number of the HDU where the column is.
     if ihdu=0 or ihdu>nhdu first binary or ASCII table is taken
     Warning: ihdu = [1,nhdu]
  blen : read buffer length
  bsens : read buffer reading direction
          (bsens>0 forward, bsens<0 backward, bsens==0 centered)
  lp : debug level
  \endverbatim
*/
FitsABTColRead::FitsABTColRead(string fname,string collabel
                            ,int ihdu,long blen,long bsens,int lp)
{
Init(fname.c_str(),collabel.c_str(),-1,ihdu,blen,bsens,lp);
}

/*!
  Constructor.
  \verbatim
  Same as before but the column is identified by its number
  colnum : number of the column to be read
  WARNING: col = [0,ncol[
  \endverbatim
*/
FitsABTColRead::FitsABTColRead(string fname,int colnum
                            ,int ihdu,long blen,long bsens,int lp)
{
Init(fname.c_str(),"",colnum,ihdu,blen,bsens,lp);
}

/*! Constructor. see below */
FitsABTColRead::FitsABTColRead(const char * cfname,const char* collabel
                            ,int ihdu,long blen,long bsens,int lp)
{
Init(cfname,collabel,-1,ihdu,blen,bsens,lp);
}

/*! Constructor. see below */
FitsABTColRead::FitsABTColRead(const char * cfname,int colnum
                            ,int ihdu,long blen,long bsens,int lp)
{
Init(cfname,"",colnum,ihdu,blen,bsens,lp);
}

/*! Constructor by copy */
FitsABTColRead::FitsABTColRead(FitsABTColRead& fbt)
{
Init(fbt.GetFileName().c_str(),fbt.GetColLabel().c_str()
    ,fbt.GetColNum(),fbt.GetHDU()
    ,fbt.GetBLen(),fbt.GetBSens(),fbt.DbgLevel);
}

/*! Init routine called by the constructor */
void FitsABTColRead::Init(const char* fname,const char* collabel,int colnum
                        ,int ihdu,long blen,long bsens,int lp)
{
 // Parametres Generaux
 FitsFN = fname;
 ColLabel = collabel;
 ColTUnit = "";
 ColTForm = "";
 ColNum = colnum;
 ColTypeCode = 0;
 IHdu = ihdu;
 NHdu = 0;
 HduType = 0;
 NBcol = 0;
 NBline = 0;
 SetNulVal();
 SetDebug(lp);
 NFitsRead = 0;
 FitsPtr = NULL;
 LineDeb = LineFin = -1;
 Buffer = NULL;
 ChangeBuffer(blen,bsens);

 //////////////////////////
 // Ouverture du fichier //
 //////////////////////////

 int sta=0;
 if(FitsFN.size() <= 0 ) {
   IHdu = -1; Delete();
   throw ParmError("FitsABTColRead::Init: Fits file name error\n");
 }
 const char * cfitsfn = FitsFN.c_str();

 // Open fits
 if(fits_open_file(&FitsPtr,cfitsfn,READONLY,&sta)) {
   printerror(sta); Delete();
   throw NullPtrError("FitsABTColRead::Init: Error opening Fits file\n");
 }

 // Get number of hdu
 if(fits_get_num_hdus(FitsPtr,&NHdu,&sta)) {
   printerror(sta); Delete();
   throw NotAvailableOperation("FitsABTColRead::Init: Error getting NHdu\n");
 }
 if(DbgLevel>1) cout<<"FitsABTColRead::Init  NHdu="<<NHdu<<endl;
 if(NHdu<=0) {
   Delete();
   throw SzMismatchError("FitsABTColRead::Init: Bad NHdu\n");
 }

 // Get HDU for bin/ascii table
 // si IHdu <=0 || >NHdu on cherche la 1ere bin/ascii table
 // sinon on se positionne sur IHdu
 if(IHdu<=0 || IHdu>NHdu)
   for(int ihdu=1;ihdu<=NHdu;ihdu++) {
     if(fits_movabs_hdu(FitsPtr,ihdu,&HduType,&sta)) printerror(sta);
     if(DbgLevel>1) cout<<"...Init ihdu="
                        <<ihdu<<" HduType="<<HduType<<endl;
     if(HduType==BINARY_TBL || HduType==ASCII_TBL) {IHdu = ihdu; break;}
   }
 if(IHdu<=0 || IHdu>NHdu) {
   cout<<"NO BINARY or ASCII hdu found"<<endl;
   IHdu = 0; Delete();
   throw TypeMismatchExc("FitsABTColRead::Init: NO BINARY or ASCII hdu found\n");
 }
 if(fits_movabs_hdu(FitsPtr,IHdu,&HduType,&sta)) {
   printerror(sta); Delete();
   throw RangeCheckError("FitsABTColRead::Init: Error moving to requested HDU\n");
 }
 if(HduType!=BINARY_TBL && HduType!=ASCII_TBL) {
   Delete();
   throw TypeMismatchExc("FitsABTColRead::Init: HDU not ASCII/BINARY table\n");
 }

 // Get number of columns
 if(fits_get_num_cols(FitsPtr,&NBcol,&sta)) {
   printerror(sta); Delete();
   throw NotAvailableOperation("FitsABTColRead::Init: Error getting number of columns\n");
 }
 if(DbgLevel>1) cout<<"...Init  NBcol="<<NBcol<<endl;
 if(NBcol<1) {
   Delete();
   throw RangeCheckError("FitsABTColRead::Init: Bad number of colums\n");
 }

 // Get number of rows
 if(fits_get_num_rows(FitsPtr,&NBline,&sta)) {
   printerror(sta); Delete();
   throw NotAvailableOperation("FitsABTColRead::Init: Error getting number of rows\n");
 }
 if(DbgLevel>1) cout<<"...Init  NBline="<<NBline<<endl;
 if(NBline<1) {
   Delete();
   throw RangeCheckError("FitsABTColRead::Init: Bad number of rows\n");
 }

 // Get column number
 char labelcol[256];
 if(ColLabel.size() > 0) {
   strcpy(labelcol,ColLabel.c_str());
   if(fits_get_colnum(FitsPtr,CASESEN,labelcol,&ColNum,&sta)) {
     printerror(sta); Delete();
     throw NotAvailableOperation("FitsABTColRead::Init: Error getting column name\n");
   }
   ColNum--;  // Convention [0,ncol[
 }
 if(DbgLevel>1) cout<<"...Init  ColNum="<<ColNum<<endl;
 if(ColNum<0 || ColNum>=NBcol) {
   Delete();
   throw RangeCheckError("FitsABTColRead::Init: Bad column number\n");
 }

 // Get column type
 int ColTypeCode;
 if(fits_get_coltype(FitsPtr,ColNum+1,&ColTypeCode,NULL,NULL,&sta)) {
   printerror(sta); Delete();
   throw ParmError("FitsABTColRead::Init: Error getting column type\n");
 }
 if(DbgLevel>1) cout<<"...Init ColTypeCode="<<ColTypeCode<<endl;
 if(ColTypeCode==TSTRING || ColTypeCode==TCOMPLEX ||  ColTypeCode==TDBLCOMPLEX) {
   Delete();
   throw ParmError("FitsABTColRead::Init: Selected column is not Numerical\n");
 }

 // Get column name back, tunit, tform
 char tunit[64], tform[64];
 int rc=0;
 if(HduType==BINARY_TBL) {
   fits_get_bcolparms(FitsPtr,ColNum+1,labelcol,tunit,tform,NULL,NULL,NULL,NULL,NULL,&sta);
 } else {
   fits_get_acolparms(FitsPtr,ColNum+1,labelcol,NULL,tunit,tform,NULL,NULL,NULL,NULL,&sta);
 }
 if(rc) {
   printerror(sta); Delete();
   throw RangeCheckError("FitsABTColRead::Init: Error getting the column caracteristics\n");
 }
 ColLabel = labelcol;
 ColTUnit = tunit;
 ColTForm = tform;

 if(DbgLevel)
   cout<<"FitsABTColRead::Init Num="<<ColNum<<" Label="<<ColLabel
       <<" TypeCode="<<ColTypeCode
       <<" TUnit="<<ColTUnit<<" TForm="<<ColTForm<<endl;

}

/*! Destructor. */
FitsABTColRead::~FitsABTColRead()
{
 Delete();
}


//////////////////////////////////////////////////////////////
/*! Change the buffer caracteristiques (see creator)  */
void FitsABTColRead::ChangeBuffer(long blen,long bsens)
{
 long oldnbuffer = NBuffer;

 // Compute buffer caracteristics
 BuffLen = (blen<=0)? 1: blen;
 BuffSens = bsens;
 NBuffer = BuffLen;
 if(bsens==0 && NBuffer%2==0) NBuffer++;

 // De-allocate if necessary
 if(Buffer!=NULL && oldnbuffer!=NBuffer)
   {delete [] Buffer; Buffer=NULL;}

 // Re-allocate
 if(Buffer==NULL) Buffer = new double[NBuffer];

 // Tell program that nothing is into buffer
 LineDeb = LineFin = -1;
}

/*! Delete called by the destructor */
void FitsABTColRead::Delete()
{
 if(Buffer!=NULL) {delete [] Buffer; Buffer=NULL;}
 LineDeb = LineFin = -1;
 int sta = 0;
 if(fits_close_file(FitsPtr,&sta)) printerror(sta);
 FitsPtr = NULL;
}

/////////////////////////////////////////////////
/*!
  Read row "n" and return the value cats into double
  \verbatim
  WARNING: row = [0,NRows[
  \endverbatim
*/
double FitsABTColRead::Read(long n)
// Attention: n=nline [0,NBline[, cfistio veut [1,NBline]
// Attention: colnum  [0,NBcol[ , cfistio veut [1,NBcol]
{
 int sta=0,anynul;
 if(n<0 || n>=NBline)
   throw RangeCheckError("FitsABTColRead::Read try to read outside line range\n");

 // Pas de bufferisation, on lit betement
 if(NBuffer==1) {
   NFitsRead++;
   fits_read_col_dbl(FitsPtr,ColNum+1,n+1,1,1,NulVal,Buffer,&anynul,&sta);
   if(sta) {
     printerror(sta);
     throw NotAvailableOperation("FitsABTColRead::Read: Error Reading Fits file\n");
   }
   return Buffer[0];
 }

 // Gestion avec bufferisation
 if(!Buffer) {
   cout<<"FitsABTNtuIntf::Read Buffer not allocated"<<endl;
   return 0;
 }
 if(n<LineDeb || n>LineFin) {
   NFitsRead++;
   long row1,row2,nrow;
   if(BuffSens>0) { // Cas remplissage forward
     row1 = n+1;
     row2 = row1+NBuffer-1; if(row2>NBline) row2 = NBline;
   } else if(BuffSens<0) { // Cas remplissage backward
     row2 = n+1;
     row1 = row2-NBuffer+1; if(row1<1) row1 = 1;
   } else { // Cas remplissage centre
     row1 = n+1 - NBuffer/2; if(row1<1) row1 = 1;
     row2 = n+1 + NBuffer/2; if(row2>NBline) row2 = NBline;
   }
   nrow = row2 - row1 + 1;
   LineDeb = row1-1; LineFin = row2-1;
   //cout<<"DBG-FitsRead: row1="<<row1<<" row2="<<row2<<" nrow="<<nrow
   //    <<" LineDeb,Fin="<<LineDeb<<","<<LineFin<<endl;
   fits_read_col_dbl(FitsPtr,ColNum+1,row1,1,nrow,NulVal,Buffer,&anynul,&sta);
   if(sta) {
     printerror(sta);
     LineDeb = LineFin = -1;
     throw NotAvailableOperation("FitsABTColRead::Read: Error Reading Fits file\n");
   }
 }

 long ibuf = n-LineDeb;
 return Buffer[ibuf];
}

/*!
  Read rows from "n1" to "n2" and return the values cats in a TVector of double
  \verbatim
  There are n2-n1+1 values returned
  WARNING: row = [0,NRows[
  \endverbatim
*/
void FitsABTColRead::Read(long n1,long n2,TVector<double>& data)
{
 if(n1<0 || n1>=NBline || n2<0 || n2>=NBline || n1>n2)
   throw RangeCheckError("FitsABTColRead::Read TVector bad requested line range\n");

 sa_size_t n = n2-n1+1;
 if(data.Size()<n) data.SetSize(n);
 // Il faut faire mieux mais comment ????
 for(long i=n1;i<=n2;i++) data(i-n1) = Read(i);
}

/////////////////////////////////////////////////
void FitsABTColRead::printerror(int sta) const
{
 int stat = sta;
 fits_report_error(stdout,stat);
 fflush(stdout);
 return;
}

/*! Print on stream os */
void FitsABTColRead::Print(ostream& os,int lp) const
{
 os<<"FitsABTColRead:Print ("<<BuffLen<<","<<BuffSens<<","<<NulVal<<")"
   <<" ncols="<<NBcol<<" nrows="<<NBline;
 if(lp>0) os<<" NRead="<<NFitsRead;
 os<<"\n... "<<FitsFN<<"["<<IHdu<<"/"<<NHdu<<"]"
   <<"\n... Label["<<ColNum<<"]="<<ColLabel
   <<" TypeCode="<<ColTypeCode
   <<" TUnit="<<ColTUnit<<" TForm="<<ColTForm
   <<endl;
}