source: Sophya/trunk/SophyaExt/FitsIOServer/fitsioserver.cc@ 687

//************************************************************************
// Class for loadind and saving from FITS-formatted file to DPC objects
//    (G. Le Meur ; Francois Touze)  OCT. 99
//
// methods 'load(X& x, char f[])' get from FITS file "f" a DPC object x 
// from DPC  class X.
// methods 'save(X& x, char f[])' save  a DPC object x from DPC // class X 
// onto a FITS file "f" .

//************************************************************************

#include "machdefs.h" 
#include <iostream.h>
#include <list>
#include <string>

#include "fitsioserver.h" 
#include "pexceptions.h" 
#include "strutil.h"

void FitsIoServer::load(TMatrix<double>& mat,char flnm[]) 
{
  int nbrows=0;
  int nbcols=0;
  FITS_tab_typ_ = TDOUBLE;
  long naxis;
  int n1, n2, n3;
  DVList dvl;
  planck_read_img(flnm, naxis, n1, n2, n3, dvl);

  nbrows=n1;
  nbcols=n2;
  if (naxis == 1) nbcols=1;
  if (naxis > 2 && n3>1) 
    {
      cout<<"  FitsIOServer : le fichier fits n'est pas une matrice, naxis= " <<naxis<< endl;
    }

  // number of components
  if (mat.NRows() != nbrows || mat.NCols() != nbcols )
    {
      //    cout << " found " << nbrows << " rows ";
      //   cout << "  expected " << mat.NRows()  << endl;
      //  cout << " found " << nbcols << " columns " ;
      // cout << "  expected " <<  mat.NCols()  << endl;
      mat.ReSize(nbrows,nbcols);
      cout << " resize the vector to nbrows=  " << nbrows  << " nbcols= " << nbcols << endl;
    }
  int ij=0;
  for (int j=0; j< nbcols; j++) 
    for (int i = 0; i < nbrows; i++)  mat(i,j) = (double)r_8tab_[ij++];
}

void FitsIoServer::load(TMatrix<float>& mat,char flnm[]) 
{
  int nbrows=0;
  int nbcols=0;
  FITS_tab_typ_ = TFLOAT;
  long naxis;
  int n1, n2, n3;
  DVList dvl;
  planck_read_img(flnm, naxis, n1, n2, n3, dvl);

  nbrows=n1;
  nbcols=n2;
  if (naxis == 1) nbcols=1;
  if (naxis > 2 && n3>1) 
    {
      cout<<"  FitsIOServer : le fichier fits n'est pas une matrice, naxis= " <<naxis<< endl;
    }

  // number of components
  if (mat.NRows() != nbrows || mat.NCols() != nbcols )
    {
      //      cout << " found " << nbrows << " rows ";
      //    cout << "  expected " << mat.NRows()  << endl;
      //  cout << " found " << nbcols << " columns " ;
      //  cout << "  expected " <<  mat.NCols()  << endl;
      mat.ReSize(nbrows,nbcols);
      cout << " resize the vector to nbrows=  " << nbrows  << " nbcols= " << nbcols << endl;
    }
  int ij=0;
  for (int j=0; j< nbcols; j++) 
    for (int i = 0; i < nbrows; i++)  mat(i,j) = (float)r_4tab_[ij++];
}

void FitsIoServer::load(TMatrix<int_4>& mat,char flnm[]) 
{
  int nbrows=0;
  int nbcols=0;
  FITS_tab_typ_ = TINT;
  long naxis;
  int n1, n2, n3;
  DVList dvl;
  planck_read_img(flnm, naxis, n1, n2, n3, dvl);

  nbrows=n1;
  nbcols=n2;
  if (naxis == 1) nbcols=1;
  if (naxis > 2 && n3>1) 
    {
      cout<<"  FitsIOServer : le fichier fits n'est pas une matrice, naxis= " <<naxis<< endl;
    }

  // number of components
  if (mat.NRows() != nbrows || mat.NCols() != nbcols )
    {
      //     cout << " found " << nbrows << " rows ";
      //   cout << "  expected " << mat.NRows()  << endl;
      //  cout << " found " << nbcols << " columns " ;
      //  cout << "  expected " <<  mat.NCols()  << endl;
      mat.ReSize(nbrows,nbcols);
      cout << " resize the vector to nbrows=  " << nbrows  << " nbcols= " << nbcols << endl;
    }
  int ij=0;
  for (int j=0; j< nbcols; j++) 
    for (int i = 0; i < nbrows; i++)  mat(i,j) = (int_4)i_4tab_[ij++];
}

void FitsIoServer::load(NTuple& ntpl,char flnm[],int hdunum) 

  //********************************************************/
  //* move to the HDU which has the specified number hdunum*/
  //* in the FITS file, perform read operations and write  */
  //* the elements in an NTuple.                           */
  //********************************************************/
{

  // pointer to the FITS file, defined in fitsio.h 
  fitsfile *fptr;
  int status = 0;
  if( fits_open_file(&fptr,flnm,READONLY,&status) ) 
    printerror( status );

  // move to the HDU
  int hdutype;
  if( fits_movabs_hdu(fptr,hdunum,&hdutype,&status) ) 
    printerror( status );

  // get number of keywords
  int nkeys,keypos;
  if( fits_get_hdrpos(fptr,&nkeys,&keypos,&status) )
    printerror( status );


  if (hdutype == BINARY_TBL)
    printf("\nReading binary table in HDU %d:\n",hdunum);
  else
    {
      printf("Error:: this HDU is not a binary table\n");
      throw IOExc("FitsIoServer::load(NTuple& ," + (string)flnm + ") Error Not a bin table !");
//      exit ( status );
    }

  // number of columns
  int tfields;
  if( fits_get_num_cols(fptr,&tfields,&status) )
    printerror( status );

  // to get table size
  long naxis[2];
  int nfound;
   if( fits_read_keys_lng(fptr, "NAXIS", 1, 2, naxis, &nfound, &status) ) 
     printerror( status );
   int nrows= naxis[1];

   //Information about each column
   char **ttype, **tform;
   ttype= new char*[tfields];
   tform= new char*[tfields];
   int ii;
   for(ii = 0; ii < tfields; ii++)
     {
       ttype[ii]= new char[FLEN_VALUE];
       tform[ii]= new char[FLEN_VALUE];
     }

   // number of TTYPEn,TFORMn and TUNITn keywords found
   int num= 0;

   // read the column names from the TTYPEn keywords
   fits_read_keys_str(fptr,"TTYPE",1,tfields,ttype,&nfound,&status);
   num += nfound;
   // read the column names from the TFORMn keywords
   fits_read_keys_str(fptr,"TFORM",1,tfields,tform,&nfound,&status);
   num += nfound;

   //for(int ii = 0; ii < tfields; ii++)
   //  printf("\nColumn name & Format %8s %8s", ttype[ii],tform[ii]);
  
   // select the columns with float data values 
   int typecode;
   long repeat,width;
   list<int> column;
   for(ii = 0; ii < tfields; ii++)
     {
       fits_binary_tform(tform[ii],&typecode, &repeat, &width, &status);
       //printf("\n%4s %3d %2ld %2ld", tform[ii], typecode, repeat, width);
       if(typecode == TFLOAT) 
         column.push_back(ii+1);
     }

   // get input elements to create the NTuple
   char **clname;
   clname= new char*[column.size()];
   for(ii = 0; ii < column.size(); ii++)
     clname[ii]= new char[FLEN_VALUE];

   list<int>::iterator itr;
   int index= 0;
   for(itr= column.begin(); itr != column.end(); itr++)
     strcpy(clname[index++],ttype[*itr-1]);

   for(ii = 0; ii < tfields; ii++)
     {
       delete [] ttype[ii];
       delete [] tform[ii];
     }
   delete [] ttype;
   delete [] tform;
   

   // check if the specified keyword BLK exists
   int blk= 512;
   if(check_keyword(fptr,nkeys,"BLK"))
     {
       if( fits_read_key(fptr,TINT,"BLK",&blk,NULL,&status) )
         printerror( status );
     }
   // create a NTuple
   NTuple nt0(column.size(),clname,blk);

   for(ii = 0; ii < column.size(); ii++)
     delete [] clname[ii];
   delete [] clname;

   float value[1];
   long felem = 1;
   char strnull[10];
   strcpy(strnull, " ");
   int anynull= 0;

   // value to represent undefined array elements
   float floatnull = FLOATNULLVALUE;

   // read elements from columns and fill NTuple
   for (int k = 0; k < nrows; k++)
     {
       int j= 0;
       float *xnt= new float[column.size()];
       list<int>::iterator itr;
       for(itr= column.begin(); itr != column.end(); itr++)
         {
           fits_read_col(fptr,TFLOAT,*itr,k+1,felem,1,&floatnull,value,
                         &anynull, &status);
           xnt[j++]= value[0];
         }
       nt0.Fill(xnt);
       delete[] xnt;
     }

   // the TUNITn keywords are optional, if they exist they are put 
   // in the DVLIst object
   char keyname[LEN_KEYWORD]= "";
   char strval[FLEN_VALUE];
   for(ii = 0; ii < tfields; ii++)
     {
       fits_make_keyn("TUNIT",ii+1,keyname,&status);
       if(check_keyword(fptr,nkeys,keyname))
         {
           num++;
           if( fits_read_key_str(fptr,keyname,strval,NULL,&status) )
             printerror(status);
           strip (keyname, 'B',' '); 
           strip(strval, 'B',' ');
           nt0.Info()[keyname]= strval;
         }
     }

   // add the number of mandatory keywords of a binary table
   num += 8;

   // put names and values of other reserved keywords in a DVList object
   char comment[FLEN_COMMENT];
   char dtype;
   for(int j = num+1; j <= nkeys; j++)
     {
       fits_read_keyn(fptr,j,keyname,strval,comment,&status);
       fits_get_keytype(strval,&dtype,&status);
       strip (keyname, 'B',' '); 
       strip(strval, 'B',' ');
       strip(strval, 'B','\'');
       
       switch( dtype ) 
         {
         case 'C':
           nt0.Info()[keyname]= strval;   
           break;
         case 'I':
           int ival;
           ctoi(strval,&ival);
           nt0.Info()[keyname]= (int_4)ival;
           break;
         case 'F':
           double dval;
           ctof(strval,&dval);
           nt0.Info()[keyname]= dval;
           break;
         }
     }

   // copy in the input NTuple
   ntpl= nt0;

   if( fits_close_file(fptr, &status) )
     printerror(status);

   printf("\n");
   return;
}


