source: Sophya/trunk/SophyaExt/FitsIOServer/fitsfile.cc@ 1193

#include "machdefs.h"
#include <stdlib.h>
#include "fitsfile.h"
#include "pexceptions.h"
#include "strutil.h"
#include "anydataobj.h" 
#include "fitsspherehealpix.h"


void   FitsIOHandler::Read(char flnm[],int hdunum)
{
  FitsInFile ifts(flnm);
  Read(ifts, hdunum);
}
void FitsIOHandler::Read(FitsInFile& is, int hdunum) 
{
    if (hdunum == 0)  is.moveToFollowingHeader();
   else is.ReadFInit(hdunum);
  ReadFromFits(is);
}


void FitsIOHandler::Write(char flnm[]) 

{
  FitsOutFile of(flnm, unknown);
  Write(of);
}

void FitsIOHandler::Write(FitsOutFile& os) 
{
  WriteToFits(os);
}



 
FitsFile::~FitsFile()
{
  int status = 0;
    if( fptr_ != NULL)
    {
      fits_close_file(fptr_,&status);
      // je ne fais pas      delete fptr_, c'est la lib. fitsio qui a fait 
      // new...
    }
   if( status )  printerror( status );
}


void FitsFile::printerror(int &status) 
  //*****************************************************/
  //* Print out cfitsio error messages and exit program */
  //*****************************************************/
{
  if( status )
    {
      fits_report_error(stderr,status);
      throw IOExc("FitsFile:: error FITSIO status");
    }
  return;
}

void FitsFile::printerror(int& status, char* texte) 
  //*****************************************************/
  //* Print out cfitsio error messages and exit program */
  //*****************************************************/
{
  // print out cfitsio error messages and exit program 
  // print error report  
  fits_report_error(stderr, status);
  cout << " erreur:: " << texte << endl;
  throw IOExc("FitsFile:: error FITSIO status");
}

void FitsFile::ResetStatus(int& status) 
{
  fits_status_ = status;
  status = 0;
}

string FitsFile::getErrStatus(int status) 
{
  char text[31];
  fits_get_errstatus(status, text);
  return string(text);
}

FitsInFile::FitsInFile()
{
  InitNull();
}
//FitsInFile::FitsInFile(char flnm[], int hdunum)
FitsInFile::FitsInFile(char flnm[])
{
    InitNull();
   int status = 0;
   fits_open_file(&fptr_,flnm,READONLY,&status);
   if( status ) printerror( status );
}


void FitsInFile::InitNull()
{
  
  bitpix_ = 0;
  naxis_ = 0;
  nbData_ = 0;
  nrows_ = 0;
  nbcols_ = 0;
  naxisn_.clear();
  repeat_.clear();
  noms_.clear();
  taille_des_chaines_.clear();
  dvl_.Clear();
  
  
}


int FitsInFile::NbBlocks(char flnm[])
{
  int status = 0;
  int nbhdu = 0; 
  fitsfile* fileptr; 
  fits_open_file(&fileptr,flnm,READONLY,&status);
  if( status ) printerror( status, "NbBlocks: erreur ouverture fichier" );
  fits_get_num_hdus(fileptr, &nbhdu, &status); 
  fits_close_file(fileptr,&status);
  return nbhdu;
}

void FitsInFile::getBlockType(char flnm[], int hdunum, string& typeOfExtension, int& naxis, vector<int>& naxisn, string& dataType, DVList& dvl )
{
  int status = 0;
  fitsfile* fileptr; 
  fits_open_file(&fileptr,flnm,READONLY,&status);
  if( status ) printerror( status, "getBlockType: erreur ouverture fichier" );
  // move to the specified HDU number
  int hdutype = 0;
  fits_movabs_hdu(fileptr,hdunum,&hdutype,&status);
  if( status ) printerror( status,"getBlockType: erreur movabs");
  if(hdutype == IMAGE_HDU) 
    {
      typeOfExtension = "IMAGE";
      int bitpix;
      GetImageParameters (fileptr, bitpix, naxis, naxisn);
      if(bitpix == DOUBLE_IMG) dataType = "double";
      else
        if(bitpix == FLOAT_IMG) dataType = "float";
        else
          if(bitpix == LONG_IMG || bitpix == SHORT_IMG ) dataType = "int";
          else 
            {
              cout << " bitpix= " << bitpix << endl;
              throw PException(" FitsFile::getBlockType : unsupprted FITS data type");
            }
      
    }
  else
    if(hdutype == ASCII_TBL  || hdutype == BINARY_TBL) 
      {
        int nrows = 0;
        vector<string> noms; 
        vector<char> types;   
        vector<int>  taille_des_chaines;
        GetBinTabParameters(fileptr, naxis, nrows, naxisn, noms, types, taille_des_chaines); 
        int k;
        for (k=0; k< naxisn.size(); k++) naxisn[k] *= nrows;
        if(hdutype == ASCII_TBL)
          {
            typeOfExtension = "ASCII_TBL"; 
            dataType = "ASCII";
          }
        else
          {
            typeOfExtension = "BINARY_TBL";
            if(types[0] == 'D') dataType = "double";
            else
              if(types[0] == 'E') dataType = "float";
              else
                if(types[0] == 'I' ) dataType = "int";
                else
                  if(types[0] == 'S' ) dataType = "char*";
                  else 
                    {
                      cout << " types[0]= " << types[0] << endl;
                      throw PException(" FitsFile::getBlockType : unsupprted FITS data type");
            }
          }
      }
    else
      {
        cout << " hdutype= " << hdutype << endl;
        throw IOExc("FitsFile::getBlockType: this HDU type is unknown");
      }
  
  KeywordsIntoDVList(fileptr, dvl, hdunum);
  fits_close_file(fileptr,&status);
}

