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

/* 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 opening a FITS file for reading

/*!
  \class SOPHYA::FitsOpenFile
  \ingroup FitsIOServer
  Class for opening a FITS file for reading
*/

/*!
  Constructor.
  \param fname : FITS file name to be opened for reading
*/
FitsOpenFile::FitsOpenFile(string fname)
{
  Init(fname.c_str());
}

/*!
  Constructor.
  \param cfname : FITS file name to be opened for reading
*/
FitsOpenFile::FitsOpenFile(const char *cfname)
{
  Init(cfname);
}

/*!
  Constructor by default.
*/
FitsOpenFile::FitsOpenFile()
{
  FitsFN = "";
  NHdu = IHdu = HduType = 0;
  HasBeenPos = false;
  FitsPtr = NULL;
}

/*!
  Constructor by copy.
*/
FitsOpenFile::FitsOpenFile(FitsOpenFile& fof)
{
  Init(fof.FileName().c_str());
}

/*!
  Destructor.
*/
FitsOpenFile::~FitsOpenFile()
{
  Delete();
  FitsFN = "";
  NHdu = IHdu = HduType = 0;
  HasBeenPos = false;
}

/*!
  Delete routine called by the destructor.
*/
void FitsOpenFile::Delete(void)
{
  if(FitsPtr != NULL) {
    int sta = 0;
    if(fits_close_file(FitsPtr,&sta)) printerror(sta);
    FitsPtr = NULL;
  }
}

/*! Init routine called by the constructor */
void FitsOpenFile::Init(const char* fname)
{
  // Parametres Generaux
  FitsFN = fname;
  NHdu = IHdu = HduType = 0;
  HasBeenPos = false;
  FitsPtr = NULL;

  // Ouverture du fichier
  if(FitsFN.size() <= 0 )
     throw ParmError("FitsOpenFile::Init: Fits file name error\n");

  int sta = 0;
  if(fits_open_file(&FitsPtr,FitsFN.c_str(),READONLY,&sta)) {
    printerror(sta);
    FitsPtr = NULL;
    throw NullPtrError("FitsOpenFile::Init: Error opening Fits file\n");
  }

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

  MoveToHDU(1);
}

/*! Move to an HDU
\param ihdu: hdu number to move
\warning ihdu = [1,nhdu]
\return 0 if positionning failed, ihdu if success
*/
int FitsOpenFile::MoveToHDU(int ihdu)
{
 if(FitsPtr==NULL)
   throw NullPtrError("FitsOpenFile::MoveToHDU: no fits file open FitsPtr==NULL\n");
 int ihdusave = IHdu;
 if(ihdu<=0) ihdu=1; if(ihdu>NHdu) ihdu=NHdu;
 int sta=0;
 if(fits_movabs_hdu(FitsPtr,ihdu,&HduType,&sta)) {
   printerror(sta);
   // On se repositionne ou on etait
   fits_movabs_hdu(FitsPtr,ihdusave,&HduType,&sta);
   IHdu = ihdusave;
 } else IHdu = ihdu;
 return IHdu;
}

/*! Move to the first HDU of a certain type
\param hdutype: type of the hdu
\param hdudeb: start at that hdu
\return the type of HDU the file is positionned
*/
int FitsOpenFile::MoveToFirst(int hdutype,int ihdudeb)
{
 if(ihdudeb<=0) ihdudeb=1; if(ihdudeb>NHdu) ihdudeb=NHdu;
 int ihdusave = IHdu;
 for(int ihdu=ihdudeb;ihdu<=NHdu;ihdu++) {
   MoveToHDU(ihdu);
   if(HduType==hdutype) break;
 }
 // Si echec, on se repositionne ou on etait
 if(HduType!=hdutype) MoveToHDU(ihdusave);
 return HduType;
}

/*! Move to the last HDU of a certain type
\param hdutype: type of the hdu
\param hdudeb: stop at that hdu
\return the type of HDU the file is positionned
*/
int FitsOpenFile::MoveToLast(int hdutype,int ihdudeb)
{
 if(ihdudeb<=0) ihdudeb=1; if(ihdudeb>NHdu) ihdudeb=NHdu;
 int ihdusave = IHdu;
 for(int ihdu=NHdu;ihdu>=ihdudeb;ihdu--) {
   MoveToHDU(ihdu);
   if(HduType==hdutype) break;
 }
 // Si echec, on se repositionne ou on etait
 if(HduType!=hdutype) MoveToHDU(ihdusave);
 return HduType;
}