void FitsIoServer::load(SphericalMap<double>& sph, char flnm[]) 
{
  int npixels=0;
  int nside=0;
  long naxis;
  int n1, n2, n3;

  FITS_tab_typ_ = TDOUBLE;

  DVList dvl;
  planck_read_img(flnm, naxis, n1, n2, n3, dvl);
  if (naxis != 1) 
    {
      cout << " le fichier fits n'est pas une sphere, naxis= " << naxis << endl;
    }
  npixels=n1;
  nside= dvl.GetI("NSIDE");
  
  // number of pixels in the sphere
  if (sph.NbPixels() != npixels)
    {
      //DBG      cout << " found " << npixels << " pixels" << endl;
      //DBG cout << " expected " << sph.NbPixels() << endl;
      if (nside <= 0 ) 
        {
          cout<<" FITSIOSERVER: no resolution parameter on fits file "<<endl;
          throw IOExc("FitsIoServer::load(SphericalMap<double>& ," + (string)flnm + ") Error NSide<0 !");
//        exit(0);
        }
      if (nside != sph.SizeIndex())
        { 
          sph.Resize(nside);
          cout << "FitsIoServer::load(SphereGorski<double> ...) ReSizing to NSide= " << nside  << endl;
        }
      //      else
      //        {
      //          cout << " FITSIOSERVER : same resolution, surprising ! " << endl;
//        exit(0);   $CHECK$   Ca peut etre OK , non ??  Reza 20/11/99
      //        }
    }
  for (int j = 0; j < sph.NbPixels(); j++) sph(j)= (double)r_8tab_[j];
}


void FitsIoServer::load(SphereGorski<float>& sph, char flnm[], int hdunum) 
{
  int npixels=0;
  int nside=0;

  FITS_tab_typ_ = TFLOAT;

  DVList dvl;
  planck_read_bntbl(flnm, hdunum, npixels, dvl);
  //dvl.Print();
  nside= dvl.GetI("NSIDE");
  const char* ordering= dvl.GetS("ORDERING").c_str();
  char* ring = "RING";
  if (  strncmp(ordering, ring,4) != 0) 

    //  if (ordering!="RING    ")
    {
      cout << " numerotation non RING" << endl;
    }
  // number of pixels in the sphere
  if (sph.NbPixels() != npixels)
    {
      //DBG      cout << " found " << npixels << " pixels" << endl;
      //DBG      cout << " expected " << sph.NbPixels() << endl;
      if (nside <= 0 )  
        {
          cout<<" FITSIOSERVER: no resolution parameter on fits file "<<endl;
          throw IOExc("FitsIoServer::load(SphericalMap<float>& ," + (string)flnm + ", ) Error No resolution parameter !");
//        exit(0);
        }
      if (nside != sph.SizeIndex())
        { 
          sph.Resize(nside);
          cout << "FitsIoServer::load(SphereGorski<float> ...) ReSizing to NSide= " << nside  << endl;
        }
      //      else
      //        {
      //          cout << " FITSIOSERVER : same resolution, surprising ! " << endl;
//        exit(0);   $CHECK$ - Il ne faut pas sortir, me semble-t-il , Reza 20/11/99
      //        }
    }
  for (int j = 0; j < sph.NbPixels(); j++) 
    {
      sph(j)= (float)r_4tab_[j];
    }
}
void FitsIoServer::load(SphereGorski<double>& sph, char flnm[], int hdunum) 
{
  int npixels=0;
  int nside=0;

  FITS_tab_typ_ = TDOUBLE;

  DVList dvl;
  planck_read_bntbl(flnm, hdunum, npixels, dvl);
  //dvl.Print();
  nside= dvl.GetI("NSIDE");
  const char* ordering= dvl.GetS("ORDERING").c_str();
  char* ring = "RING";
  if (  strncmp(ordering, ring,4) != 0) 

    //  if (ordering!="RING    ")
    {
      cout << " numerotation non RING" << endl;
    }
  // number of pixels in the sphere
  if (sph.NbPixels() != npixels)
    {
      cout << " found " << npixels << " pixels" << endl;
      cout << " expected " << sph.NbPixels() << endl;
      if (nside <= 0 ) 
        {
          cout<<" FITSIOSERVER: no resolution parameter on fits file "<<endl;
          throw IOExc("FitsIoServer::load(SphereGorski<double>& ," + (string)flnm + ", )  No resol parameter !");
//        exit(0);
        }
      if (nside != sph.SizeIndex())
        { 
          sph.Resize(nside);
          cout << " resizing the sphere to nside=  " << nside  << endl;
        }
      else
        {
          cout << " FITSIOSERVER : same resolution, surprising ! " << endl;
//        exit(0);   $CHECK$ - ne pas sortir , Reza 20/11/99
        }
    }
  for (int j = 0; j < sph.NbPixels(); j++) sph(j)= (double)r_8tab_[j];
}

void FitsIoServer::load(SphericalMap<float>& sph, char flnm[]) 
{
  int npixels=0;
  int nside=0;
  long naxis;
  int n1, n2, n3;
  DVList dvl;

  FITS_tab_typ_ = TFLOAT;

  planck_read_img(flnm, naxis, n1, n2, n3, dvl);
  if (naxis != 1) 
    {
      cout << " le fichier fits n'est pas une sphere, naxis= " << naxis << endl;
    }
  npixels=n1;
  nside= dvl.GetI("NSIDE");

  // number of pixels in the sphere
  if (sph.NbPixels() != npixels)
    {
      cout << " found " << npixels << " pixels" << endl;
      cout << " expected " << sph.NbPixels() << endl;
      if (nside <= 0 ) 
        {
          cout<<" FITSIOSERVER: no resolution parameter on fits file "<<endl;
          throw IOExc("FitsIoServer::load(SphericalMap<float>& ," + (string)flnm + ") No resolution param !");
//        exit(0);
        }
      if (nside != sph.SizeIndex())
        { 
          sph.Resize(nside);
          cout << " resizing the sphere to nside=  " << nside  << endl;
        }
      else
        {
          cout << " FITSIOSERVER : same resolution, surprising ! " << endl;
//        exit(0);  $CHECK$ - Ne pas sortir , Reza 20/11/99
        }
    }
  for (int j = 0; j < sph.NbPixels(); j++) sph(j)= (float)r_4tab_[j];
}

void FitsIoServer::load(LocalMap<double>& lcm, char flnm[]) 
{
  int nbrows=0;
  int nbcols=0;
  FITS_tab_typ_ = TDOUBLE;
  long naxis;
  int n1, n2, n3;
  DVList dvl;
  planck_read_img(flnm, naxis, n1, n2, n3, dvl);

  nbrows=n1;
  nbcols=n2;
  if (naxis != 2) 
    {
      cout<<" FitsIOServer : le fichier fits n'est pas une localmap, naxis= "<<naxis<<endl;
    }
  float theta0 = dvl.GetD("THETA0");
  float phi0 = dvl.GetD("PHI0");
  int x0 = dvl.GetI("X0");
  int y0 = dvl.GetI("Y0");
  int xo = dvl.GetI("XO");
  int yo = dvl.GetI("YO");
  float anglex=dvl.GetD("ANGLEX");
  float angley=dvl.GetD("ANGLEY");
  float angle = dvl.GetD("ANGLE");

  // number of components

  if (lcm.Size_x() != nbrows || lcm.Size_y() != nbcols )
    {
      cout << " found " << nbrows << " x-pixels ";
      cout << "  expected " << lcm.Size_x()  << endl;
      cout << " found " << nbcols << " y-pixels " ;
      cout << "  expected " <<   lcm.Size_y()  << endl;
      lcm.ReSize(nbrows,nbcols);
      cout << " resizing the map to x-size=  " << nbrows  << " y-size= " << nbcols << endl;
    }
  lcm.SetSize(anglex, angley);
  lcm.SetOrigin(theta0, phi0, x0, y0, angle);
  int ij=0;
  for (int j=0; j< nbcols; j++) 
    for (int i = 0; i < nbrows; i++)  lcm(i,j) = (double)r_8tab_[ij++];
}