void FitsInFile::GetImageParameters (fitsfile* fileptr,int& bitpix,int& naxis,vector<int>& naxisn)
{
  int hdunum=0;
  cout << " Reading a FITS image in HDU : " << fits_get_hdu_num(fileptr,&hdunum)  << endl;
  int status= 0;

  // bits per pixels
  fits_read_key(fileptr,TINT,"BITPIX",&bitpix,NULL,&status);
  if( status )  printerror( status );

  // number of dimensions in the FITS array
  naxis= 0;
  fits_read_key(fileptr,TINT,"NAXIS",&naxis,NULL,&status);
  if( status ) printerror( status );
  // read the NAXISn keywords to get image size
  long* naxes = new long[naxis] ;
  int nfound;
  fits_read_keys_lng(fileptr,"NAXIS",1,naxis,naxes,&nfound,&status);
  if( status ) printerror( status );
  if (nfound != naxis ) 
    cout << " WARNING : " << nfound << " axes found, expected naxis= " << naxis << endl;
  int k;
  for (k=0; k<naxis; k++) 
    {
      naxisn.push_back( (int)naxes[k] );
    }
  delete [] naxes;
}


void FitsInFile::ReadFInit(int hdunum) 
{
   InitNull();
  // int status = 0;
  
  //  fits_open_file(&fptr_,flnm,READONLY,&status);
  //  if( status ) printerror( status );

  int status = 0;

  if (hdunum <= 1) 
    {
      hdunum_ = 1; 
      // presence of image ?
      int naxis= 0;
      fits_read_key(fptr_,TINT,"NAXIS",&naxis,NULL,&status);
      if( status ) printerror( status );
      if (naxis > 0 )       // there is an image
        {
          hdutype_ = IMAGE_HDU;
          GetImageParameters (fptr_, bitpix_, naxis_, naxisn_);
          nbData_ =  1;
          int k;
          for (k=0; k<naxis_; k++) if (naxisn_[k] > 0) nbData_ *= naxisn_[k];
          KeywordsIntoDVList(fptr_, dvl_,hdunum_);
        }
      else
        {
          throw PException(" first header : no image, probably error in hdunum");
        }
  //
    }
  else 
    {
      hdunum_ = hdunum-1; 
      int hdutype;
      fits_movabs_hdu(fptr_,hdunum_,&hdutype,&status);
      if( status ) printerror( status,":FitsFile::ReadF : erreur movabs");
      moveToFollowingHeader();
    }
}

DVList  FitsInFile::DVListFromPrimaryHeader() const
   {
     int status;
     DVList dvl;
     KeywordsIntoDVList(fptr_, dvl, 1);
     int hdutype = 0;
     if (hdunum_ > 0) fits_movabs_hdu(fptr_,hdunum_,&hdutype,&status);
     return dvl;
   }


void FitsInFile::moveToFollowingHeader()
{
  int status = 0;
  int hdutype;
  InitNull();
  fits_movrel_hdu(fptr_, 1,&hdutype,&status); 
  if( status ) printerror( status," lecture du header suivant" );
  hdunum_++;
  hdutype_= hdutype;
  if(hdutype_ == IMAGE_HDU) 
    {
      GetImageParameters (fptr_, bitpix_, naxis_, naxisn_);
      nbData_ =  1;
      int k;
      for (k=0; k<naxis_; k++) if (naxisn_[k] > 0) nbData_ *= naxisn_[k];
      KeywordsIntoDVList(fptr_, dvl_,hdunum_);
    }
  if(hdutype_ == ASCII_TBL || hdutype_ == BINARY_TBL)
    {
      GetBinTabParameters(fptr_,nbcols_, nrows_,repeat_, noms_, types_, taille_des_chaines_);
      KeywordsIntoDVList(fptr_, dvl_, hdunum_);
    }
}


int  FitsInFile::NbColsFromFits() const
{
  if(hdutype_ ==  BINARY_TBL) return nbcols_;
  else
    if(hdutype_ == ASCII_TBL ||  hdutype_ == IMAGE_HDU) return 1;
    else
      {
        cout << " hdutype= " << hdutype_ << endl;
        throw PException("FitsFile::NbColsFromFits, this HDU is unknown");
      }
}

int FitsInFile::NentriesFromFits(int nocol) const
{
  if(hdutype_ == BINARY_TBL  )   return nrows_*repeat_[nocol];
  else
    if(hdutype_ == ASCII_TBL) return nrows_;
    else
      if(hdutype_ == IMAGE_HDU) return nbData_;
      else
        {
          cout << "hdutype= " << hdutype_ << endl;
          throw PException("FitsFile::NentriesFromFits, this HDU is unknown");
        }
}

char FitsInFile::ColTypeFromFits(int nocol) const
{
  if(hdutype_ != ASCII_TBL && hdutype_ != BINARY_TBL)
    {
      throw IOExc("FitsFile::TypeFromFits, this HDU is not an ASCII table nor a binary table");
    }
  return types_[nocol];
}
string FitsInFile::ColNameFromFits(int nocol) const
{
  if(hdutype_ != ASCII_TBL && hdutype_ != BINARY_TBL)
    {
      throw IOExc("FitsFile::TypeFromFits, this HDU is not an ASCII table nor a binary table");
    }
  return noms_[nocol];
}

