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

#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>;
  ownobj_= true;

  Read(inputfile,hdunum);
  //  dobj_->SetTemp(true);
}


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>::SetDataObj(AnyDataObj & o)
{
  SphereHEALPix<T> * po = dynamic_cast< SphereHEALPix<T> * >(&o);
  if (po == NULL) return;
  if (ownobj_ && dobj_) delete dobj_;
  dobj_ = po; 
  ownobj_ = false;
} 




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

  DVList dvl;  

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

  dvl["PIXTYPE"] = "HEALPIX";
  dvl.SetComment("PIXTYPE", "HEALPIX Pixelization");
  dvl["ORDERING"]= dobj_->TypeOfMap();
  dvl.SetComment("ORDERING", "Pixel ordering scheme, either RING or NESTED");
  
  int nSide= dobj_->SizeIndex();
  dvl["NSIDE"]= (int_4) nSide;
  dvl.SetComment("NSIDE","Resolution parameter for HEALPIX" );

  int nPix= dobj_->NbPixels();
  dvl["FIRSTPIX"]= (int_4) 0;
  dvl.SetComment("FIRSTPIX", "First pixel # (0 based)");
  dvl["LASTPIX"]= (int_4) nPix-1;
  dvl.SetComment("LASTPIX", "Last pixel # (0 based)");
  dvl["Content"]= "SphereHEALPix";
  dvl["EXTNAME"]= "SIMULATION";
  dvl["COMMENT"]= "Sky Map Pixelisation Specific Keywords";
   
  // On ecrit les dataBlocks 
  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;
  os.makeHeaderBntblOnFits(Type, Noms, nPix, 1, dvl, extname, dummy);
  delete [] Noms[0];
  delete [] Noms;
  os.putColToFits(0,  nPix, dobj_->pixels_.Data()); 
} 

template <class T> 
void FITS_SphereHEALPix<T>::ReadFromFits(FitsInFile& is)
{
  if(dobj_ == NULL) 
    {
      dobj_= new SphereHEALPix<T>;
      ownobj_= true;      
    }

 
  int nbcols, nbentries;
  nbcols = is.NbColsFromFits();
  if (nbcols != 1) 
    {
      throw IOExc("le fichier fits n'est pas une sphere Healpix");
    }
  DVList dvl=is.DVListFromFits();
   nbentries = is.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";
  string coordsys= dvl.GetS("COORDSYS");
  // $CHECK$ - Reza 2/05/2000 - compare(size_t, size_t,...) passe pas sur g++     
  //  if(coordsys.compare(0,7,"ROTATION",0,7) == 0)  $CHECK$ 
  if(coordsys.substr(0,8) == "ROTATION")
    {
      id= SphereCoordSys_ROTATION;
      description= "ROTATION SphereCoordSystem";
    } 
  //  else if(coordsys.compare(0,4,"OTHER",0,4) == 0) $CHECK$ Reza 2/05/2000 
  else if(coordsys.substr(0,5)  == "OTHER" )
    {
      id= SphereCoordSys_OTHER;
      description= "OTHER SphereCoordSystem";
    }
  dobj_->SetCoordSys(new SphereCoordSys(id,description));   

  string ordering= dvl.GetS("ORDERING");
  //  if(ordering.compare(0,3,"RING",0,3) != 0)  $CHECK$ Reza 2/05/2000
  if(ordering.substr(0,4) != "RING" )
    { 
      cerr << "FITS_SphereHEALPix/Error Not supported ORDERING= " 
           << ordering << " substr(0,4)=[" << ordering.substr(0,4) << "]" << 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);
  is.GetSingleColumn(dobj_->pixels_.Data(),nPix);

  // on effectue le decoupage en tranches 
  dobj_->SetThetaSlices();
  dobj_->Info()=is.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