void FitsIoServer::load(ImageR4& DpcImg,char flnm[])
{
  FITS_tab_typ_ = TFLOAT;
  long naxis;
  int siz_x;
  int siz_y;
  int n3;
  DVList dvl;

  planck_read_img(flnm, naxis, siz_x, siz_y, n3, dvl);

  
  if (naxis != 2) 
    {
      cout << " FitsIOServer : naxis= " << naxis << " not equal to 2 ?" << endl;
    }

  DpcImg.Allocate(siz_x, siz_y, 0., 0);
  float* pixelPtr=DpcImg.ImagePtr();
  // verifications de type
  PBaseDataTypes  dataT=DataType((r_4)0);
  int TypeDonnees=dvl.GetI("DATATYPE");
  if (int(dataT) !=  TypeDonnees)
    {
      cout << " FitsIOServer : parameter DATATYPE on file " << flnm << " is not float " << endl;
      cout << "  eventual conversion to float is achieved by cfitsio lib " << endl;
    }

  memcpy(pixelPtr, r_4tab_, siz_x*siz_y*DataSize(dataT));
  const char* nom=dvl.GetS("NAME").c_str();
  DpcImg.SetNameId(dvl.GetI("IDENT"), nom);
  DpcImg.SetOrg(dvl.GetI("ORG_X"), dvl.GetI("ORG_Y"));
  DpcImg.SetPxSize(dvl.GetD("PIXSZ_X"), dvl.GetD("PIXSZ_Y"));    
  DpcImg.SetAtt(dvl.GetI("NBNUL"), dvl.GetI("NBSAT"),
       dvl.GetD("MINPIX"), dvl.GetD("MAXPIX"), dvl.GetD("MOYPIX"), 
       dvl.GetD("SIGPIX"),dvl.GetD("FOND"), dvl.GetD("SIGFON"));

}


void FitsIoServer::load(ImageI4& DpcImg,char flnm[])
{
  FITS_tab_typ_ = TINT;
  long naxis;
  int siz_x;
  int siz_y;
  int n3;
  DVList dvl;

  // pointer to the FITS file, defined in fitsio.h
  //fitsfile *fptr;
  // initialize status before calling fitsio routines
  //int status = 0; 
        
  //fits_open_file(&fptr, flnm, READONLY, &status);
  //if( status )  printerror( status );
  planck_read_img(flnm, naxis, siz_x, siz_y, n3, dvl);
  // close the file
  //fits_close_file(fptr, &status);
  //if( status )  printerror( status );
  
  if (naxis != 2) 
    {
      cout << " FitsIOServer : naxis= " << naxis << " not equal to 2 ?" << endl;
    }

  DpcImg.Allocate(siz_x, siz_y, 0., 0);
  int_4* pixelPtr=DpcImg.ImagePtr();


  // verifications de type
  PBaseDataTypes  dataT=DataType((int_4)0);
  int TypeDonnees=dvl.GetI("DATATYPE");
  if (int(dataT) !=  TypeDonnees)
    {
      cout << " FitsIOServer : parameter DATATYPE on file " << flnm << " is not int_4 " << endl;
      cout << "  eventual conversion to int_4 is achieved by cfitsio lib " << endl;
    }
  memcpy(pixelPtr, i_4tab_, siz_x*siz_y*DataSize(dataT));
  const char* nom=dvl.GetS("NAME").c_str();
  DpcImg.SetNameId(dvl.GetI("IDENT"), nom);
  DpcImg.SetOrg(dvl.GetI("ORG_X"), dvl.GetI("ORG_Y"));
  DpcImg.SetPxSize(dvl.GetD("PIXSZ_X"), dvl.GetD("PIXSZ_Y"));    
  DpcImg.SetAtt(dvl.GetI("NBNUL"), dvl.GetI("NBSAT"),
       dvl.GetD("MINPIX"), dvl.GetD("MAXPIX"), dvl.GetD("MOYPIX"), 
       dvl.GetD("SIGPIX"),dvl.GetD("FOND"), dvl.GetD("SIGFON"));
}


void FitsIoServer::save( TMatrix<double>& mat, char filename[])   
{
  int nbrows = mat.NRows();
  int nbcols = mat.NCols();
  long naxis = nbcols > 1 ? 2 : 1; 
  cout << " nombre de lignes : " << nbrows << " colonnes " << nbcols << endl;
  FITS_tab_typ_ = TDOUBLE;
  if (r_8tab_ != NULL ) { delete[] r_8tab_;  r_8tab_ = NULL; }
  r_8tab_=new r_8[nbrows*nbcols];


  int ij=0;
  for (int j=0; j< nbcols; j++) 
    for (int i = 0; i < nbrows; i++) r_8tab_[ij++]= (r_8)mat(i,j);

  DVList dvl;

  planck_write_img(filename, naxis, nbrows, nbcols, 0, dvl);
  delete[] r_8tab_;  r_8tab_ = NULL; 
}

void FitsIoServer::save( TMatrix<float>& mat, char filename[])   
{
  int nbrows = mat.NRows();
  int nbcols = mat.NCols();
  long naxis = nbcols > 1 ? 2 : 1; 
  cout << " nombre de lignes : " << nbrows << " colonnes " << nbcols << endl;
  FITS_tab_typ_ = TFLOAT;
  if (r_4tab_ != NULL ) { delete[] r_4tab_;  r_4tab_ = NULL; }  
  r_4tab_=new r_4[nbrows*nbcols];


  int ij=0;
  for (int j=0; j< nbcols; j++) 
    for (int i = 0; i < nbrows; i++) r_4tab_[ij++]= (r_4)mat(i,j);

  DVList dvl;

  planck_write_img(filename, naxis, nbrows, nbcols, 0, dvl);
  delete[] r_4tab_;  r_4tab_ = NULL;
}

void FitsIoServer::save( TMatrix<int_4>& mat, char filename[])   
{
  int nbrows = mat.NRows();
  int nbcols = mat.NCols();
  long naxis = nbcols > 1 ? 2 : 1; 
  cout << " nombre de lignes : " << nbrows << " colonnes " << nbcols << endl;
  FITS_tab_typ_ = TINT;
  if (i_4tab_ != NULL ) { delete[]  i_4tab_;   i_4tab_ = NULL; }
  i_4tab_=new int_4[nbrows*nbcols];


  int ij=0;
  for (int j=0; j< nbcols; j++) 
    for (int i = 0; i < nbrows; i++) i_4tab_[ij++]= (int_4)mat(i,j);

  DVList dvl;

  planck_write_img(filename, naxis, nbrows, nbcols, 0, dvl);
  delete[] i_4tab_;  i_4tab_ = NULL;
}


void FitsIoServer::save(NTuple& ntpl,char flnm[])

  //****************************************************/
  //* read the elements of the NTuple ntpl, and create */ 
  //* a FITS file with a binary table extension        */
  //****************************************************/
{

  // delete old file if it already exists
  remove(flnm);

  // pointer to the FITS file, defined in fitsio.h 
  fitsfile *fptr;
  int status = 0;
  if( fits_create_file(&fptr, flnm, &status) )
    printerror( status );

  // table will have tfields columns
  int tfields= ntpl.NVar();

  // table will have nrows rows
  int nrows= ntpl.NEntry();

  // extension name
  char * extname = "NTuple_Binary_tbl";

  // define the name, and the datatype for the tfields columns
  char **ttype, **tform;
  ttype= new char*[tfields];
  tform= new char*[tfields];
  int i;
  for(i = 0; i < tfields; i++)
    {
      ttype[i]= new char[FLEN_VALUE];
      strcpy(ttype[i], ntpl.NomIndex(i));
      tform[i]= new char[FLEN_VALUE];
      strcpy(tform[i], "1E");
    }

  // 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.
  if ( fits_create_tbl(fptr,BINARY_TBL,nrows,tfields,ttype,tform,
                       NULL,extname,&status) )
    printerror( status );

  for( int ii = 0; ii < tfields; ii++)
    {
      delete [] ttype[ii];
      delete [] tform[ii];
    }
  delete [] ttype;
  delete [] tform;


  // first row in table to write
  int firstrow  = 1; 

  // first element in row  (ignored in ASCII tables)
  int firstelem = 1;

  for(i = 0; i < tfields; i++)
    {
      float *dens= new float[nrows];
      for(int j = 0; j < nrows; j++)
        {
          dens[j]= ntpl.GetVal(j,i);
        }
      fits_write_col(fptr,TFLOAT,i+1,firstrow,firstelem,nrows,dens,&status);
      delete[] dens;
    }

  // number of blocks per event
  int blk= ntpl.BLock();
  fits_write_key(fptr,TINT,"BLK",&blk,"number of blocks per evt",&status);

  // get names and values from the join DVList object
  DVList dvl= ntpl.Info();
  dvl.Print();
  DVList::ValList::const_iterator it;
  for(it = dvl.Begin(); it != dvl.End(); it++)
    {
      char keytype= (*it).second.typ;
      char keyname[10];
      strncpy(keyname,(*it).first.substr(0,64).c_str(),8);
      char comment[FLEN_COMMENT];

      switch (keytype)
        {
        case 'I' :
          {
            int ival=(*it).second.mtv.iv;
            strcpy(comment,"I entier");
            fits_write_key(fptr,TINT,keyname,&ival,comment,&status);
            break;
          }
        case 'D' :
          {
            double dval=(*it).second.mtv.dv;
            strcpy(comment,"D double");
            fits_write_key(fptr,TDOUBLE,keyname,&dval,comment,&status);
            break;
          }
        case 'S' :
          {
            char strval[128];
            strncpy(strval,(*it).second.mtv.strv,127);
            strcpy(comment,"S character string");
            fits_write_key(fptr,TSTRING,keyname,&strval,comment,&status);
            break;
          }
        } 
    }

  //close the FITS file
  if ( fits_close_file(fptr, &status) )
    printerror( status );
  return;
}  