int FitsInFile::ColStringLengthFromFits(int nocol) const
{
  if(hdutype_ != ASCII_TBL && hdutype_ != BINARY_TBL)
    {
      throw IOExc("FitsFile::TypeFromFits, this HDU is not an ASCII table nor a binary table");
    }
  int index=-1;
  int k;
  for (k=0; k<=nocol; k++) 
    {
      if (types_[k] == 'S') index++;
    }
  return  taille_des_chaines_[index];
}
void  FitsInFile::GetBinTabLine(int NoLine, double* ddata, float* fdata, int* idata, char ** cdata) 
{
  int status= 0;
  int    anull;
  double dnull= 0.;
  float fnull= 0.;
  int inull= 0;
  char* cnull= "";
  int dcount = 0.;
  int fcount = 0.;
  int icount = 0;
  int ccount =0;
  int ncol;
  long nels=1;
  for (ncol=0; ncol<nbcols_; ncol++)
    {
      switch (types_[ncol])
        { 
        case 'D' :
          fits_read_col(fptr_,TDOUBLE,ncol+1,NoLine+1,1,1,&dnull,&ddata[dcount++],&anull,&status);
            break;
        case 'E' :
          fits_read_col(fptr_,TFLOAT,ncol+1,NoLine+1,1,1,&fnull,&fdata[fcount++],&anull,&status);
          break;
        case 'I' :
          fits_read_col(fptr_,TINT,ncol+1,NoLine+1,1,1,&inull,&idata[icount++],
                        &anull,&status);
          break;
        case 'S' :
          fits_read_col(fptr_,TSTRING,ncol+1,NoLine+1,1,1,cnull,&cdata[ccount++],&anull,&status);
          break;
        }
      if (status)
        {
          ResetStatus(status);
          break;
        }
    }
}

void   FitsInFile::GetBinTabLine(long NoLine,  BnTblLine& ligne) 
{
  int status= 0;
  int    anull;
  double dnull= 0.;
  float fnull= 0.;
  int inull= 0;
  char* cnull= "";
  int dcount = 0.;
  int fcount = 0.;
  int icount = 0;
  int ccount =0;
  int ncol;
  long nels=1;
  for (ncol=0; ncol<nbcols_; ncol++)
    {
      switch (types_[ncol])
        { 
        case 'D' :
          fits_read_col(fptr_,TDOUBLE,ncol+1,NoLine+1,1,1,&dnull,&ligne.ddata_[dcount++],&anull,&status);
            break;
        case 'E' :
          fits_read_col(fptr_,TFLOAT,ncol+1,NoLine+1,1,1,&fnull,&ligne.fdata_[fcount++],&anull,&status);
          break;
        case 'I' :
          fits_read_col(fptr_,TINT,ncol+1,NoLine+1,1,1,&inull,&ligne.idata_[icount++],
                        &anull,&status);
          break;
        case 'S' :
          char* chaine = new char[taille_des_chaines_[ccount]];
          fits_read_col(fptr_,TSTRING,ncol+1,NoLine+1,1,1,cnull,&chaine,&anull,&status);
          ligne.cdata_[ccount++] = string(chaine);
          break;
        }
      if (status)
        {
          ResetStatus(status);
          break;
        }
    }
}

void FitsInFile::GetBinTabLine(int NoLine, float* fdata) 
{
  int status= 0;
  int    anull;
  float fnull= 0.;
  long nels=1;
  int ncol;
  for (ncol=0; ncol<nbcols_; ncol++)
    {
      fits_read_col(fptr_,TFLOAT,ncol+1,NoLine+1,1,1,&fnull,&fdata[ncol],&anull,&status);
      if (status)
        {
          ResetStatus(status);
          break;
        }
    }
}


void FitsInFile::GetBinTabFCol(double* valeurs,int nentries, int NoCol) const
    {
      int status= 0;
      int DTYPE;
      long repeat,width;
            fits_get_coltype(fptr_, NoCol+1,&DTYPE,&repeat,&width,&status);
      if( DTYPE != TDOUBLE) 
        {
          if (DTYPE == TFLOAT)  cout << " WARNING: reading double from float : conversion will be made by fitsio library" << endl;
          else
            throw IOExc("FitsFile::GetBinTabFCol, tentative de lecture non double");
        }
      long nels=nentries;
      int    anull;
      // no checking for undefined pixels
      double dnull= 0.;
      //     fits_read_key(fptr_,TDOUBLE,"BAD_DATA",&dnull,NULL,&status);
      //  if (status != 0) 
      //        {
      //          dnull = -1.6375e30;  // default value
      //          status = 0;
      //        }
      if (nentries != nrows_*repeat)
        {
          cout << " found " << nentries << " pixels, expected:  " << nrows_*repeat << endl;
          throw PException(" FitsFile:::GetBinTabFCol ");
        }
      fits_read_col(fptr_,TDOUBLE,NoCol+1,1,1,nels,&dnull,valeurs,
                            &anull,&status);
      if( status ) printerror( status,"erreur lecture de colonne" );
      
    }

void FitsInFile::GetBinTabFCol(float* valeurs,int nentries, int NoCol) const
    {
      int status= 0;
      int DTYPE;
      long repeat,width;
      fits_get_coltype(fptr_, NoCol+1,&DTYPE,&repeat,&width,&status);
      if( DTYPE != TFLOAT) 
        {
          if (DTYPE == TDOUBLE)  cout << " WARNING: reading float from double : conversion will be made by fitsio library" << endl;
          else
            throw IOExc("FitsFile::GetBinTabFCol, tentative de lecture non float");
        }
      long nels=nentries;
      int    anull;
      // no checking for undefined pixels
      float fnull= 0.;
      //     fits_read_key(fptr_,TFLOAT,"BAD_DATA",&fnull,NULL,&status);
      //    if (status != 0) 
      //        {
      //          fnull = -1.6375e30;  // default value
      //          status = 0;
      //        }
      if (nentries != nrows_*repeat)
        {
          cout << " found " << nentries << " pixels, expected:  " << nrows_*repeat << endl;
          throw PException(" FitsFile:::GetBinTabFCol ");
        }
      fits_read_col(fptr_,TFLOAT,NoCol+1,1,1,nels,&fnull,valeurs,
                            &anull,&status);
      if( status ) printerror( status,"erreur lecture de colonne" );
    }

