source: Sophya/trunk/SophyaExt/FitsIOServer/fitsspherethetaphi.cc@ 2103

#include "pexceptions.h"
#include "fitsspherethetaphi.h"
#include "tarray.h"
#include "fitstarray.h" 



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

/*!
  \class SOPHYA::FITS_SphereThetaPhi
  \ingroup FitsIOServer
  FITS format I/O handler for SOPHYA::SphereThetaPhi objects.
*/

template <class T>
FITS_SphereThetaPhi<T>::FITS_SphereThetaPhi()
{
  dobj_= new SphereThetaPhi<T>;
  ownobj_= true;
}

template <class T>
FITS_SphereThetaPhi<T>::FITS_SphereThetaPhi(char inputfile[],int hdunum)
{
    dobj_= new SphereThetaPhi<T>;
    ownobj_= true;
  Read(inputfile,hdunum);
}

template <class T>
FITS_SphereThetaPhi<T>::FITS_SphereThetaPhi(const SphereThetaPhi<T>& obj)
{
   dobj_= new SphereThetaPhi<T>(obj);
  ownobj_= true;
}

template <class T>
FITS_SphereThetaPhi<T>::FITS_SphereThetaPhi(SphereThetaPhi<T>* obj)
{
  dobj_= obj;
  ownobj_= false;
}

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

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

template <class T>
void FITS_SphereThetaPhi<T>::SetDataObj(AnyDataObj & o)
{
  SphereThetaPhi<T> * po = dynamic_cast< SphereThetaPhi<T> * >(&o);
  if (po == NULL) return;
  if (ownobj_ && dobj_) delete dobj_;
  dobj_ = po; 
  ownobj_ = false;
} 

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

  DVList dvl( dobj_->Info() );  

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

  dvl["PIXTYPE"] = "TETAFI";
  dvl.SetComment("PIXTYPE", "Theta-Phi-LAL Pixelization");
  dvl["ORDERING"]= dobj_->TypeOfMap();
  dvl.SetComment("ORDERING", " (unrelevant) ");

  int IndTheta= dobj_->SizeIndex();
  dvl["INDTHMAX"]= (int_4) IndTheta;
  dvl.SetComment("INDTHMAX","Resolution parameter for Theta-Phi pixelization sch." );

  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"]= "SphereThetaPhi";
  dvl.SetComment("Content", "name of SOPHYA object");
  
  // On ecrit les dataBlocks 
  vector<string> Noms(4);   
  Noms[0] = dvl.GetS("Content");
  Noms[1] = "NphiParBande";
  Noms[2] = "CumulPixParBande";
  Noms[3] = "ThetaBande";
  string extname("SIMULATION");


  string Type;
  if (typeid(T) == typeid(r_8) ) Type+='D';
  else if (typeid(T) == typeid(r_4) )  Type+='E';
  else if (typeid(T) == typeid(int_4) )  Type+='J';
    else
      {
        cout <<  "FITS_SphereThetaPhi<T>::WriteToFits - Unsupported data type (" 
             <<  typeid(T).name() << ") throwing IOExc " << endl;
        throw IOExc("FITS_SphereThetaPhi<T>::WriteToFits() unsupported data type");
      }

  Type += 'J';
  Type += 'J';
  Type += 'D';

  vector<int> dummy;
  os.makeHeaderBntblOnFits(Type, Noms, nPix, 4, &dvl, extname, dummy);
  os.PutColToFits(0,  nPix, dobj_->pixels_.Data()); 
  os.PutColToFits(1,  IndTheta+1, dobj_->NPhi_.Data()); 
  os.PutColToFits(2,  IndTheta+1, dobj_->TNphi_.Data()); 
  os.PutColToFits(3,  IndTheta+1, dobj_->Theta_.Data()); 
}