/*! Print */
void FitsOpenFile::Print(void)
{
 cout<<"FitsOpenFile::Print: "<<FitsFN
     <<" hdu="<<IHdu<<"/"<<NHdu<<" type="<<HduType
     <<" hasbeenpos="<<HasBeenPos<<endl;
}

//////////////////////////////////////////////////////////////
//// Methodes statiques
//////////////////////////////////////////////////////////////
/*!
  Read a fitsheader key into double
  \param fitsptr : cfitio pointer to Fits file
  \param keyname : name of the key
  \return value into double
*/
double FitsOpenFile::ReadKey(fitsfile *fitsptr,char *keyname)
{
 if(keyname==NULL || fitsptr==NULL) return 0.;
 int sta=0; double val=0.;
 if(fits_read_key(fitsptr,TDOUBLE,keyname,&val,NULL,&sta))
   printerror(sta);
 return val;
}

/*!
  Read a fitsheader key into long
  \param fitsptr : cfitio pointer to Fits file
  \param keyname : name of the key
  \return value into long
*/
long FitsOpenFile::ReadKeyL(fitsfile *fitsptr,char *keyname)
{
 if(keyname==NULL || fitsptr==NULL) return 0;
 int sta=0; long val=0;
 if(fits_read_key(fitsptr,TLONG,keyname,&val,NULL,&sta))
   printerror(sta);
 return val;
}

/*!
  Read a fitsheader key into string
  \param fitsptr : cfitio pointer to Fits file
  \param keyname : name of the key
  \return value into string
*/
string FitsOpenFile::ReadKeyS(fitsfile *fitsptr,char *keyname)
{
 if(keyname==NULL || fitsptr==NULL) return (string)"";
 int sta=0; char val[FLEN_VALUE];
 if(fits_read_key(fitsptr,TSTRING,keyname,val,NULL,&sta))
   printerror(sta);
 string sval = val;
 return sval;
}

/*!
  CFitsIO error printing routine
  \param sta : cfitio error return code
*/
 void FitsOpenFile::printerror(int sta)
 {
  int stat = sta;
  fits_report_error(stdout,stat);
  fflush(stdout);
  return;
 }

///////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////
//! Class for reading a column in a FITS ASCII or BINARY table

/*!
  \class SOPHYA::FitsABTColRd
  \ingroup FitsIOServer
  Class for reading a column in a FITS ASCII or BINARY table
  \verbatim
  -- Exemple:
  // Open the fits file with FitsOpenFile
  FitsOpenFile fof = new FitsOpenFile("myfits.fits");
  // Select the column to be read
  FitsABTColRd fbt(fof,"BoloMuv_28",0,1000,1,3);
  FitsABTColRd fbt2(fof,"BoloMuv_29",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);
    double y = fbt2.Read(i);
    if(i%lpmod==0) cout<<i<<": "<<x<<", "<<y<<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;
  TVector<double> data2;
  fbt2.Read(32,50,data);
  // Close the fits file
  delete fof;
  \endverbatim
*/