void FitsInFile::GetBinTabFCol(int* valeurs,int nentries, int NoCol) const
    {
      int status= 0;
      int DTYPE;
      long repeat,width;
      fits_get_coltype(fptr_, NoCol+1,&DTYPE,&repeat,&width,&status);
      if( DTYPE != TLONG &&  DTYPE != TINT &&  DTYPE != TSHORT ) 
        {
          throw IOExc("FitsFile::GetBinTabFCol, tentative de lecture non entier");
        }
      long nels=nentries;
      // no checking for undefined pixels
      int    anull;
      int inull= 0;
      //     fits_read_key(fptr_,TINT,"BAD_DATA",&inull,NULL,&status);
      //     if (status != 0) 
      //        {
      //          inull = -999999;  // default value
      //          status = 0;
      //        }
      if (nentries != nrows_*repeat)
        {
          cout << " found " << nentries << " pixels, expected:  " << nrows_*repeat << endl;
          throw PException(" FitsFile:::GetBinTabFCol ");
        }
      fits_read_col(fptr_,TINT,NoCol+1,1,1,nels,&inull,valeurs,
                            &anull,&status);
      if( status ) printerror( status,"erreur lecture de colonne" );
    }

void FitsInFile::GetBinTabFCol(char** valeurs, int nentries, int NoCol) const
    {
      int status= 0;
      int DTYPE;
      long repeat,width;
      fits_get_coltype(fptr_, NoCol+1,&DTYPE,&repeat,&width,&status);
      if( DTYPE != TSTRING ) 
        {
          throw IOExc("FitsFile::GetBinTabFCol, tentative de lecture non float");
        }
      long nels=nentries;
      // no checking for undefined pixels
      int    anull;
      char* cnull= "";
      if (nentries != nrows_*repeat/width)
        {
          cout << " found " << nentries << " pixels, expected:  " << nrows_*repeat/width << endl;
          throw PException(" FitsFile:::GetBinTabFCol ");
        }
      long frow=1;
      long felem=1;
      fits_read_col(fptr_,TSTRING,NoCol+1,frow,felem,nels,cnull,valeurs,
                    &anull,&status);
      if( status ) printerror( status,"erreur lecture de colonne" );
    }

void FitsInFile::GetSingleColumn(double* map, int nentries) const
{ 
  int status = 0;
  if(hdutype_ == IMAGE_HDU) 
    {
  
      if(bitpix_ != DOUBLE_IMG) 
        {
          cout << " The data type on fits file is not double...";
          cout << " Conversion to double achieved by cfitsio lib" << endl;
        }
  
      // no checking for undefined pixels
      int    anull;
      double dnull= 0.;
  
      long nels= nentries;
      fits_read_img(fptr_,TDOUBLE,1,nels,&dnull,map,&anull,&status);
      if( status ) printerror( status );
    }
  else
    if(hdutype_ == ASCII_TBL || hdutype_ == BINARY_TBL)
      {
        GetBinTabFCol(map,nentries, 0);
      }
    else
      {
      cout << " hdutype= " << hdutype_ << endl;
      throw IOExc("FitsFile::GetSingleColumn this HDU is unknown");
      }
}

void FitsInFile::GetSingleColumn(float* map, int nentries) const
{
  int status = 0;
  if(hdutype_ == IMAGE_HDU) 
    {
      if(bitpix_ != FLOAT_IMG) 
        {
          cout << " The data type on fits file is not float ";
          cout << " Conversion to float achieved by cfitsio lib" << endl;
        }
      // no checking for undefined pixels
      int   anull;
      float fnull= 0.;
      
      long nels= nentries;
      fits_read_img(fptr_,TFLOAT,1,nels,&fnull, map,&anull,&status);
      if( status ) printerror( status );
    }
  else
    if(hdutype_ == ASCII_TBL || hdutype_ == BINARY_TBL)
      {
        GetBinTabFCol(map,nentries, 0);
      }
    else
      {
      cout << " hdutype= " << hdutype_ << endl;
      throw IOExc("FitsFile::GetSingleColumn this HDU is unknown");
      }
}

void FitsInFile::GetSingleColumn( int* map, int nentries) const
{
  int status = 0;
  if(hdutype_ == IMAGE_HDU) 
    {
      if(bitpix_ != LONG_IMG) 
        {
          cout << " The data type on fits file is not int ";
          cout << " Conversion to float achieved by cfitsio lib" << endl;
        }
      // no checking for undefined pixels
      int   anull;
      float fnull= 0.;
      
      long nels= nentries;
      fits_read_img(fptr_,TINT,1,nels,&fnull,map,&anull,&status);
      if( status ) printerror( status );
    }
  else
    if(hdutype_ == ASCII_TBL || hdutype_ == BINARY_TBL)
      {
        GetBinTabFCol(map,nentries, 0);
      }
    else
      {
      cout << " hdutype= " << hdutype_ << endl;
      throw IOExc("FitsFile::GetSingleColumn this HDU is unknown");
      }
}