//void FitsIoServer::save(SphericalMap<double>& sph, char filename[])   
void FitsIoServer::save(SphericalMap<double>& sph, char filename[])   
{
  int npixels = sph.NbPixels();
  FITS_tab_typ_ = TDOUBLE;
  if (r_8tab_ != NULL ) { delete[] r_8tab_;  r_8tab_ = NULL; }
  r_8tab_=new r_8[npixels];


  for (int j = 0; j < npixels; j++) r_8tab_[j]= (r_8)sph(j);
  DVList dvl;
  dvl["NSIDE"] = (int_4) sph.SizeIndex();
  dvl["ORDERING"]=sph.TypeOfMap();

  planck_write_img(filename, 1, npixels, 0, 0, dvl);

  // decider ulterieurement ce qu'on fait de ce qui suit, specifique
  // pour l'instant, aux spheres gorski.

  /*
  // write supplementary keywords
  fits_write_comment(fptr, " ", &status);
  if( status ) printerror( status );


  strcpy(comment,"HEALPIX Pixelisation");
  strcpy(svalue, "HEALPIX");
  fits_write_key(fptr, TSTRING, "PIXTYPE", svalue, comment, &status);
  if( status )  printerror( status );


  strcpy(comment,"pixel ordering scheme, either RING or NESTED");
  strcpy(svalue, "RING");
  fits_write_key(fptr, TSTRING, "ORDERING", svalue, comment, &status);
  if( status )  printerror( status );


  strcpy(comment,"Random generator seed");
  int iseed= sph.iseed();
  fits_write_key(fptr, TINT, "RANDSEED", &iseed, comment, &status);
  if( status )  printerror( status );


  strcpy(comment,"--------------------------------------------");
  fits_write_comment(fptr, comment, &status);
  if( status ) printerror( status );

  strcpy(comment," Above keywords are still likely to change !");
  fits_write_comment(fptr, comment, &status);
  if( status ) printerror( status );

  strcpy(comment,"--------------------------------------------");
  fits_write_comment(fptr, comment, &status);
  if( status ) printerror( status );
  
  */

 delete[] r_8tab_;  r_8tab_ = NULL;
}

void FitsIoServer::save(SphericalMap<float>& sph, char filename[])   
{
  int npixels = sph.NbPixels();
  FITS_tab_typ_ = TFLOAT;
  if (r_4tab_ != NULL ) { delete[] r_4tab_;  r_4tab_ = NULL; }
  r_4tab_=new r_4[npixels];


  for (int j = 0; j < npixels; j++) r_4tab_[j]= (r_4)sph(j);
  DVList dvl;
  dvl["NSIDE"] = (int_4)sph.SizeIndex();
  dvl["ORDERING"]=sph.TypeOfMap();

  planck_write_img(filename, 1, npixels, 0, 0, dvl);
  delete[] r_4tab_;  r_4tab_ = NULL;

}
void FitsIoServer::save(SphereGorski<float>& sph, char filename[])   
{
  int npixels = sph.NbPixels();
  FITS_tab_typ_ = TFLOAT;
  if (r_4tab_ != NULL ) { delete[] r_4tab_;  r_4tab_ = NULL; }
  r_4tab_=new r_4[npixels];


  for (int j = 0; j < npixels; j++) r_4tab_[j]= (r_4)sph(j);
  DVList dvl;
  dvl.SetS("PIXTYPE","HEALPIX ");
  dvl["ORDERING"]=sph.TypeOfMap();
  dvl["NSIDE"] = (int_4)sph.SizeIndex();
  dvl["FIRSTPIX"]=(int_4)0;
  dvl["LASTPIX"]=(int_4)(npixels-1);
  char* typeOfContent="TEMPERATURE";
  char* extname="SIMULATION";
  char* comment1="     Sky Map Pixelisation Specific Keywords";
  planck_write_bntbl(filename, npixels, typeOfContent, extname, comment1, dvl);
  delete[] r_4tab_;  r_4tab_ = NULL;

}
void FitsIoServer::save(SphereGorski<double>& sph, char filename[])   
{
  int npixels = sph.NbPixels();
  FITS_tab_typ_ = TDOUBLE;
  if (r_8tab_ != NULL ) { delete[] r_8tab_;  r_8tab_ = NULL; }
  r_8tab_=new r_8[npixels];


  for (int j = 0; j < npixels; j++) r_8tab_[j]= (r_8)sph(j);
  DVList dvl;
  dvl.SetS("PIXTYPE","HEALPIX ");
  dvl["ORDERING"]=sph.TypeOfMap();
  dvl["NSIDE"] = (int_4)sph.SizeIndex();
  dvl["FIRSTPIX"]=(int_4)0;
  dvl["LASTPIX"]=(int_4)(npixels-1);
  char* typeOfContent="TEMPERATURE";
  char* extname="SIMULATION";
  char* comment1="     Sky Map Pixelisation Specific Keywords";
  planck_write_bntbl(filename, npixels, typeOfContent, extname, comment1, dvl);
  delete[] r_8tab_;  r_8tab_ = NULL;

}

void FitsIoServer::save(LocalMap<double>& locm, char filename[])   
{
  int nbrows = locm.Size_x();
  int nbcols = locm.Size_y();
  long naxis = 2; 
  cout << " nombre de pts en x : " << nbrows << " en y  " << nbcols << endl;
  FITS_tab_typ_ = TDOUBLE;
  if (r_8tab_ != NULL ) { delete[] r_8tab_;  r_8tab_ = NULL; }
  r_8tab_=new r_8[nbrows*nbcols];

  int ij=0;
  for (int j=0; j< nbcols; j++) 
    for (int i = 0; i < nbrows; i++) r_8tab_[ij++]= (r_8)locm(i,j);

  DVList dvl;
  dvl["NSIDE"] = (int_4) locm.SizeIndex();
  dvl["ORDERING"]=locm.TypeOfMap();
  double theta0;
  double phi0;
  int x0;
  int y0;
  double angle;
  locm.Origin(theta0,phi0,x0,y0,angle);
  double anglex;
  double angley;
  locm.Aperture(anglex,angley);
  dvl["THETA0"] = theta0;
  dvl["PHI0"] = phi0;
  dvl["X0"] = (int_4)x0;
  dvl["Y0"] = (int_4)y0;
  dvl["ANGLE"] = angle;
  dvl["ANGLEX"] = anglex;
  dvl["ANGLEY"] = angley;
  planck_write_img(filename, naxis, nbrows, nbcols, 0, dvl);
  delete[] r_8tab_;  r_8tab_ = NULL;
}


void FitsIoServer::save(const ImageR4& DpcImg,char flnm[])
{
  long naxis=2;
  int siz_x = DpcImg.XSize();
  int siz_y = DpcImg.YSize();
  FITS_tab_typ_ = TFLOAT;



  // write FITS image
  DVList dvl;

  dvl["DATATYPE"] = (int_4)DpcImg.PixelType();
  dvl["ORG_X"] = DpcImg.XOrg();
  dvl["ORG_Y"] = DpcImg.YOrg();
  dvl["PIXSZ_X"] = DpcImg.XPxSize();
  dvl["PIXSZ_Y"] = DpcImg.YPxSize();
  dvl["IDENT"] = DpcImg.Ident();

  //dvl["NAME"] = DpcImg.Nom();

  // j utilise la methode SetS parce que ses parametres sont const et
  // que l'argument DpcImg est const dans la presente method.
  // (dans dvlist, l'operateur surcharge [] renvoie  MuTyV&, ce qui 
  // est non const et provoque un warning sur mac (CodeWarrior)
  dvl.SetS("NAME",DpcImg.Nom());

  dvl["NBSAT"] = DpcImg.nbSat;
  dvl["NBNUL"] = DpcImg.nbNul;
  dvl["MINPIX"] = DpcImg.minPix;
  dvl["MAXPIX"] = DpcImg.maxPix;
  dvl["MOYPIX"] = DpcImg.moyPix;
  dvl["SIGPIX"] = DpcImg.sigPix;
  dvl["FOND"] = DpcImg.fond;
  dvl["SIGFON"] = DpcImg.sigmaFond;

  // get the values of the DpcImage
  if (r_4tab_ != NULL) { delete [] r_4tab_;  r_4tab_ = NULL; }
  r_4tab_=new r_4[siz_x*siz_y];
  PBaseDataTypes  dataT=DataType((r_4)0);
  memcpy( r_4tab_, DpcImg.ImagePtr(), siz_x*siz_y*DataSize(dataT));
  planck_write_img(flnm, naxis, siz_x, siz_y, 0, dvl);


   delete [] r_4tab_;  r_4tab_ = NULL;
}