//////////////////////////////////////////////////////////////
/*!
  Constructor.
  \param fof : Pointer to the Class for opening the FITS file
  \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 first BINARY or ASCII table is taken
  - if ihdu>nhdu ihdu is set to nhdu
  - bsens>0    read forward
    bsens<0    read backward
    bsens==0   read centered
  \endverbatim
  \warning ihdu = [1,nhdu]
*/
FitsABTColRd::FitsABTColRd(FitsOpenFile* fof,string collabel
                          ,int ihdu,long blen,long bsens,int lp)
{
  Init(fof,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[
*/
FitsABTColRd::FitsABTColRd(FitsOpenFile* fof,int colnum
                          ,int ihdu,long blen,long bsens,int lp)
{
  Init(fof,"",colnum,ihdu,blen,bsens,lp);
}

/*! Constructor by copy */
FitsABTColRd::FitsABTColRd(FitsABTColRd& fbt)
{
  Init(fbt.GetFitsOpenFile(),fbt.GetColLabel().c_str()
      ,fbt.GetColNum(),fbt.HDU()
      ,fbt.GetBLen(),fbt.GetBSens(),fbt.DbgLevel);
}

/*! Constructor by default */
FitsABTColRd::FitsABTColRd()
{
 ColLabel = ""; ColTUnit = ""; ColTForm = "";
 ColNum = -1; ColTypeCode = 0;
 NBcol = 0; NBline = 0;
 SetNulVal(); SetDebug(0);
 NFitsRead = 0;
 FitsOF = NULL; FitsPtr = NULL;
 LineDeb = LineFin = -1;
 Buffer = NULL;
}

/*! Init routine called by the constructor */
void FitsABTColRd::Init(FitsOpenFile* fof,const char* collabel,int colnum
                        ,int ihdu,long blen,long bsens,int lp)
{
 // Initialisation des Parametres Generaux
 ColLabel=collabel; ColTUnit=""; ColTForm=""; ColNum=colnum; ColTypeCode=0;
 NBcol = 0; NBline = 0;
 SetNulVal(); SetDebug(lp);
 NFitsRead = 0;
 FitsOF = NULL; FitsPtr = NULL;
 LineDeb = LineFin = -1;
 Buffer = NULL;

 // Caracteristiques du FitsOpenFile
 FitsOF = fof;
 if(FitsOF==NULL)
   throw NullPtrError("FitsABTColRd::Init: FitsOpenFile pointer is NULL\n");

 FitsPtr = FitsOF->GetFitsPtr();
 if(FitsPtr==NULL)
   throw NullPtrError("FitsABTColRd::Init: FitsPtr pointer is NULL\n");

 int sta = 0;
 if(ihdu<0) ihdu=0; if(ihdu>NHDU()) ihdu=NHDU();

 // Get HDU for bin/ascii table
 // ATTENTION: le fichier est ouvert mais non positionne sur un HDU,
 // une classe utilisant ce fichier doit le positionner sur un HDU.
 // Par contre, si une autre classe utilise ce meme FitsOpenFile,
 // elle ne peut le positionner que sur ce meme HDU !
 if(FitsOF->GetPosStatus()==false) {
   if(ihdu==0) { // find the first BINARY then the first ASCII
     int rc = FitsOF->MoveToFirst(BINARY_TBL);
     if(rc!=BINARY_TBL) FitsOF->MoveToFirst(ASCII_TBL);
   } else {
     int rc = FitsOF->MoveToHDU(ihdu);
     if(rc!=ihdu)
       throw RangeCheckError("FitsABTColRd::Init: Error moving to requested HDU\n");
   }
 } else { // Fits file has already been positionned
   if(ihdu>0 && ihdu!=HDU())
     throw RangeCheckError("FitsABTColRd::Init: file already posit. at another HDU\n");
 }

 // Check HDUType and set position status to TRUE
 if(HDUType()!=BINARY_TBL && HDUType()!=ASCII_TBL)
   throw TypeMismatchExc("FitsABTColRd::Init: HDU not ASCII/BINARY table\n");
 if(DbgLevel>1) cout<<"...Init ihdu="<<ihdu<<" HduType="<<HDUType()<<endl;
 FitsOF->SetPosStatus(true);

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

 // Get number of rows
 if(fits_get_num_rows(FitsPtr,&NBline,&sta)) {
   FitsOpenFile::printerror(sta);
   throw NotAvailableOperation("FitsABTColRd::Init: Error getting number of rows\n");
 }
 if(DbgLevel>1) cout<<"...Init  NBline="<<NBline<<endl;
 if(NBline<1)
   throw RangeCheckError("FitsABTColRd::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)) {
     FitsOpenFile::printerror(sta);
     throw NotAvailableOperation("FitsABTColRd::Init: Error getting column name\n");
   }
   ColNum--;  // Convention [0,ncol[
 }
 if(DbgLevel>1) cout<<"...Init  ColNum="<<ColNum<<endl;
 if(ColNum<0 || ColNum>=NBcol)
   throw RangeCheckError("FitsABTColRd::Init: Bad column number\n");

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

 // Get column name back, tunit, tform
 char tunit[64], tform[64], tdisp[64];
 long repeat=0; double tscale=1., tzero=0.;
 int rc=0;
 if(HDUType()==BINARY_TBL) {
   fits_get_bcolparms(FitsPtr,ColNum+1,labelcol,tunit,tform
                     ,&repeat,&tscale,&tzero,NULL,tdisp,&sta);
 } else {
   fits_get_acolparms(FitsPtr,ColNum+1,labelcol,&repeat,tunit,tform
                     ,&tscale,&tzero,NULL,tdisp,&sta);
 }
 if(rc) {
   FitsOpenFile::printerror(sta);
   throw RangeCheckError("FitsABTColRd::Init: Error getting the column caracteristics\n");
 }
 ColLabel = labelcol;
 ColTUnit = tunit;
 ColTForm = tform;

 // Set the buffer for reading
 ChangeBuffer(blen,bsens);

 if(DbgLevel)
   cout<<"FitsABTColRd::Init Num="<<ColNum<<" Label="<<ColLabel
       <<" TypeCode="<<ColTypeCode<<" TUnit="<<ColTUnit<<" TForm="<<ColTForm<<endl;
 if(DbgLevel>1)
   cout<<"      (repeat="<<repeat<<",tscale="<<tscale<<",tzero="<<tzero
       <<",tdisp="<<tdisp<<")"<<endl;

}

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

/*! Delete called by the destructor */
void FitsABTColRd::Delete(void)
{
 if(Buffer!=NULL) {delete [] Buffer; Buffer=NULL;}
 LineDeb = LineFin = -1;
 //--- Surtout on ne "fits_close_file" pas le fichier FITS !!!
}

//////////////////////////////////////////////////////////////
/*! Change the buffer caracteristiques (see creator) */
void FitsABTColRd::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;
}

//////////////////////////////////////////////////////////////
/*!
  Read a fitsheader key into double
  \param keyname : name of the key
  \return value into double
*/
double FitsABTColRd::ReadKey(char *keyname)
{
 return FitsOpenFile::ReadKey(FitsPtr,keyname);
}

/*!
  Read a fitsheader key into long
  \param keyname : name of the key
  \return value into long
*/
long FitsABTColRd::ReadKeyL(char *keyname)
{
 return FitsOpenFile::ReadKeyL(FitsPtr,keyname);
}

/*!
  Read a fitsheader key into string
  \param keyname : name of the key
  \return value into string
*/
string FitsABTColRd::ReadKeyS(char *keyname)
{
 return FitsOpenFile::ReadKeyS(FitsPtr,keyname);
}

/////////////////////////////////////////////////
/*!
  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 FitsABTColRd::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("FitsABTColRd::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) {
     FitsOpenFile::printerror(sta);
     throw NotAvailableOperation("FitsABTColRd::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("FitsABTColRd::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) {
     FitsOpenFile::printerror(sta);
     LineDeb = LineFin = -1;
     throw NotAvailableOperation("FitsABTColRd::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 FitsABTColRd::Read(long n1,long n2,TVector<double>& data)
{
 if(n1<0 || n1>=NBline)
   throw RangeCheckError("FitsABTColRd::Read TVector bad requested 1srt line \n");
 if(data.Size()<=0 && n2<n1)
   throw RangeCheckError("FitsABTColRd::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) {
   FitsOpenFile::printerror(sta);
   throw NotAvailableOperation("FitsABTColRd::Read_TVector<double>: Error Reading Fits file\n");
 }

 return nread;
}

/*! idem before but for TVector of float */
long FitsABTColRd::Read(long n1,long n2,TVector<float>& data)
{
 if(n1<0 || n1>=NBline)
   throw RangeCheckError("FitsABTColRd::Read TVector bad requested 1srt line \n");
 if(data.Size()<=0 && n2<n1)
   throw RangeCheckError("FitsABTColRd::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) {
   FitsOpenFile::printerror(sta);
   throw NotAvailableOperation("FitsABTColRd::Read_TVector<float>: Error Reading Fits file\n");
 }

 return nread;
}

/*! idem before but for TVector of unsigned short */
long FitsABTColRd::Read(long n1,long n2,TVector<uint_2>& data)
{
 if(n1<0 || n1>=NBline)
   throw RangeCheckError("FitsABTColRd::Read TVector bad requested 1srt line \n");
 if(data.Size()<=0 && n2<n1)
   throw RangeCheckError("FitsABTColRd::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);

 int sta=0;
 fits_read_col(FitsPtr,TUSHORT,ColNum+1,n1+1,1,nread,NULL,data.Data(),NULL,&sta);
 if(sta) {
   FitsOpenFile::printerror(sta);
   throw NotAvailableOperation("FitsABTColRd::Read_TVector<uint_2>: Error Reading Fits file\n");
 }

 return nread;
}

/*! idem before but for TVector of int_4 */
long FitsABTColRd::Read(long n1,long n2,TVector<int_4>& data)
{
 if(n1<0 || n1>=NBline)
   throw RangeCheckError("FitsABTColRd::Read TVector bad requested 1srt line \n");
 if(data.Size()<=0 && n2<n1)
   throw RangeCheckError("FitsABTColRd::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) {
   FitsOpenFile::printerror(sta);
   throw NotAvailableOperation("FitsABTColRd::Read_TVector<int_4>: Error Reading Fits file\n");
 }

 return nread;
}

/*! idem before but for TVector of int_8 */
long FitsABTColRd::Read(long n1,long n2,TVector<int_8>& data)
{
#ifdef TLONGLONG
 if(n1<0 || n1>=NBline)
   throw RangeCheckError("FitsABTColRd::Read TVector bad requested 1srt line \n");
 if(data.Size()<=0 && n2<n1)
   throw RangeCheckError("FitsABTColRd::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);

 int sta=0;
 fits_read_col(FitsPtr,TLONGLONG,ColNum+1,n1+1,1,nread,NULL,data.Data(),NULL,&sta);
 if(sta) {
   FitsOpenFile::printerror(sta);
   throw NotAvailableOperation("FitsABTColRd::Read_TVector<int_8>: Error Reading Fits file\n");
 }

 return nread;
#else
  throw PException("FitsABTColRd::Read(..,TVector<int_8>&) Not in that cfitsio version");
#endif
}

/////////////////////////////////////////////////
/*!
  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 FitsABTColRd::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 FitsABTColRd::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;
}

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

///////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////

//! Class for reading a column in a FITS ASCII or BINARY table with fits file opening

/*!
  \class SOPHYA::FitsABTColRead
  \ingroup FitsIOServer
  Class for reading a column in a FITS ASCII or BINARY table with fits file opening
  \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 first BINARY or ASCII table is taken
  - if ihdu>nhdu ihdu is set to nhdu
  - 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)
: FitsABTColRd(new FitsOpenFile(fname),collabel,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)
: FitsABTColRd(new FitsOpenFile(fname),colnum,ihdu,blen,bsens,lp)
{
}

/*! Constructor. see below */
FitsABTColRead::FitsABTColRead(const char * cfname,const char* collabel
                              ,int ihdu,long blen,long bsens,int lp)
: FitsABTColRd(new FitsOpenFile(cfname),collabel,ihdu,blen,bsens,lp)
{
}

/*! Constructor. see below */
FitsABTColRead::FitsABTColRead(const char * cfname,int colnum
                              ,int ihdu,long blen,long bsens,int lp)
: FitsABTColRd(new FitsOpenFile(cfname),colnum,ihdu,blen,bsens,lp)
{
}
/*! Constructor by default */
FitsABTColRead::FitsABTColRead()
{
}

/*! Constructor by copy */
FitsABTColRead::FitsABTColRead(FitsABTColRead& fbt)
{
 // --- ATTENTION ---
 // FitsABTColRead ferme le fichier FITS: il faut dupliquer le FitsOpenFile
 FitsOpenFile* fof = new FitsOpenFile(*fbt.GetFitsOpenFile());
 Init(fof,fbt.GetColLabel().c_str()
     ,fbt.GetColNum(),fbt.HDU()
     ,fbt.GetBLen(),fbt.GetBSens(),fbt.DbgLevel);
}

/*! Destructor. */
FitsABTColRead::~FitsABTColRead()
{
 Delete();
 if(FitsOF!=NULL) delete FitsOF;
}

///////////////////////////////////////////////////////////////////
//! Class for reading a 2D image from a FITS file

/*!
  \class SOPHYA::FitsImg2DRd
  \ingroup FitsIOServer
  Class for reading a 2D image from a FITS file
*/

//////////////////////////////////////////////////////////////
/*!
  Constructor.
  \param fof : Pointer to the Class for opening the FITS file
  \param ihdu : number of the HDU where the image is.
  \param lp : debug level
  \verbatim
  - if ihdu<=0 first IMAGE hdu is taken
  - if ihdu>nhdu ihdu is set to nhdu
  \endverbatim
  \warning ihdu = [1,nhdu]
*/
FitsImg2DRd::FitsImg2DRd(FitsOpenFile* fof,int ihdu,int lp)
{
 Init(fof,ihdu,lp);
}

/*! Constructor by copy */
FitsImg2DRd::FitsImg2DRd(FitsImg2DRd& fbt)
{
 Init(fbt.GetFitsOpenFile(),fbt.HDU(),fbt.DbgLevel);
}

/*! Constructor by default */
FitsImg2DRd::FitsImg2DRd()
{
 Naxis[0] = Naxis[1] = 0;
 SetNulVal(); SetDebug(0);
 FitsOF = NULL; FitsPtr = NULL;
}

/*! Init routine called by the constructor */
void FitsImg2DRd::Init(FitsOpenFile* fof,int ihdu,int lp)
{
 // Initialisation des Parametres Generaux
 Naxis[0] = Naxis[1] = 0;
 SetNulVal(); SetDebug(lp);
 FitsOF = NULL; FitsPtr = NULL;

 // Caracteristiques du FitsOpenFile
 FitsOF = fof;
 if(FitsOF==NULL)
   throw NullPtrError("FitsImg2DRd::Init: FitsOpenFile pointer is NULL\n");

 FitsPtr = FitsOF->GetFitsPtr();
 if(FitsPtr==NULL)
   throw NullPtrError("FitsImg2DRd::Init: FitsPtr pointer is NULL\n");

 int sta = 0;
 if(ihdu<0) ihdu=0; if(ihdu>NHDU()) ihdu=NHDU(); 

 // Get HDU 2D image
 // ATTENTION: ... cf blabla equivalent dans FitsABTColRd::Init()
 if(FitsOF->GetPosStatus()==false) {
   if(ihdu==0) { // find the first IMAGE_HDU
     FitsOF->MoveToFirst(IMAGE_HDU);
   } else {
     int rc = FitsOF->MoveToHDU(ihdu);
     if(rc!=ihdu)
       throw RangeCheckError("FitsABTColRd::Init: Error moving to requested HDU\n");
   }
 } else { // Fits file has already been positionned
   if(ihdu>0 && ihdu!=HDU())
     throw RangeCheckError("FitsABTColRd::Init: file already posit. at another HDU\n");
 }

 // Check HDUType and set position status to TRUE
 if(HDUType()!=IMAGE_HDU)
   throw TypeMismatchExc("FitsImg2DRd::Init: HDU not IMAGE_HDU\n");
 FitsOF->SetPosStatus(true);

 // Get NAXIS 1 et 2
 int nfound=0;
 if(fits_read_keys_lng(FitsPtr,"NAXIS",1,2,Naxis,&nfound,&sta)) {
   FitsOpenFile::printerror(sta);
   throw RangeCheckError("FitsImg2DRd::Init: Error reading NAXIS cards\n");
 }
 if(DbgLevel>1)
   cout<<"...Init(hdu="<<HDU()<<")  NAXIS1="<<Naxis[0]<<" NAXIS2="
       <<Naxis[1]<<" (nfound="<<nfound<<")"<<endl;
 if(nfound!=2 || Naxis[0]<=0 || Naxis[1]<=0)
   throw NotAvailableOperation("FitsImg2DRd::Init: bad Naxis[0-1] value\n");

}

/*! Destructor. */
FitsImg2DRd::~FitsImg2DRd()
{
 //--- Surtout on ne "fits_close_file" pas le fichier FITS !!!
 Naxis[0] = Naxis[1] = 0;
}

//////////////////////////////////////////////////////////////
/*!
  Read a fitsheader key into double
  \param keyname : name of the key
  \return value into double
*/
double FitsImg2DRd::ReadKey(char *keyname)
{
 return FitsOpenFile::ReadKey(FitsPtr,keyname);
}

/*!
  Read a fitsheader key into long
  \param keyname : name of the key
  \return value into long
*/
long FitsImg2DRd::ReadKeyL(char *keyname)
{
 return FitsOpenFile::ReadKeyL(FitsPtr,keyname);
}

/*!
  Read a fitsheader key into string
  \param keyname : name of the key
  \return value into string
*/
string FitsImg2DRd::ReadKeyS(char *keyname)
{
 return FitsOpenFile::ReadKeyS(FitsPtr,keyname);
}

//////////////////////////////////////////////////////////////
/* REMARQUE:
 * Si une image FITS a NAXIS1=100 et NAXIS2=50
 * alors un tableau 2D juste assez grand pour contenir l'image
 * doit etre declare array[50][100] (et non pas array[100][50])
 *                   array[NAXIS2][NAXIS1]
 */
/*!
Read image into a TMatrix<uint_2>
\warning TMatrix data(Naxis2,Naxis1)
*/
long FitsImg2DRd::Read(TMatrix<uint_2>& data)
{
 int sta=0;
 uint_2* arr = new uint_2[Naxis[0]];
 data.ReSize(Naxis[1],Naxis[0]);

 for(int j=0;j<Naxis[1];j++) {
   long deb = j*Naxis[0]+1, nel = Naxis[0];
   fits_read_img(FitsPtr,TUSHORT,deb,nel,&NulVal,arr,NULL,&sta);
   if(sta) {
     FitsOpenFile::printerror(sta); delete [] arr;
     throw
       NotAvailableOperation("FitsImg2DRd::Read(TMatrix<uint_2>): Error Reading Fits file\n");
   }
   for(int i=0;i<Naxis[0];i++) data(j,i) = arr[i];
 }

 delete [] arr;
 return Naxis[0]*Naxis[1];
 }

/*! Read image into a TMatrix<int_4> */
long FitsImg2DRd::Read(TMatrix<int_4>& data)
{
 int sta=0;
 int_4* arr = new int_4[Naxis[0]];
 data.ReSize(Naxis[1],Naxis[0]);
 int T = (sizeof(long)==4) ? TLONG: TINT;

 for(int j=0;j<Naxis[1];j++) {
   long deb = j*Naxis[0]+1, nel = Naxis[0];
   fits_read_img(FitsPtr,T,deb,nel,&NulVal,arr,NULL,&sta);
   if(sta) {
     FitsOpenFile::printerror(sta); delete [] arr;
     throw
       NotAvailableOperation("FitsImg2DRd::Read(TMatrix<int_4>): Error Reading Fits file\n");
   }
   for(int i=0;i<Naxis[0];i++) data(j,i) = arr[i];
 }

 delete [] arr;
 return Naxis[0]*Naxis[1];
 }

/*! Read image into a TMatrix<int_8> */
long FitsImg2DRd::Read(TMatrix<int_8>& data)
{
 int sta=0;
 int_8* arr = new int_8[Naxis[0]];
 data.ReSize(Naxis[1],Naxis[0]);

 for(int j=0;j<Naxis[1];j++) {
   long deb = j*Naxis[0]+1, nel = Naxis[0];
   fits_read_img(FitsPtr,TLONGLONG,deb,nel,&NulVal,arr,NULL,&sta);
   if(sta) {
     FitsOpenFile::printerror(sta); delete [] arr;
     throw
       NotAvailableOperation("FitsImg2DRd::Read(TMatrix<int_8>): Error Reading Fits file\n");
   }
   for(int i=0;i<Naxis[0];i++) data(j,i) = arr[i];
 }

 delete [] arr;
 return Naxis[0]*Naxis[1];
 }

/*! Read image into a TMatrix<float> */
long FitsImg2DRd::Read(TMatrix<float>& data)
{
 int sta=0;
 float* arr = new float[Naxis[0]];
 data.ReSize(Naxis[1],Naxis[0]);

 for(int j=0;j<Naxis[1];j++) {
   long deb = j*Naxis[0]+1, nel = Naxis[0];
   fits_read_img(FitsPtr,TFLOAT,deb,nel,&NulVal,arr,NULL,&sta);
   if(sta) {
     FitsOpenFile::printerror(sta); delete [] arr;
     throw
       NotAvailableOperation("FitsImg2DRd::Read(TMatrix<float>): Error Reading Fits file\n");
   }
   for(int i=0;i<Naxis[0];i++) data(j,i) = arr[i];
 }

 delete [] arr;
 return Naxis[0]*Naxis[1];
 }

/*! Read image into a TMatrix<double> */
long FitsImg2DRd::Read(TMatrix<double>& data)
{
 int sta=0;
 double* arr = new double[Naxis[0]];
 data.ReSize(Naxis[1],Naxis[0]);

 for(int j=0;j<Naxis[1];j++) {
   long deb = j*Naxis[0]+1, nel = Naxis[0];
   fits_read_img(FitsPtr,TDOUBLE,deb,nel,&NulVal,arr,NULL,&sta);
   if(sta) {
     FitsOpenFile::printerror(sta); delete [] arr;
     throw
       NotAvailableOperation("FitsImg2DRd::Read(TMatrix<double>): Error Reading Fits file\n");
   }
   for(int i=0;i<Naxis[0];i++) data(j,i) = arr[i];
 }

 delete [] arr;
 return Naxis[0]*Naxis[1];
 }