void FitsInFile::GetBinTabParameters(fitsfile* fileptr, int& nbcols, int& nrows,
                                  vector<int>& repeat,
                                  vector<string>& noms, 
                                  vector<char>& types,   
                                  vector<int>&  taille_des_chaines)
{
  int status= 0;
  int hdunum=0;
  int hdutype=0;
  fits_get_hdu_num(fileptr,&hdunum);
  fits_get_hdu_type(fileptr, &hdutype, &status);

  if(hdutype != ASCII_TBL && hdutype != BINARY_TBL)
    {
      throw IOExc("FitsFile::GetBinTabParameters this HDU is not an ASCII table nor a binary table");
    }
  if(hdutype == ASCII_TBL)
    cout << " Reading a FITS ascii table in HDU : " << hdunum << endl;  
  if(hdutype == BINARY_TBL)
    cout << " Reading a FITS binary table in HDU : " << hdunum << endl;

  // get the number of columns
  fits_get_num_cols(fileptr, &nbcols,&status);
  if( status ) printerror( status );

  // get the number of rows
  long naxis2= 0;
  fits_get_num_rows(fileptr,&naxis2,&status);
  if( status ) printerror( status );
  nrows = (int)naxis2;

  // get the datatype,  names and the repeat count
  noms.clear();
  noms.reserve(nbcols);
  types.clear();
  types.reserve(nbcols);
  repeat.clear();
  repeat.reserve(nbcols);
  taille_des_chaines.clear();
  char **ttype = new char*[nbcols];
  int ii;
  // 
  //
  for (ii=0; ii < nbcols; ii++) ttype[ii]=new char[FLEN_VALUE];
  int nfound;
  fits_read_keys_str(fileptr, "TTYPE",1,nbcols,ttype,&nfound, &status);
  if( status ) printerror( status,"erreur lecture des noms de colonne");
  int rept=0;
  for(ii = 0; ii <  nbcols; ii++)
    {
      int DTYPE;
      long width;
      long repete = 0;
      fits_get_coltype(fileptr,ii+1,&DTYPE,&repete,&width,&status);
      if( status ) printerror( status,"erreur lecture type de colonne");
      rept = repete;
      noms.push_back(string(ttype[ii]));
      switch (DTYPE)
        {
        case TDOUBLE :
          types.push_back('D');
          break;
        case TFLOAT :
          types.push_back('E');
          break;
        case TLONG :
          types.push_back('I');
          break;
        case TINT :
          types.push_back('I');
          break;
        case TSHORT :
          types.push_back('I');
          break;
        case TSTRING :
          types.push_back('S');
          taille_des_chaines.push_back(width);
          rept/=width;
          break;
        default :
          cout << " field " << ii+1 << " DTYPE= " << DTYPE << endl;
          throw IOExc("FitsFile:: unknown type of field");
        }
      repeat.push_back(rept);
    }
  for (ii=0; ii < nbcols; ii++) delete [] ttype[ii];
  delete [] ttype;
}

void FitsInFile::KeywordsIntoDVList(fitsfile* fileptr, DVList& dvl, int hdunum)
{
  int status = 0;
  int hdutype;
  fits_movabs_hdu(fileptr,hdunum,&hdutype,&status);
  if( status ) printerror( status,":KeywordsIntoDVList : erreur movabs");
  // get number of keywords
  int nkeys,keypos;
  fits_get_hdrpos(fileptr,&nkeys,&keypos,&status);
  if( status ) printerror( status );

  // put  keywords in a DVList object
  char keyname[LEN_KEYWORD]= "";
  char strval[FLEN_VALUE]= "";
  char dtype;
  char card[FLEN_CARD];  
  char *comkey = "COMMENT";
  char comment[FLEN_COMMENT];

  // shift with the number of mandatory keywords
  //  int num= 8;
    int num= 0;
    // primary header
    if (hdunum == 1) 
      {
        // read NAXIS
        int naxis=0;
        fits_read_key(fileptr,TINT,"NAXIS",&naxis,NULL,&status);
        // number of mandatory keywords
        num = naxis+3;
      }
    // extensions
    else
      {
        if (hdutype == IMAGE_HDU)
          {
            // read NAXIS
            int naxis=0;
            fits_read_key(fileptr,TINT,"NAXIS",&naxis,NULL,&status);
            // number of mandatory keywords
            num = naxis+5;
          }
        else
          if(hdutype == ASCII_TBL || hdutype == BINARY_TBL)
            {
              // number of mandatory keywords
              num = 8;
            }
      }
  int j;
  for(j = num+1; j <= nkeys; j++)
    {
      fits_read_keyn(fileptr,j,card,strval,NULL,&status);
      if(status) printerror(status);

      strncpy(keyname,card,LEN_KEYWORD-1);
      if(strncmp(keyname,comkey,LEN_KEYWORD-1) != 0 && strlen(keyname) != 0 
         && strlen(strval) != 0)
        {
          fits_get_keytype(strval,&dtype,&status);
          if(status)  printerror(status);
          
          strip(keyname, 'B',' '); 
          strip(strval, 'B',' ');
          strip(strval, 'B','\'');
          
          switch( dtype )  
            { 
            case 'C':
              fits_read_key(fileptr,TSTRING,keyname,strval,comment,&status);
              dvl[keyname]= strval;
              dvl.SetComment(keyname, comment);
              break;
            case 'I':
              int ival;
              fits_read_key(fileptr,TINT,keyname,&ival,comment,&status);
                dvl[keyname]= (int_4) ival;             // Portage mac DY
                dvl.SetComment(keyname, comment);
              break;
            case 'L':
              int ilog;
              fits_read_key(fileptr,TLOGICAL,keyname,&ilog,comment,&status);
              dvl[keyname]= (int_4) ilog;
              dvl.SetComment(keyname, comment);
              break;
            case 'F':
              double dval;
              fits_read_key(fileptr,TDOUBLE,keyname,&dval,comment,&status);
              dvl[keyname]= dval;
              dvl.SetComment(keyname, comment);
              break;
            }
          
        }
    } 
  //  dvl_.Print();
} 