void FitsIoServer::save(const ImageI4& DpcImg,char flnm[])
{
  long naxis=2;
  int siz_x = DpcImg.XSize();
  int siz_y = DpcImg.YSize();
  FITS_tab_typ_ = TINT;

  // pointer to the FITS file, defined in fitsio.h
  //fitsfile *fptr;

  // initialize status before calling fitsio routines
  //int status = 0;         

  // delete old file if it already exists
  //remove(flnm);

  // create new FITS file  
  //fits_create_file(&fptr, flnm, &status); 
  //if( status )  printerror( status );


  // write FITS image
  DVList dvl;

  dvl["DATATYPE"] = (int_4)DpcImg.PixelType();
  dvl["ORG_X"] = DpcImg.XOrg();
  dvl["ORG_Y"] = DpcImg.YOrg();
  dvl["PIXSZ_X"] = DpcImg.XPxSize();
  dvl["PIXSZ_Y"] = DpcImg.YPxSize();
  dvl["IDENT"] = DpcImg.Ident();

  //dvl["NAME"] = DpcImg.Nom();
  // j utilise la methode SetS parce que ses parametres sont const et
  // que l'argument DpcImg est const dans la presente method.
  // (dans dvlist, l'operateur surcharge [] renvoie  MuTyV&, ce qui 
  // est non const et provoque un warning sur mac (CodeWarrior)
  dvl.SetS("NAME",DpcImg.Nom());

  dvl["NBSAT"] = DpcImg.nbSat;
  dvl["NBNUL"] = DpcImg.nbNul;
  dvl["MINPIX"] = DpcImg.minPix;
  dvl["MAXPIX"] = DpcImg.maxPix;
  dvl["MOYPIX"] = DpcImg.moyPix;
  dvl["SIGPIX"] = DpcImg.sigPix;
  dvl["FOND"] = DpcImg.fond;
  dvl["SIGFON"] = DpcImg.sigmaFond;

  // get the values of the DpcImage
  if (i_4tab_ != NULL) { delete [] i_4tab_;  i_4tab_ = NULL; }
  i_4tab_=new int_4[siz_x*siz_y];
  PBaseDataTypes  dataT=DataType((int_4)0);
  memcpy( i_4tab_, DpcImg.ImagePtr(), siz_x*siz_y*DataSize(dataT));


  planck_write_img(flnm, naxis, siz_x, siz_y, 0, dvl);

  // close the file
  //fits_close_file(fptr, &status);
  //if( status )  printerror( status );

   delete [] i_4tab_;  i_4tab_ = NULL;
}




void FitsIoServer::planck_write_img(char flnm[], int naxis,int n1, int n2, int n3, DVList& dvl)
{


  // pointer to the FITS file, defined in fitsio.h
  fitsfile *fptr;

  // initialize status before calling fitsio routines
  int status = 0;         

  // delete old file if it already exists
  remove(flnm);

  // create new FITS file  
  fits_create_file(&fptr, flnm, &status); 
  if( status )  printerror( status );

  int bitpix=0;
  if ( FITS_tab_typ_ == TDOUBLE) bitpix = DOUBLE_IMG;
  if ( FITS_tab_typ_ == TFLOAT)  bitpix = FLOAT_IMG;
  if ( FITS_tab_typ_ == TINT)   bitpix = LONG_IMG;
  long naxes[3];
  naxes[0] = n1;
  naxes[1] = n2;
  naxes[2] = n3;
  if (n2 > 0 && naxis < 2) cout << "  FitsIoServer:: n2 = " << n2 << " seems to be incompatible with naxis = " << naxis << endl;
  if (n3 > 0 && naxis < 3) cout << "  FitsIoServer:: n3 = " << n3 << " seems to be incompatible with naxis = " << naxis << endl;
  fits_create_img(fptr, bitpix, naxis, naxes, &status);
  if( status )  printerror( status );



  long nelements= naxes[0];
  if (naxis >=2) nelements*=naxes[1];
  if (naxis == 3)  nelements*=naxes[2];
  switch (FITS_tab_typ_) 
    {
    case TDOUBLE :
      fits_write_img(fptr, TDOUBLE, 1, nelements, r_8tab_,  &status); 
      if( status ) printerror( status );
      break;
    case TFLOAT :
      fits_write_img(fptr, TFLOAT, 1, nelements, r_4tab_,  &status); 
      if( status ) printerror( status );
      break;
    case TINT :
      fits_write_img(fptr, TINT, 1, nelements, i_4tab_,  &status); 
      if( status ) printerror( status );
      break;
    default :
      cout << " FitsIOServer : type de tableau non traite en ecriture " << endl;
      break;
    }

  // write the current date
  fits_write_date(fptr, &status);
  if( status )  printerror( status );
  // on convient d ecrire apres la date, les mots cles utilisateur.
  // la date servira de repere pour relire ces mots cles.

  //dvl.Print();
  char keyname[16];
  int flen_keyword = 9;
  char comment[FLEN_COMMENT];
  char strval[128];

  DVList::ValList::const_iterator it;
  for(it = dvl.Begin(); it != dvl.End(); it++)  
    {
      int datatype=  key_type_PL2FITS( (*it).second.typ);
      // FLEN_KEYWORD est la longueur max d'un mot-cle. Il doit y avoir une 
      // erreur dans la librairie fits qui donne FLEN_KEYWORD=72 
      // contrairement a la notice qui donne FLEN_KEYWORD=9 (ch. 5, p.39)
      //$CHECK$ Reza 20/11/99
      //      strncpy(keyname, (*it).first.substr(0,64).c_str(),flen_keyword-1);
      strncpy(keyname, (*it).first.substr(0,64).c_str(),flen_keyword);  
      keyname[flen_keyword] = '\0';
      int ival=0;
      double dval=0.;
      switch (datatype)
        {
        case TINT :
          ival=(*it).second.mtv.iv;
          strcpy(comment,"I entier");
//DBG          cerr << " Writing I " << (string)keyname << " = " << ival << endl;
          fits_write_key(fptr, datatype, keyname, &ival, comment, &status);
          break;
        case TDOUBLE :
          dval=(*it).second.mtv.dv;
          strcpy(comment,"D double");
//DBG          cerr << " Writing D " << (string)keyname << " = " << dval << endl;
          fits_write_key(fptr, datatype, keyname, &dval, comment, &status);
          break;
        case TSTRING :
          strncpy(strval, (*it).second.mtv.strv, 128); strval[127] = '\0';
          strcpy(comment,"S character string");
//DBG          cerr << " Writing S " << (string)keyname << " = " << (string)strval << endl;
          fits_write_key(fptr, datatype, keyname, &strval, comment, &status);
          break;
        default :
          cout << " FitsIOServer : probleme dans type mot cle optionnel" << endl;
          break;
        }
      if( status ) printerror( status );
    }

  strncpy(keyname, "CREATOR",flen_keyword); keyname[flen_keyword] = '\0';
  strcpy(strval, "SOPHYA");
  strcpy(comment," SOPHYA Package - FITSIOServer ");
  fits_write_key(fptr, TSTRING, keyname, &strval, comment, &status);
  if( status ) printerror( status );

  fits_write_comment(fptr,"..............................................", &status);
  fits_write_comment(fptr, " SOPHYA package - FITSIOSever ", &status);
  fits_write_comment(fptr, " (C) LAL/IN2P3-CNRS Orsay,  FRANCE 1999", &status);
  fits_write_comment(fptr, " (C) DAPNIA/CEA     Saclay, FRANCE 1999", &status);
  fits_write_comment(fptr,"..............................................", &status);
  if( status ) printerror( status );

  // close the file
  fits_close_file(fptr, &status);
  if( status )  printerror( status );
}

