#include "nbmath.h"
#include <math.h>
#include <iostream.h>
#include "fitsioserver.h" 
#include "dvlist.h"
#include "strutil.h"

  typedef map<string, MuTyV, less<string> >  ValList; 
void FitsIoServer::load(TMatrix<double>& mat, char flnm[]) 
{
  int nbrows=0;
  int nbcols=0;
  FITS_tab_typ_ = TDOUBLE;
  int status=0;
  long naxis;
  int n1, n2, n3;
  DVList dvl;
  // pointer to the FITS file, defined in fitsio.h 
  fitsfile *fptr;
  fits_open_file(&fptr, flnm, READONLY, &status);
  if( status )  printerror( status );
  planck_read_img(fptr, naxis, n1, n2, n3, dvl);
  // close the file
  fits_close_file(fptr, &status);
  if( status )  printerror( status );

  nbrows=n1;
  nbcols=n2;
  if (naxis == 1) nbcols=1;
  // cout << " lecture  du dvlist par load " << endl;
  dvl.Print();
  ValList::const_iterator it;
  for(it = dvl.Begin(); it != dvl.End(); it++)  
    {
      char datatype= (*it).second.typ;
      char keyname[]="";
      strcpy(keyname, (*it).first.substr(0,64).c_str());   
      int ival=0;
      double dval=0.;
      char strval[]="";
      switch (datatype)
	{
	case 'I' :
	  ival=(*it).second.mtv.iv;
	  break;
	case 'D' :
	  dval=(*it).second.mtv.dv;
	  break;
	case 'S' :
	  strcpy(strval, (*it).second.mtv.strv); 
	  break;
	default :
	  cout << " probleme dans type mot cle optionnel" << endl;
	  break;
	}
    }

  if (naxis > 2) 
    {
      cout << " 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) = dtab_[ij++];
  // cout << " chargement termine " << endl;
}

void FitsIoServer::load(SphericalMap<double>& sph, char flnm[]) 
{
  int npixels=0;
  int nside=0;
  FITS_tab_typ_ = TDOUBLE;
  read_sphere(flnm, npixels, 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;
	  exit(0);
	}
      if (nside != sph.SizeIndex())
	{ 
	  sph.Resize(nside);
	  cout << " resize the sphere to nside=  " << nside  << endl;
	}
      else
	{
	  cout << " FITSIOSERVER : same resolution, surprising ! " << endl;
	  exit(0);
	}
    }
  for (int j = 0; j < sph.NbPixels(); j++) sph(j)= dtab_[j];
  // cout << " chargement termine " << endl;
}
void FitsIoServer::load(SphericalMap<float>& sph, char flnm[]) 
{
  int npixels=0;
  int nside=0;
  FITS_tab_typ_ = TFLOAT;
  read_sphere(flnm, npixels, 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;
	  exit(0);
	}
      if (nside != sph.SizeIndex())
	{ 
	  sph.Resize(nside);
	  cout << " resize the sphere to nside=  " << nside  << endl;
	}
      else
	{
	  cout << " FITSIOSERVER : same resolution, surprising ! " << endl;
	  exit(0);
	}
    }
  for (int j = 0; j < sph.NbPixels(); j++) sph(j)= ftab_[j];
  //  cout << " chargement termine " << endl;
}