FitsOutFile::FitsOutFile()
{
  InitNull(); 
}

FitsOutFile::FitsOutFile(char flnm[], WriteMode wrm)
{

  InitNull();
  int status = 0; 

  // create new FITS file

  fits_create_file(&fptr_,flnm,&status);
  if( status )
    {

      switch (wrm)
        {
      // si on veut ecrire a la fin de ce fichier
        case append :
          status = 0;
          fits_open_file(&fptr_,flnm,READWRITE,&status);
          if( status )
            {
              cout << " error opening file: " << flnm << endl;
              printerror(status, "failure opening a file supposed to exist");
            }
          else cout << " file " << flnm << " opened, new objects will be appended " << endl;
          fits_get_num_hdus(fptr_, &hdunum_, &status);
          int hdutype;
          fits_movabs_hdu(fptr_,hdunum_,&hdutype,&status);
          if( status ) printerror( status,":FitsFile::WriteF : erreur movabs");
          break;
       
        case clear :
          {
            status = 0;
            char* newname = new char[strlen(flnm)+1];
            //
            newname[0] = '!';
            newname[1] = '\0';
            strcat(newname, flnm);
            fits_create_file(&fptr_,newname,&status);
            delete [] newname;
            if (status)
              {
                cout << " error opening file: " << flnm << endl;
                printerror(status, "unable to open file, supposed to exist");
              }
            else  cout << "  WARNING : file " << flnm << " is overwritten " << endl;
            break;
          }
        case unknown : 
          printerror(status, " file seems already to exist");
          break;
      
        }
    }
}



void FitsOutFile::makeHeaderImageOnFits(char type, int nbdim, int* naxisn,  DVList &dvl) 
{
  int status = 0;
  long naxis = nbdim;
  long* naxes = new long[nbdim];
  if (hdunum_ == 1)
    {
      if (imageOnPrimary_ == false)
        { 
          fits_create_img(fptr_,FLOAT_IMG,0,naxes,&status);
          writeSignatureOnFits();       }
      else  hdunum_--;
    }
  int k;
  for (k=0; k< nbdim; k++) naxes[k] = (long)naxisn[k];
  if (type == 'D')
    fits_create_img(fptr_,DOUBLE_IMG,naxis,naxes,&status);
  else
    if  (type == 'E')
      fits_create_img(fptr_,FLOAT_IMG,naxis,naxes,&status);
    else
      if  (type == 'I')
        fits_create_img(fptr_,LONG_IMG,naxis,naxes,&status);
      else
        {
          cout << " type of data: " << type << endl;
          throw PException("FitsFile:::makeHeaderImageOnFits:unprogrammed type of data ");
        }

  hdunum_++;

  // write supplementary keywords
  //  dvl.Print();
  addKeywordsOfDVList(dvl);

  delete [] naxes;
  if( status ) printerror( status, "erreur creation HDU IMAGE"  );

}
void FitsOutFile::putImageToFits(int nbData, double* map) const
{
  int status = 0;
  long npix= nbData;
  fits_write_img(fptr_,TDOUBLE,1,npix,map,&status); 
  if( status ) printerror( status, "erreur ecriture putImageToFits" );
  //  writeSignatureOnFits();
}

void FitsOutFile::putImageToFits(int nbData, float* map) const
{
  int status = 0;
  long npix= nbData;
  fits_write_img(fptr_,TFLOAT,1,npix, map,&status); 
  if( status ) printerror( status, "erreur ecriture putImageToFits" );
  //  writeSignatureOnFits();

}
void FitsOutFile::putImageToFits( int nbData, int* map) const
{
  int status = 0;

  long npix= nbData;
  fits_write_img(fptr_,TINT,1,npix,map,&status); 
  if( status ) printerror( status, "erreur ecriture putImageToFits" );
  //  writeSignatureOnFits();
}



void FitsOutFile::makeHeaderBntblOnFits( string fieldType, vector<string> Noms, int nentries, int tfields, DVList &dvl, string extname, vector<int> taille_des_chaines) 
{ 
  int status = 0;
  long nrows;
  if (fieldType.length() != tfields)
    {
      cout << " nombre de champs :" << tfields << "nombre de types: " << fieldType.length() << endl;
      throw ParmError("FitsFile:: fields and types don't match");

    }
  char ** ttype= new char*[tfields];
  char ** tform= new char*[tfields];
  char largeur[FLEN_VALUE];
  int noColString=0;
  int k;
  for (k=0; k<tfields;k++)
    {
      char  format[FLEN_VALUE];
      
      if(nentries < 1024) 
        {
          nrows= nentries;
          if (fieldType[k] == 'A')
            {
              sprintf(largeur,"%d",taille_des_chaines[noColString++]);
              strcpy(format,largeur);
            }
          else strcpy(format,"1");
        }
      else
        {
          nrows = nentries/1024;
          if(nentries%1024 != 0) nrows++;
          if (fieldType[k] == 'A') 
            {
              char largaux[FLEN_VALUE];
              sprintf(largeur,"%d",taille_des_chaines[noColString]);
              sprintf(largaux,"%d",1024*taille_des_chaines[noColString]);
              noColString++;
              strcpy(format, largaux);    
            }
          else strcpy(format,"1024");
        }
      strncat(format,&fieldType[k],1);
      if (fieldType[k] == 'A') 
        {
          strcat(format,largeur);
        }
      ttype[k] =  const_cast<char*>(Noms[k].c_str());
      tform[k]= new char[FLEN_VALUE];
      strcpy(tform[k],format);
    }
  char* extn = const_cast<char*>(extname.c_str());

  // create a new empty binary table onto the FITS file
  // physical units if they exist, are defined in the DVList object
  // so the NULL pointer is given for the tunit parameters.
  nrows=0;
  fits_create_tbl(fptr_,BINARY_TBL,nrows,tfields,ttype,tform,
                  NULL,extn,&status);
  if( status ) printerror( status );
  if ( hdunum_ == 0 ) hdunum_ = 2;
  else  hdunum_++;
  int ii;
  for(ii = 0; ii < tfields; ii++)
    {
      delete [] tform[ii];
    }
  delete [] ttype;
  delete [] tform;
  //
  // write supplementary keywords
  addKeywordsOfDVList(dvl);
}