template <class T> 
void FITS_SphereThetaPhi<T>::ReadFromFits(FitsInFile& is)
{

  if(dobj_ == NULL) 
    {
      dobj_= new SphereThetaPhi<T>;
      ownobj_= true;      
    }
  int nbcols, nbentries;

  nbcols = is.NbColsFromFits();
  if (nbcols != 4) 
    {
      throw IOExc("le fichier fits n'est pas une sphere ThetaPhi");
    }
  DVList dvl=is.DVListFromFits();
  nbentries = is.NentriesFromFits(0);
  int lastpix=dvl.GetI("LASTPIX");
  if (lastpix>0)
    {
      if (nbentries!=lastpix+1) 
        {
          nbentries=lastpix+1;
        }
    }
  
  // Let's Read the SphereCoordSys object
  int id= SphereCoordSys_NEUTRAL;
  string description= "NEUTRAL SphereCoordSystem";
  string coordsys= dvl.GetS("COORDSYS");
  if(coordsys.substr(0,8) == "ROTATION")
    {
      id= SphereCoordSys_ROTATION;
      description= "ROTATION SphereCoordSystem";
    } 
  else if(coordsys.substr(0,5)  == "OTHER" )
    {
      id= SphereCoordSys_OTHER;
      description= "OTHER SphereCoordSystem";
    }
  dobj_->SetCoordSys(new SphereCoordSys(id,description));   
  
  int IndThetaMax= dvl.GetI("INDTHMAX");
  int nPix = 0;
  if(IndThetaMax <= 0)
    {
      throw IOExc("FITS_SphereTheta:: No resol parameter INDTHMAX");
    }
  nPix = nbentries;
  double Omega= 4.0*Pi/nPix;
  dobj_->setParameters(IndThetaMax,nPix,Omega);
  
  // On lit les DataBlocks;
  //
  dobj_->pixels_.ReSize(nPix);
  is.GetBinTabFCol(dobj_->pixels_.Data(),nPix,0);
  dobj_->NPhi_.ReSize(IndThetaMax+1);
  is.GetBinTabFCol(dobj_->NPhi_.Data(),IndThetaMax+1,1);
  dobj_->TNphi_.ReSize(IndThetaMax+1);
  is.GetBinTabFCol(dobj_->TNphi_.Data(),IndThetaMax+1,2);
  dobj_->Theta_.ReSize(IndThetaMax+1);
  is.GetBinTabFCol(dobj_->Theta_.Data(),IndThetaMax+1,3);
}

template <class T>
void FITS_SphereThetaPhi<T>::Mollweide_picture_projection(char filename[])
{
  int ni = 300;
  int nj = 600;
  TMatrix<float> map(ni, nj);

   int i,j;
   for(i=0; i < ni; i++) 
     {
       double yd = (i+0.5)/ni-0.5;
       double facteur=2.*Pi/sin(acos(yd*2));
       double theta = (0.5-yd)*Pi;
       for(j=0; j < nj; j++)  
         {
           double xa = (j+0.5)/nj-0.5;
           double phi =  xa*facteur+Pi;
           float th=float(theta);
           float ff=float(phi);
           if (phi<2*Pi && phi>= 0) 
             {
               map(i,j) = (float)dobj_->PixValSph(th, ff);
             }
         }
     }
   FitsOutFile  fout(filename);
   fout.firstImageOnPrimaryHeader(true);
   FITS_TArray<r_4> fta(map);
   fta.Write(fout);
   
}
template <class T>
void FITS_SphereThetaPhi<T>::sinus_picture_projection(char filename[])
{
  int ni = 300;
  int nj = 600;
  TMatrix<float> map(ni, nj);

  int i,j;
   for(i=0; i < ni; i++) 
     {
       double yd = (i+0.5)/ni-0.5;
       double theta = (0.5-yd)*Pi;
        double facteur=1./sin(theta);
       for(j=0; j < nj; j++)  
         {
           double xa = (j+0.5)/nj-0.5;
           double phi = 2.*Pi*xa*facteur+Pi;
           float th=float(theta);
           float ff=float(phi);
           if (phi<2*Pi && phi>= 0) 
             {
               map(i,j) = (float)dobj_->PixValSph(th, ff);
             }
         }
     }
   FitsOutFile  fout(filename);
   fout.firstImageOnPrimaryHeader(true);
   FITS_TArray<r_4> fta(map);
   fta.Write(fout);

   
}


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