source: Sophya/trunk/SophyaExt/FitsIOServer/fitsspherehealpix.cc@ 903

#include "pexceptions.h"
#include "fitsspherehealpix.h"




///////////////////////////////////////////////////////////////////////
// ----objets delegues  pour la gestion de persistance I/O format fits--
// SphereHealpix
////////////////////////////////////////////////////////////////////


template <class T>
FITS_SphereHEALPix<T>::FITS_SphereHEALPix()
{
  dobj_= new SphereHEALPix<T>;
  ownobj= true;
}
  
template <class T>
FITS_SphereHEALPix<T>::FITS_SphereHEALPix(char inputfile[],int hdunum)
{
  dobj_= new SphereHEALPix<T>;
  dobj_->SetTemp(true);
  ownobj= true;

  ReadF(inputfile,hdunum);
}

template <class T>
void FITS_SphereHEALPix<T>::Write(char outputfile[],int hdunum)
{
  WriteF(outputfile, hdunum);
}

template <class T>
FITS_SphereHEALPix<T>::FITS_SphereHEALPix(const SphereHEALPix<T>& obj)
{
  dobj_= new SphereHEALPix<T>(obj, true);
  dobj_->SetTemp(true);
  ownobj= true;
}

template <class T>
FITS_SphereHEALPix<T>::FITS_SphereHEALPix(SphereHEALPix<T>* obj)
{
  dobj_= obj;
  ownobj= false;
}

template <class T>
FITS_SphereHEALPix<T>::~FITS_SphereHEALPix()
{
  if (ownobj && dobj_) delete dobj_;
}

template <class T>
AnyDataObj* FITS_SphereHEALPix<T>::DataObj()
{
  return(dobj_);
}

template <class T>
void FITS_SphereHEALPix<T>::WriteToFits(const FitsFile& fn)
{
  if(dobj_ == NULL) 
    { 
      cout << " WriteTo:: dobj_= null " << endl;
      return;
    }

  DVList dvl;  

  SphereCoordSys* cs= dobj_->GetCoordSys();
  string description= cs->description();
  dvl["COORDSYS"]= description;

  dvl["ORDERING"]= dobj_->TypeOfMap();
  
  int nSide= dobj_->SizeIndex();
  dvl["NSIDE"]= nSide;

  int nPix= dobj_->NbPixels();
  dvl["FIRSTPIX"]= 0;
  dvl["LASTPIX"]= nPix-1;
  dvl["Content"]= "SphereHEALPix";
  dvl["EXTNAME"]= "SIMULATION";
  dvl["COMMENT"]= "Sky Map Pixelisation Specific Keywords";
   
  // On ecrit les dataBlocks 
  //  fn.put_fits_bintab(dobj_->pixels_,dvl);
  char** Noms = new char*[1];   
  Noms[0]= new char[15];
  strncpy(Noms[0],dvl.GetS("Content").c_str(),15);
  char extname[15];
  strncpy(extname,dvl.GetS("EXTNAME").c_str(),15);

  char Type[2];
  if (typeid(T) == typeid(r_8) ) Type[0]='D';
  else 
    if (typeid(T) == typeid(r_4) )  Type[0]='E';
    else
      {
        cout <<  " type de la sphere= " <<  typeid(T).name() << endl;
        throw IOExc("FITS_SphereHEALPix:: unknown type");
      }
  Type[1]='\0'; 
  vector<int> dummy;
  fn.makeHeaderBntblOnFits(Type, Noms, nPix, 1, dvl, extname, dummy);
  delete [] Noms[0];
  delete [] Noms;
  putColToFits(0,  nPix, dobj_->pixels_.Data()); 
} 

template <class T> 
void FITS_SphereHEALPix<T>::ReadFromFits(const FitsFile& fn)
{
  if(dobj_ == NULL) 
    {
      dobj_= new SphereHEALPix<T>;
      dobj_->SetTemp(true);
      ownobj= true;      
    }

 
  int nbcols, nbentries;
  nbcols = fn.NbColsFromFits();
  if (nbcols != 1) 
    {
      throw IOExc("le fichier fits n'est pas une sphere Healpix");
    }
  //   const DVList* dvl = &fn.DVListFromFits();
  DVList dvl=fn.DVListFromFits();
  //  dvl.Print();
   nbentries = fn.NentriesFromFits(0);
  int lastpix=dvl.GetI("LASTPIX");
  if (lastpix>0)
    {
      cout << " lastpix trouve= " << lastpix << endl;
      if (nbentries!=lastpix+1) 
        {
          cout << " nb pixels from LASTPIX = "  << lastpix+1 << endl;
          nbentries=lastpix+1;
        }
    }
  // Let's Read the SphereCoordSys object
  int id= SphereCoordSys_NEUTRAL;
  string description= "NEUTRAL SphereCoordSystem";
  const char* coordsys= dvl.GetS("COORDSYS").c_str();    
  if(strncmp(coordsys,"ROTATION",8) == 0)
    {
      id= SphereCoordSys_ROTATION;
      description= "ROTATION SphereCoordSystem";
    } 
  else if(strncmp(coordsys,"OTHER",5) == 0)
    {
      id= SphereCoordSys_OTHER;
      description= "OTHER SphereCoordSystem";
    }
  dobj_->SetCoordSys(new SphereCoordSys(id,description));   

  const char* ordering= dvl.GetS("ORDERING").c_str();
  if(strncmp(ordering,"RING",4) != 0)
    { 
      cout << " ORDERING= " << ordering << endl;
      throw IOExc(" FITS_SphereHEALPix:: numerotation non RING");
    }

  int nside= dvl.GetI("NSIDE");
  if(nside <= 0)
    throw IOExc("FITS_SphereHEALPix:: No resol parameter nside");
  int nPix = 12*nside*nside;
  if (nPix != nbentries)
    {
      cout << "WARNING: le nombre de pixels relu " << nbentries << " est incoherent avec la pixelisation Healpix qui donne nPix= " << nPix << endl; 
      }
  double Omega= 4.0*Pi/nPix;
  dobj_->setParameters(nside,nPix,Omega);
  // On lit les DataBlocks;
  dobj_->pixels_.ReSize(nPix);
  fn.GetSingleColumn(dobj_->pixels_.Data(),nPix);

  // on effectue le decoupage en tranches 
  dobj_->SetThetaSlices();
  dobj_->Info()=fn.DVListFromFits();


}


#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template FITS_SphereHEALPix<r_8>
#pragma define_template FITS_SphereHEALPix<r_4>
#endif
#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
template class FITS_SphereHEALPix<r_8>;
template class FITS_SphereHEALPix<r_4>;
#endif