void FitsIoServer::load(LocalMap<double>& lcm, char flnm[]) 
{
  int nside;
  int nbrows=0;
  int nbcols=0;
  FITS_tab_typ_ = TDOUBLE;
  int status=0;
  long naxis;
  int n1, n2, n3;
  DVList dvl;
  // pointer to the FITS file, defined in fitsio.h 
  fitsfile *fptr;
  fits_open_file(&fptr, flnm, READONLY, &status);
  if( status )  printerror( status );
  planck_read_img(fptr, naxis, n1, n2, n3, dvl);
  // close the file
  fits_close_file(fptr, &status);
  if( status )  printerror( status );

  nbrows=n1;
  nbcols=n2;
  if (naxis != 2) 
    {
      cout << " le fichier fits n'est pas une localmap, naxis= " << naxis << endl;
    }
  // cout << " lecture  du dvlist par load " << endl;
  dvl.Print();
  ValList::const_iterator it;
  for(it = dvl.Begin(); it != dvl.End(); it++)  
    {
      char datatype= (*it).second.typ;
      char keyname[]="";
      strcpy(keyname, (*it).first.substr(0,64).c_str());   
      int ival=0;
      double dval=0.;
      char strval[]="";
      switch (datatype)
	{
	case 'I' :
	  ival=(*it).second.mtv.iv;
	  break;
	case 'D' :
	  dval=(*it).second.mtv.dv;
	  break;
	case 'S' :
	  strcpy(strval, (*it).second.mtv.strv); 
	  break;
	default :
	  cout << " probleme dans type mot cle optionnel" << endl;
	  break;
	}
      if (!strcmp(keyname, "NSIDE") )
	{
	  // cout << " j'ai trouve NSIDE " << endl;
	  nside = ival;
	} 
      //if (!strcmp(keyname, "ORDERING") )
      //	{
      //	  cout << " j'ai trouve ORDERING " << 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");
  //cout << " theta0= " << theta0 << endl;
  //cout << " phi0= " << phi0 << endl;
  //cout << " x0= " << x0 << endl;
  //cout << " y0= " << y0 << endl;
  //cout << " angle= " << angle << endl;
  //cout << " anglex= " << anglex << endl;
  //cout << " angley= " << angley << endl;

  // 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 << " resize 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) = dtab_[ij++];
  //cout << " chargement termine " << endl;
}

void FitsIoServer::printerror(int status) const
{
  //Print out cfitsio error messages and exit program 
  if( status )
    {
      //print error report  
      fits_report_error(stderr, status);
      //terminate the program, returning error status 
      exit( status ); 
    }
  return;
}
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 (dtab_ != NULL ) delete[] dtab_;
  dtab_=new double[nbrows*nbcols];


  //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(filename);

  // create new FITS file  
  fits_create_file(&fptr, filename, &status); 
  if( status )  printerror( status );
  int ij=0;
  for (int j=0; j< nbcols; j++) 
    for (int i = 0; i < nbrows; i++) dtab_[ij++]= mat(i,j);

  DVList dvl;

  //  cout << " save  : debut ecriture " << endl;
  planck_write_img(fptr, naxis, nbrows, nbcols, 0, dvl);
  //  cout << " save  : fin ecriture " << endl;


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

}

//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 (dtab_ != NULL ) delete[] dtab_;
  dtab_=new double[npixels];

  //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(filename);

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

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

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

  // decider ulterieurement ce qu'on fait de ce qui suit, specifique
  // pour l'instant, aux spheres gorski. Actuellement, on ne sauve pas 
  // les informations d'analyse harmonique
  /*
  // write supplementary keywords
  fits_write_comment(fptr, " ", &status);
  if( status ) printerror( status );

  char comment[] = "Generated by synfast";
  char svalue [] = "SYNFAST ";
  strcat(svalue, version);
  fits_write_key(fptr, TSTRING, "ORIG-MAP", svalue, comment, &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,"Maximum multipole l used in map synthesis");
  int lmax= sph.nlmax();
  fits_write_key(fptr, TINT, "MAX-LPOL", &lmax, 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,"FWHM in DEGREES");
  float fwhm_deg= sph.fwhm()/60.0;
  fits_write_key(fptr, TFLOAT, "FWHM", &fwhm_deg, 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 );
  
  //   write information about the power spectrum
  char card[]="";
  sph.powfile(card);
  strcpy(comment,"Input power spectrum in : ");
  strcat(comment, card);
  fits_write_comment(fptr, comment, &status);
  if( status ) printerror( status );

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

  strcpy(comment,"C(2)= ");
  sprintf(card,"%g",sph.quadrupole());
  strcat(comment, card);
  fits_write_comment(fptr, comment, &status);
  if( status ) printerror( status );

  strcpy(comment,"keywords relative to input power spectrum ");
  fits_write_comment(fptr, comment, &status);
  if( status ) printerror( status );
  strcpy(comment,"have not yet been defined for PLANCK Surveyor Project");
  fits_write_comment(fptr, comment, &status);
  if( status ) printerror( status );
  */

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

}

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


  //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(filename);

  // create new FITS file  
  fits_create_file(&fptr, filename, &status); 
  if( status )  printerror( status );
  for (int j = 0; j < npixels; j++) ftab_[j]= sph(j);
  DVList dvl;
  dvl["NSIDE"] = sph.SizeIndex();
  dvl["ORDERING"]=sph.TypeOfMap();

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


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

}

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 (dtab_ != NULL ) delete[] dtab_;
  dtab_=new double[nbrows*nbcols];

  //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(filename);

  // create new FITS file  
  fits_create_file(&fptr, filename, &status); 
  if( status )  printerror( status );
  int ij=0;
  for (int j=0; j< nbcols; j++) 
    for (int i = 0; i < nbrows; i++) dtab_[ij++]= locm(i,j);

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


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

}


