/* Lecteur de colonne de table Fits (binaire ou ASCII) avec buffer */
#include "sopnamsp.h"
#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.
  You can read many columns of the same FITS table by instanciating
  many FitsABTColRd on the same FitsOpenFile. So, the FITS file is
  opened only once. Of course the various FitsABTColRd must read
  the same FITS file HDU.
  \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;
 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;
 LineDeb = LineFin = -1;
 Buffer = NULL;

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

 if(GetFitsPtr()==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(GetFitsPtr(),&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(GetFitsPtr(),&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(GetFitsPtr(),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(GetFitsPtr(),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(GetFitsPtr(),ColNum+1,labelcol,tunit,tform
                     ,&repeat,&tscale,&tzero,NULL,tdisp,&sta);
 } else {
   fits_get_acolparms(GetFitsPtr(),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(GetFitsPtr(),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(GetFitsPtr(),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(GetFitsPtr(),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(GetFitsPtr(),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(GetFitsPtr(),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(GetFitsPtr(),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(GetFitsPtr(),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(GetFitsPtr(),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(GetFitsPtr(),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(GetFitsPtr(),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.
  The FITS file is opened each time you instanciate a FitsABTColRead.
  So reading "n" columns of the same FITS table by instanciating "n"
  FitsABTColRead, will open "n" times te FITS file.
  Use FitsABTColRd if you want to open the FITS file only once.
  \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()
: FitsABTColRd()
{
}

/*! 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(); // ?? inutile ??
 // On detruit le FitsOpenFile, cad qu'on ferme (fits_file_close) le fichier FITS
 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;
}

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

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

 if(GetFitsPtr()==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("FitsImg2DRd::Init: Error moving to requested HDU\n");
   }
 } else { // Fits file has already been positionned
   if(ihdu>0 && ihdu!=HDU())
     throw RangeCheckError("FitsImg2DRd::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(GetFitsPtr(),"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(GetFitsPtr(),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(GetFitsPtr(),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(GetFitsPtr(),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)
*/
size_t 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(GetFitsPtr(),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> */
size_t 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(GetFitsPtr(),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> */
size_t 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(GetFitsPtr(),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> */
size_t 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(GetFitsPtr(),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> */
size_t 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(GetFitsPtr(),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];
 }

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

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

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

//////////////////////////////////////////////////////////////
/*!
  Constructor.
  \param fname : name of 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]
*/
FitsImg2DRead::FitsImg2DRead(string fname,int ihdu,int lp)
: FitsImg2DRd(new FitsOpenFile(fname),ihdu,lp)
{
}

/*! Constructor. see below */
FitsImg2DRead::FitsImg2DRead(const char * cfname,int ihdu,int lp)
: FitsImg2DRd(new FitsOpenFile(cfname),ihdu,lp)
{
}

/*! Constructor by default */
FitsImg2DRead::FitsImg2DRead()
: FitsImg2DRd()
{
}

/*! Constructor by copy */
FitsImg2DRead::FitsImg2DRead(FitsImg2DRead& fimg)
{
 // --- ATTENTION ---
 // FitsImg2DRead ferme le fichier FITS: il faut dupliquer le FitsOpenFile
 FitsOpenFile* fof = new FitsOpenFile(*fimg.GetFitsOpenFile());
 Init(fof,fimg.HDU(),fimg.DbgLevel);
}

/*! Destructor. */
FitsImg2DRead::~FitsImg2DRead()
{
 // On detruit le FitsOpenFile, cad qu'on ferme (fits_file_close) le fichier FITS
 if(FitsOF!=NULL) delete FitsOF;
}



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

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

/*!
  \class SOPHYA::FitsImg3DRd
  \ingroup FitsIOServer
  Class for reading a 3D 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 3D 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]
*/
FitsImg3DRd::FitsImg3DRd(FitsOpenFile* fof,int ihdu,int lp)
{
 Init(fof,ihdu,lp);
}

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

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

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

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

 if(GetFitsPtr()==NULL)
   throw NullPtrError("FitsImg3DRd::Init: FitsPtr pointer is NULL\n");

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

 // Get HDU 3D 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("FitsImg3DRd::Init: Error moving to requested HDU\n");
   }
 } else { // Fits file has already been positionned
   if(ihdu>0 && ihdu!=HDU())
     throw RangeCheckError("FitsImg3DRd::Init: file already posit. at another HDU\n");
 }

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

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

}

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

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

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

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

//////////////////////////////////////////////////////////////
/* REMARQUE:
 * Dans TArray A(naxis1,naxis2,naxis3);
 *      A(i,j,k) -> i varie le plus vite et k le moins vite
 */
/*!
Read 3D image into a TArray<uint_2>
*/
size_t FitsImg3DRd::Read(TArray<uint_2>& data)
{
 int sta=0;
 uint_2* arr = new uint_2[Naxis[0]];
 sa_size_t ndim[3] = {Naxis[0],Naxis[1],Naxis[2]}; data.ReSize(3,ndim);

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

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

/*! Read 3D image into a TArray<int_4> */
size_t FitsImg3DRd::Read(TArray<int_4>& data)
{
 int sta=0;
 int_4* arr = new int_4[Naxis[0]];
 sa_size_t ndim[3] = {Naxis[0],Naxis[1],Naxis[2]}; data.ReSize(3,ndim);
 int T = (sizeof(long)==4) ? TLONG: TINT;

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

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

/*! Read 3D image into a TArray<int_8> */
size_t FitsImg3DRd::Read(TArray<int_8>& data)
{
 int sta=0;
 int_8* arr = new int_8[Naxis[0]];
 sa_size_t ndim[3] = {Naxis[0],Naxis[1],Naxis[2]}; data.ReSize(3,ndim);

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

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

/*! Read 3D image into a TArray<float> */
size_t FitsImg3DRd::Read(TArray<float>& data)
{
 int sta=0;
 float* arr = new float[Naxis[0]];
 sa_size_t ndim[3] = {Naxis[0],Naxis[1],Naxis[2]}; data.ReSize(3,ndim);

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

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

/*! Read 3D image into a TArray<double> */
size_t FitsImg3DRd::Read(TArray<double>& data)
{
 int sta=0;
 double* arr = new double[Naxis[0]];
 sa_size_t ndim[3] = {Naxis[0],Naxis[1],Naxis[2]}; data.ReSize(3,ndim);

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

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

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

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

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

//////////////////////////////////////////////////////////////
/*!
  Constructor.
  \param fname : name of the FITS file
  \param ihdu : number of the HDU where the 3D 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]
*/
FitsImg3DRead::FitsImg3DRead(string fname,int ihdu,int lp)
: FitsImg3DRd(new FitsOpenFile(fname),ihdu,lp)
{
}

/*! Constructor. see below */
FitsImg3DRead::FitsImg3DRead(const char * cfname,int ihdu,int lp)
: FitsImg3DRd(new FitsOpenFile(cfname),ihdu,lp)
{
}

/*! Constructor by default */
FitsImg3DRead::FitsImg3DRead()
: FitsImg3DRd()
{
}

/*! Constructor by copy */
FitsImg3DRead::FitsImg3DRead(FitsImg3DRead& fimg)
{
 // --- ATTENTION ---
 // FitsImg3DRead ferme le fichier FITS: il faut dupliquer le FitsOpenFile
 FitsOpenFile* fof = new FitsOpenFile(*fimg.GetFitsOpenFile());
 Init(fof,fimg.HDU(),fimg.DbgLevel);
}

/*! Destructor. */
FitsImg3DRead::~FitsImg3DRead()
{
 // On detruit le FitsOpenFile, cad qu'on ferme (fits_file_close) le fichier FITS
 if(FitsOF!=NULL) delete FitsOF;
}