void FitsOutFile::putColToFits(int nocol, int nentries, double* donnees) const
{
  int status = 0;
  int hdutype;
  fits_movabs_hdu(fptr_,hdunum_,&hdutype,&status);
  if( status ) printerror(status,"putColToFits: le movabs a foire");
  fits_get_hdu_type(fptr_, &hdutype, &status);
    if(hdutype != ASCII_TBL && hdutype != BINARY_TBL)
        {
            cout << " hdunum= " << hdunum_ << " hdutype= " << hdutype << endl;
     throw IOExc("FitsFile::putColToFits, this HDU is not an ASCII table nor a binary table");
          }
  int code; 
  long repeat, width;
  fits_get_coltype(fptr_, nocol+1, &code, &repeat,&width, &status);
  if( code != TDOUBLE)
    {
      cout << " WARNING : types don't match (putColToFits) : on fits file= " << code << " to be written= DOUBLE " << endl;
    } 
  fits_write_col(fptr_,TDOUBLE,nocol+1,1,1,nentries, donnees ,&status);
  if( status )  printerror( status,"erreur ecriture du fichier fits" );
}
void FitsOutFile::putColToFits(int nocol, int nentries, float* donnees) const
{
  int status = 0;
  int hdutype;
  fits_movabs_hdu(fptr_,hdunum_,&hdutype,&status);
  if( status ) printerror(status,"putColToFits: le movabs a foire");
  fits_get_hdu_type(fptr_, &hdutype, &status);
  if(hdutype != ASCII_TBL && hdutype != BINARY_TBL)
    {
      cout << " hdunum= " << hdunum_ << " hdutype= " << hdutype << endl;
      throw IOExc("FitsFile::putColToFits, this HDU is not an ASCII table nor a binary table");
    }
  if(hdutype == ASCII_TBL && nocol>0)
    {
      throw IOExc("FitsFile::putColToFits, this HDU is an ASCII table, nocol>0  forbidden");
    }
  int code; 
  long repeat, width;
  fits_get_coltype(fptr_, nocol+1, &code, &repeat,&width, &status);
  if( code != TFLOAT)
    {
      cout << " WARNING : types don't match (putColToFits) : on fits file= " << code << " (FITS code), to be written= FLOAT " << endl;
    } 
  fits_write_col(fptr_,TFLOAT,nocol+1,1,1,nentries, donnees ,&status);
  if( status )  printerror( status,"erreur ecriture du fichier fits" );
}
void FitsOutFile::putColToFits(int nocol, int nentries, int* donnees) const
{
  int status = 0;
  int hdutype;
  fits_movabs_hdu(fptr_,hdunum_,&hdutype,&status);
  if( status ) printerror(status,"putColToFits: le movabs a foire");
  fits_get_hdu_type(fptr_, &hdutype, &status);
  if(hdutype != ASCII_TBL && hdutype != BINARY_TBL)
    {
      cout << " hdunum= " << hdunum_ << " hdutype= " << hdutype << endl;
      throw IOExc("FitsFile::putColToFits, this HDU is not an ASCII table nor a binary table");
    }
  if(hdutype == ASCII_TBL && nocol>0)
    {
      throw IOExc("FitsFile::putColToFits, this HDU is an ASCII table, nocol>0  forbidden");
    }
  int code; 
  long repeat, width;
  fits_get_coltype(fptr_, nocol+1, &code, &repeat,&width, &status);
  if( code != TLONG &&  code != TINT &&  code != TSHORT ) 
    {
      cout << " WARNING : types don't match (putColToFits) : on fits file= " << code << " (FITS code), to be written= INT " << endl;
    } 
  fits_write_col(fptr_,TINT,nocol+1,1,1,nentries, donnees ,&status);
  if( status )  printerror( status," ecriture du fichier fits" );
}
void FitsOutFile::putColToFits(int nocol, int nentries, char** donnees) const
{
  int status = 0;
  int hdutype;
  fits_movabs_hdu(fptr_,hdunum_,&hdutype,&status);
  if( status ) printerror(status,"putColToFits: le movabs a foire");
  fits_get_hdu_type(fptr_, &hdutype, &status);
  if(hdutype != ASCII_TBL && hdutype != BINARY_TBL)
    {
      cout << " hdunum= " << hdunum_ << " hdutype= " << hdutype << endl;
      throw IOExc("FitsFile::putColToFits, this HDU is not an ASCII table nor a binary table");
    }
  if(hdutype == ASCII_TBL && nocol>0)
    {
      throw IOExc("FitsFile::putColToFits, this HDU is an ASCII table, nocol>0  forbidden");
    }
  int code; 
  long repeat, width;
  fits_get_coltype(fptr_, nocol+1, &code, &repeat,&width, &status);
  if( code != TSTRING)
    {
      cout << " WARNING : types don't match (putColToFits) : on fits file= " << code << " (FITS code), to be written= char** " << endl;
    } 
  fits_write_col(fptr_,TSTRING,nocol+1,1,1,nentries, donnees ,&status);
  if( status )  printerror( status,"erreur ecriture du fichier fits" );
}