void FitsIoServer::planck_write_img( fitsfile *fptr, int naxis,int n1, int n2, int n3, DVList& dvl)
{

  int status=0;
  int bitpix=0;
  if ( FITS_tab_typ_ == TDOUBLE) bitpix = DOUBLE_IMG;
  if ( FITS_tab_typ_ == TFLOAT)  bitpix = FLOAT_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, dtab_,  &status); 
      if( status ) printerror( status );
      break;
    case TFLOAT :
      fits_write_img(fptr, TFLOAT, 1, nelements, ftab_,  &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();
  ValList::const_iterator it;
  for(it = dvl.Begin(); it != dvl.End(); it++)  
    {
      int datatype=  key_type_PL2FITS( (*it).second.typ);
      char keyname[]="";
      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)
      strncpy(keyname, (*it).first.substr(0,64).c_str(),flen_keyword-1);
      char comment[FLEN_COMMENT]="";
      int ival=0;
      double dval=0.;
      char strval[]="";
      switch (datatype)
	{
	case TINT :
	  ival=(*it).second.mtv.iv;
	  strcpy(comment,"I entier");
	  fits_write_key(fptr, datatype, keyname, &ival, comment, &status);
	  break;
	case TDOUBLE :
	  dval=(*it).second.mtv.dv;
	  strcpy(comment,"D double");
	  fits_write_key(fptr, datatype, keyname, &dval, comment, &status);
	  break;
	case TSTRING :
	  strcpy(strval, (*it).second.mtv.strv);
	  strcpy(comment,"S character string");
	  fits_write_key(fptr, datatype, keyname, &strval, comment, &status);
	  break;
	default :
	  cout << " probleme dans type mot cle optionnel" << endl;
	  break;
	}
      if( status ) printerror( status );
    }
}
void FitsIoServer::planck_read_img( fitsfile *fptr, long& naxis,int& n1, int& n2, int& n3, DVList& dvl)
{
  int status=0;
  long bitpix;
  long naxes[3]={0,0,0};
  char* comment=NULL;

  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;
  double dnullval=0.;
  float fnullval=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 (dtab_ != NULL) delete [] dtab_;
      dtab_=new double[nelements];
      fits_read_img(fptr, TDOUBLE, 1, nelements,  &dnullval, dtab_, 
		    &anynull, &status);
      if( status )  printerror( status );
      break;
    case TFLOAT :
      if (ftab_ != NULL) delete [] ftab_;
      ftab_=new float[nelements];
      fits_read_img(fptr, TFLOAT, 1, nelements,  &fnullval, ftab_, 
		    &anynull, &status);
      if( status )  printerror( status );
      break;
    default :
      cout << " fitsio::read_img : type non traite: " << FITS_tab_typ_  << endl;
      break;
    }
  status = 0;
  char card[FLEN_CARD];
  int num = 0;
  char comment2[FLEN_COMMENT] = "x";
  char keyname[]= "";
  char datekey[]= "DATE";
  char endkey[] = "END";
  char typ='x';
  int ival;
  double dval;
  char strval[]="";
  // on a convenu que les mots cles utilisateur sont apres le mot cle DATE
  // on va jusqu'au mot cle DATE
  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)
  while (status == 0 &&  strncmp(keyname, datekey,4) != 0 ) 
    {
      num++;
      fits_read_record(fptr, num, card, &status);
      strncpy(keyname,card,flen_keyword-1);
    }
  if (status != 0 ) 
    {
      cout << " fitsio::planck_read_img : erreur, mot cle DATE absent " << endl;
    }
  // on recupere la liste des mots-cles utilisateurs
  while (status == 0) 
    {
      num++;
      // on lit un record pour recuperer le nom du mot-cle
      fits_read_record(fptr, num, card, &status);
      strncpy(keyname,card,flen_keyword-1);
      char value[FLEN_VALUE];
      // on recupere le premier caractere du commentaire, qui contient 
      // le code du type de la valeur
      // (tant que ce n est pas le mot cle END)
      fits_read_keyword(fptr, keyname, value, comment2, &status);
      if (  strncmp(keyname, endkey,flen_keyword-1) != 0) 
	{
	  typ = comment2[0];
	  // quand le type est connu, on lit la valeur
	  switch (typ) 
	    {
	    case 'I' :
	      fits_read_key(fptr, TINT, keyname, &ival, comment2, &status);
	      if( status )  printerror( status );
	      strip (keyname, 'B',' '); 
	      dvl[keyname] = ival;
	      break;
	    case 'D' : 
	      fits_read_key(fptr, TDOUBLE, keyname, &dval, comment2, &status);
	      if( status )  printerror( status );
	      strip (keyname, 'B',' '); 
	      dvl[keyname] = dval;
	      break;
	    case 'S' : 
	      fits_read_key(fptr, TSTRING, keyname, strval, comment2, &status);
	      if( status )  printerror( status );
	      strip (keyname, 'B',' '); 
	      strip(strval, 'B',' ');
	      dvl[keyname]=strval;
	      break;
	    default :
	      cout << " FITSIOSERVER::planck_read_img : type de donnee non prevu " << endl;
	      break;
	    }
	}
    }

}