void FitsIoServer::planck_write_bntbl(char flnm[], int npixels, char typeOfContent[], char extname[], char comment1[], DVList& dvl)
{


  // pointer to the FITS file, defined in fitsio.h
  fitsfile *fptr;

  // initialize status before calling fitsio routines
  int status = 0;         

  // delete old file if it already exists
  remove(flnm);

  // create new FITS file  
  fits_create_file(&fptr, flnm, &status);
  //DBG  cerr << " DBG - Apres fits_create_file status = " << status << endl;
  if( status )  printerror( status );

  // primary header
  int bitpix=LONG_IMG;
  int naxis=0;
  fits_create_img(fptr, bitpix, naxis, NULL, &status);
  // write the current date
  fits_write_date(fptr, &status);
  //DBG  cerr << " DBG - Apres write_date status = " << status << endl;
  if( status )  printerror( status );


  long nrows=npixels/1024;
  if (nrows==0) nrows=1;
  if (npixels%1024 !=0) nrows++;
  int tfields=1;
  char **ttype, **tform;
  ttype= new char*[tfields];
  tform= new char*[tfields];
  char* format;
  if ( FITS_tab_typ_ == TDOUBLE) format="1024D";
  if ( FITS_tab_typ_ == TFLOAT)  format="1024E";
  if ( FITS_tab_typ_ == TINT)    format="1024I";
  for(int i = 0; i < tfields; i++)
    {
      ttype[i]= new char[FLEN_VALUE];
      strcpy(ttype[i],  typeOfContent);
      tform[i]= new char[FLEN_VALUE];
      strcpy(tform[i], format);
    }
  // 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.
 fits_create_tbl(fptr,BINARY_TBL,nrows,tfields,ttype,tform,
                 NULL,extname,&status);
  //DBG  cerr << " DBG - Apres create_tbl status = " << status << endl;
  if( status )  printerror( status );
  for( int ii = 0; ii < tfields; ii++)
    {
      delete [] ttype[ii];
      delete [] tform[ii];
    }
  delete [] ttype;
  delete [] tform;
  long nelements= npixels;
  switch (FITS_tab_typ_) 
    {
    case TDOUBLE :
      fits_write_col(fptr, TDOUBLE, 1, 1, 1, nelements, r_8tab_, &status);
      if( status ) 
        printerror( status, "planck_write_bntbl: Error writing double table" );
      break;

    case TFLOAT :
 //DBG!!! $CHECK$ Reza      for (int kk=0; kk<10; kk++) cout << r_4tab_[kk] << endl;
      fits_write_col(fptr, TFLOAT, 1, 1, 1, nelements, r_4tab_, &status);
      if( status ) 
        printerror( status, "planck_write_bntbl: Error writing float table" );
      break;
    case TINT :
      fits_write_col(fptr, TINT, 1, 1, 1, nelements, i_4tab_, &status);
      if( status ) 
        printerror( status, "planck_write_bntbl: Error writing int table");
      break;
    default :
      cout << " FitsIOServer : type de tableau non traite en ecriture " << endl;
      break;
    }
  // write supplementary keywords
  fits_write_comment(fptr, " ", &status);
  fits_write_comment(fptr, " ", &status);
  //DBG  cerr << " DBG - Apres write_comment1 status = " << status << endl;  
  if( status ) printerror( status );
  fits_write_comment(fptr,"--------------------------------------------", &status);
  fits_write_comment(fptr, comment1, &status);
  fits_write_comment(fptr,"--------------------------------------------", &status);
  //DBG  cerr << " DBG - Apres write_comment2 status = " << status << endl;  
  if( status ) printerror( status );


  int flen_keyword = 9;
  // FLEN_KEYWORD est la longueur max d'un mot-cle. Il doit y avoir une 
  // erreur dans la librairie fits qui donne FLEN_KEYWORD=72 
  // contrairement a la notice qui donne FLEN_KEYWORD=9 (ch. 5, p.39)
  DVList::ValList::const_iterator it;
  for(it = dvl.Begin(); it != dvl.End(); it++)  
    {
      int datatype=  key_type_PL2FITS( (*it).second.typ);
      char keyname[16];
      strncpy(keyname, (*it).first.substr(0,64).c_str(),flen_keyword);
      keyname[flen_keyword] = '\0';
      char comment[FLEN_COMMENT];
      int ival=0;
      double dval=0.;
      char strval[128];
      switch (datatype)
        {
        case TINT :
          ival=(*it).second.mtv.iv;
          strcpy(comment," ");
//DBG          cerr << " Writing I " << (string)keyname << " = " << ival << endl;
          fits_write_key(fptr, datatype, keyname, &ival, comment, &status);
          break;
        case TDOUBLE :
          dval=(*it).second.mtv.dv;
          strcpy(comment," ");
//DBG          cerr << " Writing D " << (string)keyname << " = " << dval << endl;
          fits_write_key(fptr, datatype, keyname, &dval, comment, &status);
          break;
        case TSTRING :
          strncpy(strval, (*it).second.mtv.strv, 128); strval[127] = '\0';
          strcpy(comment," ");
//DBG          cerr << " Writing S " << (string)keyname << " = " << (string)strval << endl;
          fits_write_key(fptr, datatype, keyname, &strval, comment, &status);
          break;
        default :
          cout << " FitsIOServer : probleme dans type mot cle optionnel" << endl;
          break;
        }
      if( status ) printerror( status );
    }
  char keyname[16];
  char strval[32];
  char comment[FLEN_COMMENT];
  strncpy(keyname, "CREATOR",flen_keyword); keyname[flen_keyword] = '\0';
  strcpy(strval, "SOPHYA");
  strcpy(comment," SOPHYA Package - FITSIOServer ");
  fits_write_key(fptr, TSTRING, keyname, &strval, comment, &status);
  if( status ) printerror( status );
  // close the file
  fits_close_file(fptr, &status);
  //DBG  cerr << " DBG - Apres close status = " << status << endl;  
  if( status )  printerror( status );
}

void FitsIoServer::planck_read_img(char flnm[], long& naxis,int& n1, int& n2, int& n3, DVList& dvl)
{
  int status=0;
  long bitpix;
  long naxes[3]={0,0,0};
  char* comment=NULL;

  // pointer to the FITS file, defined in fitsio.h
  fitsfile *fptr;
  // initialize status before calling fitsio routines
  fits_open_file(&fptr, flnm, READONLY, &status);
  if( status )  printerror( status );


  fits_read_key_lng(fptr, "BITPIX", &bitpix, comment, &status);
  if( status )  printerror( status );
  fits_read_key_lng(fptr, "NAXIS", &naxis, comment, &status);
  if( status )  printerror( status );
  int nfound;
  int nkeys=(int)naxis;
  fits_read_keys_lng(fptr, "NAXIS", 1, nkeys, naxes, &nfound, &status); 
  if( status )  printerror( status );

  n1 = naxes[0] ;
  n2 = naxes[1] ;
  n3 = naxes[2] ;


  long nelements= naxes[0];
  if (naxis >=2) nelements*=naxes[1];
  if (naxis == 3)  nelements*=naxes[2];
  int anynull=0;
  r_8 dnullval=DOUBLENULLVALUE;
  r_4 fnullval=FLOATNULLVALUE;
  int_4 inullval=0;
  // on laisse a fits le soin de convertir le type du tableau lu vers
  // le type suppose par l'utilisateur de fitsioserver 
  // 
  switch ( FITS_tab_typ_) 
    {
    case TDOUBLE :
      if (bitpix !=  DOUBLE_IMG) 
        {
          cout << " FitsIOServer : the data type on fits file " << flnm << " is not double, "
               << " conversion to double will be achieved by cfitsio lib " << endl; 
        }
      if (r_8tab_ != NULL) { delete [] r_8tab_;  r_8tab_ = NULL; }
      r_8tab_=new r_8[nelements];
      fits_read_img(fptr, TDOUBLE, 1, nelements,  &dnullval, r_8tab_, 
                    &anynull, &status);
      if( status )  printerror( status );
      break;
    case TFLOAT :
      if (bitpix !=  FLOAT_IMG) 
        {
          cout << " FitsIOServer : the data type on fits file " << flnm << " is not float, "
               << " conversion to float will be achieved by cfitsio lib " << endl; 
        }
      if (r_4tab_ != NULL) { delete [] r_4tab_;  r_4tab_ = NULL; } 
      r_4tab_=new r_4[nelements];
      fits_read_img(fptr, TFLOAT, 1, nelements,  &fnullval, r_4tab_, 
                    &anynull, &status);
      //SHV: remove useless print...
      //      for (int kk=0; kk<nelements; kk++) cout << r_4tab_[kk] << endl;
      cout << " anynull= " << anynull << endl;
      if( status )  printerror( status );
      break;


    case TINT :
      if (bitpix !=  LONG_IMG) 
        {
          cout << " FitsIOServer : the data type on fits file " << flnm << " is not long, "
               << " conversion to long will be achieved by cfitsio lib " << endl; 
        }
      if (i_4tab_ != NULL) { delete [] i_4tab_;  i_4tab_ = NULL; }
      i_4tab_=new int_4[nelements];
      fits_read_img(fptr, TINT, 1, nelements,  &inullval, i_4tab_, 
                    &anynull, &status);
      if( status )  printerror( status );
      break;


    default :
      cout << " FitsIOServer::read_img : type non traite: " << FITS_tab_typ_  << endl;
      break;
    }
  status = 0;
  char card[FLEN_CARD];
  char keyname[LEN_KEYWORD]= "";
  char strval[FLEN_VALUE];
  char comment1[FLEN_COMMENT];
  char dtype;
  char * comkey = "COMMENT";
  int keypos;
   fits_get_hdrpos(fptr,&nkeys,&keypos,&status);
   if( status )  printerror( status );
   //   cout << " planck_read_img : nombre de mots-cles : " << nkeys << endl;
   for(int j = 1; j <= nkeys; j++)
     {
       fits_read_keyn(fptr,j,card,strval,comment1,&status);
       strncpy(keyname,card,LEN_KEYWORD-1);
        if ( strncmp(keyname,comkey ,LEN_KEYWORD-1) != 0)
          {
            fits_get_keytype(strval,&dtype,&status);
            //          cout<<" keyname= "<< keyname <<" dtype= "<<dtype <<endl; 
            strip (keyname, 'B',' '); 
            strip(strval, 'B',' ');
            strip(strval, 'B','\'');
            switch( dtype ) 
              {
              case 'C':
                dvl[keyname]= strval;     
                break;
              case 'I':
                int ival;
                ctoi(strval,&ival);
                dvl[keyname]= (int_4)ival;
                break;
              case 'L':
                int ilog;
                if (strncmp(strval,"T" ,1) ==0) ilog=1;
                else ilog=0;
                dvl[keyname]= (int_4)ilog;
                break;
              case 'F':
                double dval;
                ctof(strval,&dval);
                dvl[keyname]= dval;
                break;
              default :
                cout << " FitsIOServer : type de mot-cle bizarre: " << keyname  <<" dtype= "<<dtype << endl;
                break;
              }
          }
     }

  // close the file
  status=0;
  fits_close_file(fptr, &status);
  if( status )  printerror( status );
}