void FitsOutFile::putBinTabLine(int NoLine, double* ddata, float* fdata, int* idata, char ** cdata) const
{
  int status= 0;
  int    anull;
  int dcount = 0.;
  int fcount = 0.;
  int icount = 0;
  int ccount =0;
  int ncol;
  long nels=1;
  int nbcols;
  fits_get_num_cols(fptr_, &nbcols,&status);
  for (ncol=0; ncol<nbcols; ncol++)
    {
      int code; 
      long repeat, width;
      fits_get_coltype(fptr_, ncol+1, &code, &repeat,&width, &status);
      switch (code)
        { 
        case TDOUBLE :
          fits_write_col(fptr_,TDOUBLE,ncol+1,NoLine+1,1,1, &ddata[dcount++] ,&status);
            break;
        case TFLOAT :
          fits_write_col(fptr_,TFLOAT,ncol+1,NoLine+1,1,1, &fdata[fcount++] ,&status);
          break;
        case TINT :
          fits_write_col(fptr_,TINT,ncol+1,NoLine+1,1,1, &idata[icount++] ,&status);
          break;
        case TLONG :
          fits_write_col(fptr_,TINT,ncol+1,NoLine+1,1,1, &idata[icount++] ,&status);
          break;
        case TSHORT :
          fits_write_col(fptr_,TINT,ncol+1,NoLine+1,1,1, &idata[icount++] ,&status);
          break;
        case TSTRING :
          fits_write_col(fptr_,TSTRING,ncol+1,NoLine+1,1,1, &cdata[ccount++] ,&status);
          break;
        }
      if (status)
        {
          cout << " WARNING : unknown type (putBinTabLine) : on fits file= " << code << " (FITS code) " << endl;
          status = 0;
          break;
        }
    }
}


void  FitsOutFile::DVListIntoPrimaryHeader(DVList& dvl) const
{
  int status = 0;
  int hdutype;
  fits_movabs_hdu(fptr_,1,&hdutype,&status);
  addKeywordsOfDVList(dvl);
  fits_movabs_hdu(fptr_,hdunum_,&hdutype,&status);
}


void FitsOutFile::writeSignatureOnFits() const
{
  int status = 0;
  char keyname[LEN_KEYWORD];
  char strval[FLEN_VALUE];
  char comment[FLEN_COMMENT];
  strncpy(keyname, "CREATOR", LEN_KEYWORD); 
  keyname[LEN_KEYWORD-1] = '\0';
  strcpy(strval, "SOPHYA");
  strcpy(comment," SOPHYA Package - FITSIOServer ");
  fits_write_key(fptr_, TSTRING, keyname, &strval, comment, &status);
  if( status ) printerror( status );
  fits_write_date(fptr_, &status);
  fits_write_comment(fptr_,"..............................................", &status);
  fits_write_comment(fptr_, " SOPHYA package - FITSIOSever ", &status);
  fits_write_comment(fptr_, " (C) LAL/IN2P3-CNRS Orsay,  FRANCE 2000", &status);
  fits_write_comment(fptr_, " (C) DAPNIA/CEA     Saclay, FRANCE 2000", &status);
  fits_write_comment(fptr_,"..............................................", &status);
  if( status ) printerror( status, "erreur writeSignatureOnFits"  );
}


void FitsOutFile::addKeywordsOfDVList(DVList& dvl) const
{
  int status = 0;
  fits_write_comment(fptr_,"---------- keywords from SOPHYA ---------", &status);
  if (hdunum_ == 1) writeSignatureOnFits();
  DVList::ValList::const_iterator it;
  for(it = dvl.Begin(); it != dvl.End(); it++)
    {
      char keytype= (*it).second.elval.typ;
      char keyname[10];
      strncpy(keyname,(*it).first.substr(0,64).c_str(),10);
      string key(keyname);
      char comment[FLEN_COMMENT];
      char strval[FLEN_VALUE]= "";
      char *comkey = "COMMENT";
      fits_read_keyword(fptr_, keyname, strval, NULL, &status);
      if (status != 0 || strncmp(keyname,comkey,LEN_KEYWORD-1) == 0 ) 
        {
          string coco = dvl.GetComment(key);
          coco.copy( comment, FLEN_COMMENT-1);
          int bout = (coco.length() < FLEN_COMMENT) ? coco.length() : FLEN_COMMENT-1;
          comment[bout]= '\0';
          status = 0;
          switch (keytype)
            {
            case 'I' :
              {
                int ival = (int)dvl.GetI(key);
                fits_write_key(fptr_,TINT,keyname,&ival, comment,&status);
                break;
              }
            case 'D' :
              {
                double dval= (double)dvl.GetD(key);
                fits_write_key(fptr_,TDOUBLE,keyname,&dval,comment,&status);
                break;
              }
            case 'S' :
              {
                char strvaleur[FLEN_VALUE]= "";
                string valChaine = dvl.GetS(key);
                valChaine.copy(strvaleur, FLEN_VALUE-1);
                int fin =  (valChaine.length() < FLEN_VALUE) ? valChaine.length() : FLEN_VALUE-1;
                strvaleur[fin]= '\0';

                fits_write_key(fptr_,TSTRING,keyname,&strvaleur,comment,&status);
                break;
              }
            } 
        }
  if( status )  printerror( status,"fitsfile: probleme ecriture mot-cle du dvlist" );
    }
  fits_write_comment(fptr_,"--------------------------------------", &status);
}