void FitsIoServer::sinus_picture_projection(SphericalMap<double>& sph, char filename[])
{
    
  long naxes[2]={600, 300};
  float map [ 600*300 ];
  int npixels=  naxes[0]*naxes[1];

   cout << " image FITS en projection SINUS" << endl;
  // table will have npixels rows
  for(int j=0; j < npixels; j++) map[j]=0.;
    for(int 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(int 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);

}
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 [ 600*300 ];
  int npixels=  naxes[0]*naxes[1];

   cout << " image FITS en projection SINUS" << endl;
  // table will have npixels rows
  for(int j=0; j < npixels; j++) map[j]=0.;
    for(int 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(int 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);


}

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
  for(int j=0; j < npixels; j++) map[j]=0.;
    for(int j=0; j<naxes[1]; j++) 
      {
	for(int i=0; i<naxes[0]; i++)  
	  {
		map[j*naxes[0]+i] = lcm(i, j);
	  }
      }
 
    write_picture(naxes, map,  filename);
    delete [] map;
}

void  FitsIoServer::read_sphere(char flnm[], int& npixels, int& nside)
{
  int status=0;
  long naxis;
  nside=0;
  int n1, n2, n3;
  DVList dvl;
  // pointer to the FITS file, defined in fitsio.h 
  fitsfile *fptr;
  fits_open_file(&fptr, flnm, READONLY, &status);
  if( status )  printerror( status );
  planck_read_img(fptr, naxis, n1, n2, n3, dvl);
  // close the file
  fits_close_file(fptr, &status);
  if( status )  printerror( status );

  npixels=n1;
  dvl.Print();
  ValList::const_iterator it;
  for(it = dvl.Begin(); it != dvl.End(); it++)  
    {
      char datatype= (*it).second.typ;
      char keyname[]="";
      strcpy(keyname, (*it).first.substr(0,64).c_str());   
      int ival=0;
      double dval=0.;
      char strval[]="";
      switch (datatype)
	{
	case 'I' :
	  ival=(*it).second.mtv.iv;
	  break;
	case 'D' :
	  dval=(*it).second.mtv.dv;
	  break;
	case 'S' :
	  strcpy(strval, (*it).second.mtv.strv); 
	  break;
	default :
	  cout << " probleme dans type mot cle optionnel" << endl;
	  break;
	}
      if (!strcmp(keyname, "NSIDE") )
	{
	  //  cout << " j'ai trouve NSIDE " << endl;
	  nside = ival;
	} 
      // if (!strcmp(keyname, "ORDERING") )
      ///	{
      //	  cout << " j'ai trouve ORDERING " << endl;
      //	} 
    }

  if (naxis != 1) 
    {
      cout << " le fichier fits n'est pas une sphere, naxis= " << naxis << endl;
    }
}


void  FitsIoServer::write_picture(long* naxes, float* map, char* filename) const
{
  //pointer to the FITS file, defined in fitsio.h 
  fitsfile *fptr;  
  int bitpix = FLOAT_IMG;
  long naxis = 2;   

  // 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;
  // first element in row   
  long firstelem = 1;  
  int colnum = 1;
  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 );
}