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

/* 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);
  // Read element by element
  for(long i=0;i<fbt.GetNbLine();i++) {
    double x = fbt.Read(i);
    if(i%lpmod==0) cout<<i<<": "<<x<<endl;
  }
  // Read into a vector
  TVector<double> data;
  long n = fbt.Read(32,50,data);
    cout<<"Number of values read: "<<n<<endl;
  data.ReSize(100);
  n = fbt.Read(10,-1,data);
    cout<<"Number of values read: "<<n<<endl;
  \endverbatim
*/

//////////////////////////////////////////////////////////////
/*!
  Constructor.
  \param fname : FITS file name to be read
  \param collabel : label of the column to be read
  \param ihdu : number of the HDU where the column is.
  \param blen : read buffer length
  \param bsens : buffer reading direction
  \param lp : debug level
  \verbatim
  - if ihdu=0 or ihdu>nhdu first binary or ASCII table is taken
  - bsens>0    read forward
    bsens<0    read backward
    bsens==0   read centered
  \endverbatim
  \warning ihdu = [1,nhdu]
*/
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.
  Same as before but the column is identified by its column number
  \param colnum : number of the column to be read
  \warning col = [0,ncol[
*/
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[128];
 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
                         || ColTypeCode<0 ) {
   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();
}

//////////////////////////////////////////////////////////////
double FitsABTColRead::ReadKey(char *keyname)
{
 if(keyname==NULL) return 0.;
 int sta=0; double val;
 if(fits_read_key(FitsPtr,TDOUBLE,keyname,&val,NULL,&sta))
   printerror(sta);
 return val;
}

//////////////////////////////////////////////////////////////
/*! 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) {
   // On des-alloue si pas assez de place
   // ou si l'ancienne place est beaucoup trop grande (>25%)
   if(oldnbuffer<NBuffer || (oldnbuffer>NBuffer+long(0.25*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 into a double
  \warning be carefull for the range: row = [0,NRows[
  \return value in double
  \param n : number of the row to be read.
  \verbatim
  usebuffer == true  : use read optimisation with bufferisation
            == false : no optimisation with bufferisation
                       just read one value
  \endverbatim
*/
double FitsABTColRead::Read(long n,bool usebuffer)
// Attention: n=nline [0,NBline[, cfistio veut [1,NBline]
// Attention: colnum  [0,NBcol[ , cfistio veut [1,NBcol]
{
 int sta=0;
 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 || !usebuffer) {
   NFitsRead++;
   double val;
   fits_read_col(FitsPtr,TDOUBLE,ColNum+1,n+1,1,1,NULL,&val,NULL,&sta);
   if(sta) {
     printerror(sta);
     throw NotAvailableOperation("FitsABTColRead::Read: Error Reading Fits file\n");
   }
   // On ne remplit Buffer[0] que si on a choisit
   // un mode de lecture non bufferise (n==1) DES LE DEBUT.
   // Si on a initialement choisit un mode bufferise (avec n>1),
   // Buffer contient les valeurs chargees auparavent.
   // Il ne faut pas faire {Buffer[0]=val; LineDeb=LineFin=n;}
   // car on perd l'info de ces valeurs.
   if(NBuffer==1) {Buffer[0]=val; LineDeb=LineFin=n;}
   return val;
 }

 // Gestion avec bufferisation
 if(!Buffer)
   throw RangeCheckError("FitsABTColRead::Read: Buffer not allocated\n");
 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(FitsPtr,TDOUBLE,ColNum+1,row1,1,nrow,NULL,Buffer,NULL,&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 into TVector of double
  \return NREAD the number of values read (n2-n1+1).
  \warning row = [0,NRows[, the routine read [n1,n2]
  \verbatim
  - if n2<0 then read [n1,n2] where "n2=min(n1+vector_size-1,nrows-1)"
  - Last row read is ALWAYS: "n2 = n1 + NREAD -1"
  - The TVector is never resized if not necessary
  -------------------------------------------------------------------------
  - ex: suppose the column table contains 10 elements: nrows=10, rows=[0,9]

    TVector<double> V(5);
      bt.Read(3,5,V)  -> read rows=3,4,5     -> V.Size()==5 -> return 3
      bt.Read(3,-1,V) -> read rows=3,4,5,6,7 -> V.Size()==5 -> return 5
      bt.Read(7,-1,V) -> read rows=7,8,9     -> V.Size()==5 -> return 3
      bt.Read(2,-1,V) -> read rows=2,3,4,5,6 -> V.Size()==5 -> return 5
      bt.Read(-1,5,V) -> throw exception

    TVector<double> V(5);
      bt.Read(3,99,V) -> read rows=3,4,5,6,7,8,9 -> V.Size()==7 -> return 7

    TVector<double> V(5);
      bt.Read(2,8,V)  -> read rows=2,3,4,5,6,7,8 -> V.Size()==7 -> return 7

    TVector<double> V;
      bt.Read(3,5,V)  -> read rows=3,4,5 -> V.Size()==3 -> return 3

    TVector<double> V;
      bt.Read(3,-1,V) -> throw exception
  -------------------------------------------------------------------------
  \endverbatim
*/
long FitsABTColRead::Read(long n1,long n2,TVector<double>& data)
{
 if(n1<0 || n1>=NBline)
   throw RangeCheckError("FitsABTColRead::Read TVector bad requested 1srt line \n");
 if(data.Size()<=0 && n2<n1)
   throw RangeCheckError("FitsABTColRead::Read TVector bad requested 2sd line \n");
 if(n2<0) n2 = n1 + data.Size()-1;
 if(n2>=NBline) n2 = NBline-1;

 sa_size_t nread = n2-n1+1;
 if(data.Size()<nread) data.SetSize(nread);

 //for(long i=n1;i<=n2;i++) data(i-n1) = Read(i);
 int sta=0;
 fits_read_col(FitsPtr,TDOUBLE,ColNum+1,n1+1,1,nread,NULL,data.Data(),NULL,&sta);
 if(sta) {
   printerror(sta);
   throw NotAvailableOperation("FitsABTColRead::Read_TVector<double>: Error Reading Fits file\n");
 }

 return nread;
}

/*! idem before but for TVector of float */
long FitsABTColRead::Read(long n1,long n2,TVector<float>& data)
{
 if(n1<0 || n1>=NBline)
   throw RangeCheckError("FitsABTColRead::Read TVector bad requested 1srt line \n");
 if(data.Size()<=0 && n2<n1)
   throw RangeCheckError("FitsABTColRead::Read TVector bad requested 2sd line \n");
 if(n2<0) n2 = n1 + data.Size()-1;
 if(n2>=NBline) n2 = NBline-1;

 sa_size_t nread = n2-n1+1;
 if(data.Size()<nread) data.SetSize(nread);

 //for(long i=n1;i<=n2;i++) data(i-n1) = Read(i);
 int sta=0;
 fits_read_col(FitsPtr,TFLOAT,ColNum+1,n1+1,1,nread,NULL,data.Data(),NULL,&sta);
 if(sta) {
   printerror(sta);
   throw NotAvailableOperation("FitsABTColRead::Read_TVector<float>: Error Reading Fits file\n");
 }

 return nread;
}

/*! idem before but for TVector of int_4 */
long FitsABTColRead::Read(long n1,long n2,TVector<int_4>& data)
{
 if(n1<0 || n1>=NBline)
   throw RangeCheckError("FitsABTColRead::Read TVector bad requested 1srt line \n");
 if(data.Size()<=0 && n2<n1)
   throw RangeCheckError("FitsABTColRead::Read TVector bad requested 2sd line \n");
 if(n2<0) n2 = n1 + data.Size()-1;
 if(n2>=NBline) n2 = NBline-1;

 sa_size_t nread = n2-n1+1;
 if(data.Size()<nread) data.SetSize(nread);

 //for(long i=n1;i<=n2;i++) data(i-n1) = Read(i);
 int sta=0;
 int T = (sizeof(long)==4) ? TLONG: TINT;
 fits_read_col(FitsPtr,T,ColNum+1,n1+1,1,nread,NULL,data.Data(),NULL,&sta);
 if(sta) {
   printerror(sta);
   throw NotAvailableOperation("FitsABTColRead::Read_TVector<float>: Error Reading Fits file\n");
 }

 return nread;
}

/////////////////////////////////////////////////
/*!
  Return the number of the first row where "val1"<=val<="val2" starting at row "rowstart"
  \verbatim
  - The search is performed from "rowstart" to the end
      in ascending order (from "rowstart" to nrows).
  - Warning: "rowstart<0" means "rowstart==0" (search all the table column)
             That is the default
  \endverbatim
  \return <0 means not found
*/
long FitsABTColRead::FirstRow(double val1,double val2,long rowstart)
{
 long row = -1;
 if(NBline==0) return row;
 // Change buffer for efficiency
 long bsens=BuffSens; bool bchange=false;
 if(bsens<=0) {ChangeBuffer(BuffLen,1); bchange=true;}
 if(rowstart<0) rowstart = 0;
 if(rowstart>=NBline) rowstart = NBline-1;
 for(long i=rowstart;i<NBline;i++) {
   double val = Read(i);
   if(val<val1 || val>val2) continue;
   row = i;
   break;
 }
 if(bchange) ChangeBuffer(BuffLen,bsens);
 return row;
}

/*!
  Return the number of the first row where val1<=val<=val2 starting at row rowstart
  \return <0 means not found
  \verbatim
  - The search is performed from "rowstart" to the beginning
      in descending order (from "rowstart" to 0).
  - Warning: "rowstart<0" means "rowstart==nrows-1" (search all the table column)
             That is the default
  \endverbatim
*/
long FitsABTColRead::LastRow(double val1,double val2,long rowstart)
{
 long row = -1;
 if(NBline==0) return row;
 // Change buffer for efficiency
 long bsens=BuffSens; bool bchange=false;
 if(bsens>=0) {ChangeBuffer(BuffLen,-1); bchange=true;}
 if(rowstart<0 || rowstart>=NBline) rowstart = NBline-1;
 for(long i=rowstart;i>=0;i--) {
   double val = Read(i);
   if(val<val1 || val>val2) continue;
   row = i;
   break;
 }
 if(bchange) ChangeBuffer(BuffLen,bsens);
 return row;
}

/////////////////////////////////////////////////
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;
}