void FitsIoServer::planck_read_bntbl(char flnm[], int hdunum, int& npixels, DVList& dvl)
{
  int status=0;
  int nkeys,keypos;
  int hdutype;
  int tfields;
  int datype;
  long lastpix;
  long repeat, width;
  long nrows;
  long extend;
  char* comment=NULL;

  // pointer to the FITS file, defined in fitsio.h
  fitsfile *fptr;
  // initialize status before calling fitsio routines
  fits_open_file(&fptr, flnm, READONLY, &status);
  if( status )  printerror( status );
  fits_read_key_lng(fptr, "EXTEND", &extend, comment, &status);
  if( status )  printerror( status );
  if (extend!=1) 
    {
      cout << "FitsIoServer:: le fichier fits ne contient pas d'extension binary table" << endl;
      throw IOExc("FitsIoServer::planck_read_bntbl(" + (string)flnm + ") Error No bin table extension !");
//      return;
    }
  fits_movabs_hdu(fptr, hdunum,&hdutype,&status);
  if( status )  printerror( status );
  if (hdutype!=BINARY_TBL) 
    {
      cout << "FitsIoServer:: this HDU is not a binary table " << endl;
      throw IOExc("FitsIoServer::planck_read_bntbl(" + (string)flnm + ") Error Not a bin table (1) !");
//      exit(status);
    }
   char xtension[FLEN_VALUE];
   fits_read_key_str(fptr,"XTENSION",xtension,NULL,&status);
   if( status )  printerror( status );

  char * binta = "BINTABLE";
  if (  strncmp(xtension, binta,8) != 0) 
    //   if (xtension !="BINTABLE")
     {
      cout << "FitsIoServer::  not a binary table " << endl;
      throw IOExc("FitsIoServer::planck_read_bntbl(" + (string)flnm + ") Error Not a bin table (2) !");
//      exit(status);
     }
   fits_get_hdrpos(fptr,&nkeys,&keypos,&status);
   if( status )  printerror( status );
   //cout << " nombre de mots-cles : " << nkeys << endl;
  fits_get_num_cols(fptr, &tfields, &status);
  if (tfields != 1)
    {
      cout << "FitsIoServer:: il y a plus d'une colonne" << endl;
      throw IOExc("FitsIoServer::planck_read_bntbl(" + (string)flnm + ") Error >1 column !");
//      return;
    }
  fits_get_num_rows(fptr, &nrows, &status);
  //cout << "nblignes= " << nrows << endl;
  fits_get_coltype(fptr, 1, &datype, &repeat, &width, &status);
  if( status )  printerror( status );
  //cout << " type de donnees : " << datype << endl;
  //cout << " repeat : " << repeat << endl;
  //cout << " width : " << width << endl;
  fits_read_key_lng(fptr, "LASTPIX", &lastpix, comment, &status);
  int IsLASTPIX=1;
  if( status ) 
    { 
      IsLASTPIX=0;
      printerror( status," mot cle LASTPIX" );
    }

  long nelements= nrows*repeat;
  npixels=nelements;
  if (IsLASTPIX == 1) 
    {
      if (nelements!=lastpix+1) 
        {
          cout << " probleme sur longueur du vecteur " << endl;
          cout << " nb elements prevu = " << nelements << " lastpix+1=" << lastpix+1 << endl;
          npixels=lastpix+1;
        }
    }
  int anynull;
  r_8 dnullval=DOUBLENULLVALUE;
  r_4 fnullval=FLOATNULLVALUE;
  int_4 inullval=0;
  // on laisse a fits le soin de convertir le type du tableau lu vers
  // le type suppose par l'utilisateur de fitsioserver 
  // 
  switch ( FITS_tab_typ_) 
    {
    case TDOUBLE :
      if (datype !=  TDOUBLE) 
        {
          cout << " FitsIOServer : the data type on fits file " << flnm << " is not double, "
               << " conversion to double will be achieved by cfitsio lib " << endl; 
        }
      if (r_8tab_ != NULL) { delete [] r_8tab_;  r_8tab_ = NULL; }
      r_8tab_=new r_8[ npixels];
      fits_read_col(fptr, TDOUBLE, 1, 1, 1, npixels,  &dnullval, 
                    r_8tab_, 
                    &anynull, &status);
      if( status )  printerror( status, "probleme dans lecture du tableau de doubles" );
      break;
    case TFLOAT :
      if (datype !=  TFLOAT) 
        {

          cout << " FitsIOServer : the data type on fits file " << flnm << " is not float, "
               << " conversion to float will be achieved by cfitsio lib " << endl; 
        }
      if (r_4tab_ != NULL) { delete [] r_4tab_;  r_4tab_ = NULL; }
      r_4tab_=new r_4[npixels];
      fits_read_col(fptr, TFLOAT, 1, 1, 1, npixels,  &fnullval, 
                    r_4tab_, &anynull, &status);
      if( status )  printerror( status,"probleme dans lecture du tableau de floats" );
      break;


    case TINT :
      if (datype !=  TLONG) 
        {
          cout << " FitsIOServer : the data type on fits file " << flnm << " is not long, "
               << " conversion to long will be achieved by cfitsio lib " << endl; 
        }
      if (i_4tab_ != NULL) { delete [] i_4tab_;  i_4tab_ = NULL; }
      i_4tab_=new int_4[npixels];
      fits_read_col(fptr, TLONG, 1, 1, 1, npixels,  &inullval, 
                    i_4tab_, &anynull, &status);
      //fits_read_img(fptr, TINT, 1, nelements,  &inullval, i_4tab_, 
      //                    &anynull, &status);
      if( status )  printerror( status,"probleme dans lecture du tableau de ints" );
      break;


    default :
      cout << " FitsIOServer::read_bntbl : type non traite: " << FITS_tab_typ_  << endl;
      break;
    }
  char card[FLEN_CARD];  
  char keyname[LEN_KEYWORD]= "";
  char strval[FLEN_VALUE];
  char comment1[FLEN_COMMENT];
  char dtype;
  //char bidon[LEN_KEYWORD];
  char * comkey = "COMMENT";

   for(int j = 1; j <= nkeys; j++)
     {
       //       fits_read_record(fptr, j, card, &status);
       // strncpy(keyname,card,LEN_KEYWORD-1);
       // cout << " bidon= " << keyname << endl;
       // if ( strncmp(keyname,comkey ,LEN_KEYWORD-1) != 0)
       fits_read_keyn(fptr,j,card,strval,comment1,&status);
       strncpy(keyname,card,LEN_KEYWORD-1);

        if ( strncmp(keyname,comkey ,LEN_KEYWORD-1) != 0)
          {
            fits_get_keytype(strval,&dtype,&status);
            //  cout<<" keyname= "<< keyname <<" dtype= "<<dtype <<endl; 
            strip (keyname, 'B',' '); 
            strip(strval, 'B',' ');
            strip(strval, 'B','\'');
            switch( dtype )  
              { 
              case 'C':
                dvl[keyname]= strval;     
                break;
              case 'I':
                int ival;
                ctoi(strval,&ival);
                dvl[keyname]= (int_4)ival;
                break;
              case 'L':
                int ilog;
                if (strncmp(strval,"T" ,1) ==0) ilog=1;
                else ilog=0;
                dvl[keyname]= (int_4)ilog;
                break;
              case 'F':
                double dval;
                ctof(strval,&dval);
                dvl[keyname]= dval;
                break;
              default :
                cout << " FitsIOServer : type de mot-cle bizarre: "  << keyname <<" dtype= "<<dtype << endl;
                break;
              }
          }
     }

   // close the file
   status=0;
   fits_close_file(fptr, &status);
   if( status )  printerror( status );
}


// projects a SphericalMap<double>, according sinus-method, and saves onto 
// a FITS-file
void FitsIoServer::sinus_picture_projection(SphericalMap<double>& sph, char filename[])
{
    
  long naxes[2]={600, 300};
  float* map =new float[ 600*300 ];
  int npixels=  naxes[0]*naxes[1];

   cout << " image FITS en projection SINUS" << endl;
  // table will have npixels rows
  int i,j;
  for(j=0; j < npixels; j++) map[j]=0.;

  for(j=0; j<naxes[1]; j++) 
      {
        double yd = (j+0.5)/naxes[1]-0.5;
        double theta = (0.5-yd)*Pi;
        double facteur=1./sin(theta);
        for(i=0; i<naxes[0]; i++)  
          {
            double xa = (i+0.5)/naxes[0]-0.5;
            double phi = 2.*Pi*xa*facteur+Pi;
            float th=float(theta);
            float ff=float(phi);
            if (phi<2*Pi && phi>= 0) 
              {
                map[j*naxes[0]+i] = sph.PixValSph(th, ff);
              }
          }
      }
 
    write_picture(naxes, map,  filename);
    delete [] map;
}

// projects a SphericalMap<double>, according sinus-method, and saves onto 
// a FITS-file
void FitsIoServer::sinus_picture_projection(SphericalMap<float>& sph, char filename[])
{
  // le code de cete methode duplique celui de la precedente, la seule 
  //difference etant le type de sphere en entree. Ces methodes de projection
  // sont provisoires, et ne servent que pour les tests. C est pourquoi je 
  // ne me suis pas casse la tete, pour l instant
    
  long naxes[2]={600, 300};
  float* map = new float[ 600*300 ];
  int npixels=  naxes[0]*naxes[1];

   cout << " image FITS en projection SINUS" << endl;
  // table will have npixels rows
  int i,j;
  for(j=0; j < npixels; j++) map[j]=0.;
  for(j=0; j<naxes[1]; j++) 
      {
        double yd = (j+0.5)/naxes[1]-0.5;
        double theta = (0.5-yd)*Pi;
        double facteur=1./sin(theta);
        for(i=0; i<naxes[0]; i++)  
          {
            double xa = (i+0.5)/naxes[0]-0.5;
            double phi = 2.*Pi*xa*facteur+Pi;
            float th=float(theta);
            float ff=float(phi);
            if (phi<2*Pi && phi>= 0) 
              {
                map[j*naxes[0]+i] = sph.PixValSph(th, ff);
              }
          }
      }
 
    write_picture(naxes, map,  filename);
    delete [] map;

}

// projects a SphericalMap<float>, according Mollweide-method, and saves onto 
// a FITS-file
void FitsIoServer::Mollweide_picture_projection(SphericalMap<float>& sph, char filename[])
{
  // le code de cete methode duplique celui de la precedente, la seule 
  //difference etant le type de sphere en entree. Ces methodes de projection
  // sont provisoires, et ne servent que pour les tests. C est pourquoi je 
  // ne me suis pas casse la tete, pour l instant
    
  long naxes[2]={600, 300};
  float* map = new float[ 600*300 ];
  int npixels=  naxes[0]*naxes[1];

   cout << " image FITS en projection MOLLWEIDE" << endl;
  // table will have npixels rows
   int i,j;
   for(j=0; j < npixels; j++) map[j]=0.;
   for(j=0; j<naxes[1]; j++) 
     {
       double yd = (j+0.5)/naxes[1]-0.5;
       double facteur=2.*Pi/sin(acos(yd*2));
       double theta = (0.5-yd)*Pi;
       for(i=0; i<naxes[0]; i++)  
         {
           double xa = (i+0.5)/naxes[0]-0.5;
           double phi =  xa*facteur+Pi;
           float th=float(theta);
           float ff=float(phi);
           if (phi<2*Pi && phi>= 0) 
             {
               map[j*naxes[0]+i] = sph.PixValSph(th, ff);
             }
         }
     }
   
   write_picture(naxes, map,  filename);
   delete [] map;
   
}
// projects a SphericalMap<double>, according Mollweide-method, and saves onto 
// a FITS-file
void FitsIoServer::Mollweide_picture_projection(SphericalMap<double>& sph, char filename[])
{
  // le code de cete methode duplique celui de la precedente, la seule 
  //difference etant le type de sphere en entree. Ces methodes de projection
  // sont provisoires, et ne servent que pour les tests. C est pourquoi je 
  // ne me suis pas casse la tete, pour l instant
    
  long naxes[2]={600, 300};
  float* map = new float[ 600*300 ];
  int npixels=  naxes[0]*naxes[1];

   cout << " image FITS en projection MOLLWEIDE" << endl;
  // table will have npixels rows
  int i,j;
   for(j=0; j < npixels; j++) map[j]=0.;
   for(j=0; j<naxes[1]; j++) 
     {
       double yd = (j+0.5)/naxes[1]-0.5;
       double facteur=2.*Pi/sin(acos(yd*2));
       double theta = (0.5-yd)*Pi;
       for(i=0; i<naxes[0]; i++)  
         {
           double xa = (i+0.5)/naxes[0]-0.5;
           double phi =  xa*facteur+Pi;
           float th=float(theta);
           float ff=float(phi);
           if (phi<2*Pi && phi>= 0) 
             {
               map[j*naxes[0]+i] = sph.PixValSph(th, ff);
             }
         }
     }
   
   write_picture(naxes, map,  filename);
   delete [] map;
   
}



// saves a (LocalMap<double> on a FITS-file in order to be visualized
// (for tests)
void FitsIoServer::picture(LocalMap<double>& lcm, char filename[])
{
    
  long naxes[2];
  naxes[0] = lcm.Size_x();
  naxes[1] = lcm.Size_y();
  int npixels=  naxes[0]*naxes[1];
  float* map = new float[npixels];

  // table will have npixels rows
  int i,j;
  for(j=0; j < npixels; j++) map[j]=0.;
    for(j=0; j<naxes[1]; j++) 
      {
        for(i=0; i<naxes[0]; i++)  
          {
                map[j*naxes[0]+i] = lcm(i, j);
          }
      }
 
    write_picture(naxes, map,  filename);
    delete [] map;
}



void  FitsIoServer::write_picture(long* naxes, float* map, char* filename) const
{
 
  int bitpix = FLOAT_IMG;
  long naxis = 2;

  //pointer to the FITS file, defined in fitsio.h 
  fitsfile *fptr; 
  // delete old file if it already exists
  remove(filename);
    // initialize status before calling fitsio routines
  int status = 0;         

  // create new FITS file  
  fits_create_file(&fptr, filename, &status); 
  if( status )  printerror( status );

  // write the required header keywords
  fits_create_img(fptr, bitpix, naxis, naxes, &status);
  if( status )  printerror( status );

  // write the current date
  fits_write_date(fptr, &status);
  if( status )  printerror( status );


  // first row in table to write 
  //  long firstrow  = 1;  Not referenced 
  // first element in row   
  long firstelem = 1;  
  //  int colnum = 1;     Not referenced 
  int nelements=naxes[0]*naxes[1];
  fits_write_img(fptr, TFLOAT, firstelem, nelements, map, &status);
  if( status )  printerror( status );

  // close the file
  fits_close_file(fptr, &status);
  if( status )  printerror( status );
}


bool FitsIoServer::check_keyword(fitsfile *fptr,int nkeys,char keyword[])

  //*****************************************************/
  //* check if the specified keyword exits in the CHU   */
  //*****************************************************/
{

  bool KEY_EXIST = false;
  int status = 0;
  char strbide[FLEN_VALUE];
  char keybide[LEN_KEYWORD]= "";
  for(int jj = 1; jj <= nkeys; jj++)
    {
      if( fits_read_keyn(fptr,jj,keybide,strbide,NULL,&status) )
        printerror( status );
      if( !strcmp(keybide,keyword) )
        { 
          KEY_EXIST= true;
          break;
        }
    }
  return(KEY_EXIST);
}

void FitsIoServer::readheader ( char filename[] )

  //**********************************************************************/
  //* Print out all the header keywords in all extensions of a FITS file */
  //**********************************************************************/
{

  // standard string lengths defined in fitsioc.h
  char card[FLEN_CARD];

  // pointer to the FITS file, defined in fitsio.h
  fitsfile *fptr;

  int status = 0;
  if ( fits_open_file(&fptr, filename, READONLY, &status) ) 
    printerror( status );
  
  // attempt to move to next HDU, until we get an EOF error
  int hdutype;
  for (int ii = 1; !(fits_movabs_hdu(fptr,ii,&hdutype,&status));ii++) 
    {
      if (hdutype == ASCII_TBL)
        printf("\nReading ASCII table in HDU %d:\n", ii);
      else if (hdutype == BINARY_TBL)
        printf("\nReading binary table in HDU %d:\n", ii);
      else if (hdutype == IMAGE_HDU)
        printf("\nReading FITS image in HDU %d:\n", ii);
      else
        {
          printf("Error: unknown type of this HDU \n");
          printerror( status );
        }

      // get the number of keywords
      int nkeys, keypos;
      if ( fits_get_hdrpos(fptr, &nkeys, &keypos, &status) )
        printerror( status );
      
      printf("Header listing for HDU #%d:\n", ii);
      for (int jj = 1; jj <= nkeys; jj++)  
        {
          if ( fits_read_record(fptr, jj, card, &status) )
            printerror( status );

          // print the keyword card
          printf("%s\n", card); 
        }
      printf("END\n\n");
    }

  // got the expected EOF error; reset = 0
  if (status == END_OF_FILE)   
    status = 0;             
  else
    printerror( status );

  if ( fits_close_file(fptr, &status) )
    printerror( status );
  
  return;
}

void FitsIoServer::printerror(int& status) const

  //*****************************************************/
  //* Print out cfitsio error messages and exit program */
  //*****************************************************/
{

  // print out cfitsio error messages and exit program 
      // print error report  
      fits_report_error(stderr, status);
      // terminate the program, returning error status 
      //      exit( status ); 
      status=0;
}
void FitsIoServer::printerror(int& status, char* texte) const

  //*****************************************************/
  //* 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;
      status=0;
}