source: Sophya/trunk/SophyaPI/PIext/nomskymapadapter.cc@ 2683

#include "sopnamsp.h"
#include "machdefs.h"
#include <stdlib.h>
#include <math.h>
#include <typeinfo>
#include <iostream>
#include <string>
#include <complex>

#include "datatype.h"

#include "nomskymapadapter.h"
#include "skymap.h"
#include "pitvmaad.h"
#include "complexios.h"

#include "fitsspherehealpix.h"
#include "fitslocalmap.h"

// Valeur par defaut pour la projection MolleWeide, en dehors de la sphere 
static double _prjmol_defval = 0.;
void Set_Project_Mol_DefVal(double def)
{
  _prjmol_defval = def;
  return;
}

//  Classe array adapter pour localMap
template <class T> 
class LocalMapArrAdapter : public P2DArrayAdapter {
public:
                LocalMapArrAdapter(LocalMap<T>* lm, bool d=false) :
                    P2DArrayAdapter(lm->SizeX(), lm->SizeY())
                    { ad = d; map = lm; }

  virtual       ~LocalMapArrAdapter() { if (ad) delete map; }
  virtual double  Value(int ix, int iy) { return((*map)(ix, iy)); }

protected :
  bool ad;
  LocalMap<T>* map;
};

/* --Methode-- */
DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
double LocalMapArrAdapter< complex<float> >::Value(int ix, int iy) 
{
double re,im;
re = (*map)(iy, ix).real();
im = (*map)(iy, ix).imag();
return(sqrt(re*re+im*im));
}
/* --Methode-- */
DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
double LocalMapArrAdapter< complex<double> >::Value(int ix, int iy) 
{
double re,im;
re = (*map)(iy, ix).real();
im = (*map)(iy, ix).imag();
return(sqrt(re*re+im*im));
}

//----------------------------------------------------------------
// Class Adaptateur d'objet (Pour NamedObjMgr) d'objet PixelMap<T>
//----------------------------------------------------------------


/* --Methode-- */
template <class T> 
NOMAdapter_PixelMap<T>::NOMAdapter_PixelMap(PixelMap<T> * o)
  : NObjMgrAdapter(o)
{
mMap = o;
}

/* --Methode-- */
template <class T> 
NOMAdapter_PixelMap<T>::~NOMAdapter_PixelMap()
{
}

/* --Methode-- */
template <class T> 
NObjMgrAdapter* NOMAdapter_PixelMap<T>::Clone(AnyDataObj* o)
{
PixelMap<T>* m = dynamic_cast<PixelMap<T> *>(o);
if (m) return ( new NOMAdapter_PixelMap<T>(m) );
return ( new NObjMgrAdapter(o) );
}

/* --Methode-- */
template <class T> 
string NOMAdapter_PixelMap<T>::GetDataObjType()
{
string type = "PixelMap< ";
LocalMap<T>* lm = dynamic_cast< LocalMap<T> * >(mMap); 
if (lm != NULL) type = "LocalMap< ";
SphereThetaPhi<T>* st = dynamic_cast< SphereThetaPhi<T> * >(mMap); 
if (st != NULL) type = "SphereThetaPhi< ";
SphereHEALPix<T>* sg = dynamic_cast< SphereHEALPix<T> * >(mMap); 
if (sg != NULL) type = "SphereHEALPix< ";
SphereECP<T>* se = dynamic_cast< SphereECP<T> * >(mMap); 
if (se != NULL) type = "SphereECP< ";

// type +=  DecodeTypeIdName(typeid(T).name());
type += DataTypeInfo<T>::getTypeName();
type +=  " > ";
return(type);
}

/* --Methode-- */
template <class T> 
AnyDataObj* NOMAdapter_PixelMap<T>::CloneDataObj(bool share)
{
LocalMap<T>* lm = dynamic_cast< LocalMap<T> * >(mMap); 
if (lm != NULL) return( new LocalMap<T>(*lm, share) );
SphereThetaPhi<T>* st = dynamic_cast< SphereThetaPhi<T> * >(mMap); 
if (st != NULL) return( new SphereThetaPhi<T>(*st, share) );
SphereHEALPix<T>* sg = dynamic_cast< SphereHEALPix<T> * >(mMap); 
if (sg != NULL) return( new SphereHEALPix<T>(*sg, share) );
SphereECP<T>* se = dynamic_cast< SphereECP<T> * >(mMap); 
if (se != NULL) return( new SphereECP<T>(*se, share) );
return(NULL);
}

/* --Methode-- */
template <class T> 
void NOMAdapter_PixelMap<T>::ReadFits(string const & flnm)
{
LocalMap<T>* lm = dynamic_cast< LocalMap<T> * >(mMap); 
if (lm != NULL) {
  FitsInFile fis(flnm);
  fis >> (*lm);
  return;
}
SphereHEALPix<T>* sg = dynamic_cast< SphereHEALPix<T> * >(mMap); 
if (sg != NULL) {
  FitsInFile fis(flnm);
  fis >> (*sg);
  return;
}
string s = typeid(*mMap).name();
cout << " NOMAdapter_PixelMap<T>::ReadFits() - Error " 
     << " Not supported for " << s << endl;
return;
}


/* --Methode-- */
template <class T> 
void NOMAdapter_PixelMap<T>::SaveFits(string const & flnm)
{
LocalMap<T>* lm = dynamic_cast< LocalMap<T> * >(mMap); 
if (lm != NULL) {
  FitsOutFile fos(flnm);
  fos << (*lm);
  return;
}
SphereHEALPix<T>* sg = dynamic_cast< SphereHEALPix<T> * >(mMap); 
if (sg != NULL) {
  FitsOutFile fos(flnm);
  fos << (*sg);
  return;
}
string s = typeid(*mMap).name();
cout << " NOMAdapter_PixelMap<T>::SaveFits() - Error " 
     << " Not supported for " << s << endl;
return;
}

// ---- Specialisation pour complexes -----
DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
void NOMAdapter_PixelMap< complex<r_4> >::ReadFits(string const & flnm)
{
cout << " void NOMAdapter_PixelMap< complex<r_4> >::ReadFits() - Error " 
     << " Not supported (complex data type)" << endl;
}
DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
void NOMAdapter_PixelMap< complex<r_4> >::SaveFits(string const & flnm)
{
cout << " void NOMAdapter_PixelMap< complex<r_4> >::SaveFits() - Error " 
     << " Not supported (complex data type)" << endl;
}
DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
void NOMAdapter_PixelMap< complex<r_8> >::ReadFits(string const & flnm)
{
cout << " void NOMAdapter_PixelMap< complex<r_4> >::ReadFits() - Error " 
     << " Not supported (complex data type)" << endl;
}
DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
void NOMAdapter_PixelMap< complex<r_8> >::SaveFits(string const & flnm)
{
cout << " void NOMAdapter_PixelMap< complex<r_8> >::SaveFits() - Error " 
     << " Not supported (complex data type)" << endl;
}


/* --Methode-- */
template <class T> 
void NOMAdapter_PixelMap<T>::SavePPF(POutPersist& pos, string const & nom)
{
LocalMap<T>* lm = dynamic_cast< LocalMap<T> * >(mMap); 
if (lm != NULL) {
  FIO_LocalMap<T> fio(lm);
  fio.Write(pos, nom);
  return;
  }
SphereThetaPhi<T>* st = dynamic_cast< SphereThetaPhi<T> * >(mMap); 
if (st != NULL) {
  FIO_SphereThetaPhi<T> fio(st);
  fio.Write(pos, nom);
  return;
  }
SphereHEALPix<T>* sg = dynamic_cast< SphereHEALPix<T> * >(mMap); 
if (sg != NULL) {
  FIO_SphereHEALPix<T> fio(sg);
  fio.Write(pos, nom);
  return;
  }
SphereECP<T>* se = dynamic_cast< SphereECP<T> * >(mMap); 
if (se != NULL) {
  FIO_SphereECP<T> fio(se);
  fio.Write(pos, nom);
  return;
  }
string s = typeid(*mMap).name();
cout << "NOMAdapter_PixelMap<T>::SavePPF() - Error : Not supported for " << s << endl;
}

/* --Methode-- */
template <class T> 
void NOMAdapter_PixelMap<T>::Print(ostream& os)
{
string s = typeid(*mMap).name();
T moy, sig;
MeanSig(moy, sig);
cout << "SkyMap Type: " << s << " NbPixels= " << mMap->NbPixels() << endl;
cout << "  Mean= " << moy << "  Sig2= " << sig << endl;
}

// -- Fonction pour convertir un X,Y en Theta-Phi sur projection MolleWeide --
char * _SphMollweideProj_XY2ThetaPhi(P2DArrayAdapter * aa, int ix, int iy);

/* --Methode-- */
template <class T> 
P2DArrayAdapter* NOMAdapter_PixelMap<T>::Get2DArray(string & opt)
{
LocalMap<T>* lm = dynamic_cast< LocalMap<T> * >(mMap); 
if (lm != NULL) return(new LocalMapArrAdapter<T>(lm, false));
int nr = 250;
int nc = 500;
SphereECP<T>* se = dynamic_cast< SphereECP<T> * >(mMap); 
size_t olen = opt.length();
if ((se != NULL) && (opt.find("ecparray") < olen)) {
  // On veut la carte bi-dim complete de la sphere ECP
  TMatrix<T> * mtxecp = new TMatrix<T>(se->GetPixelArray(), true);  // on partage les donnees
  POTMatrixAdapter<T> * potmtxecp = new POTMatrixAdapter<T>(mtxecp, true);
  if (opt.find("ecpscaledeg") < olen)  // Echelle angulaire en degre
    potmtxecp->DefineXYCoordinates(se->MinPhi()*180./M_PI, se->MinTheta()*180./M_PI, 
                                   se->DeltaPhi()*180./M_PI, se->DeltaTheta()*180./M_PI);
  else if (opt.find("ecpscalerad") < olen)  // Echelle angulaire en radian
    potmtxecp->DefineXYCoordinates(se->MinPhi(), se->MinTheta(), se->DeltaPhi(), se->DeltaTheta());

  return ( potmtxecp );
}
SphericalMap<T>* sm = dynamic_cast< SphericalMap<T> *>(mMap);
if (sm != NULL) { nr = sqrt(0.75*mMap->NbPixels());  nc = 2*nr; }
TMatrix<T> * mtx = new TMatrix<T>(nr, nc);
Project_Mol(*mtx, (T)_prjmol_defval);
POTMatrixAdapter<T> * potma = new POTMatrixAdapter<T>(mtx, true);
potma->SetInfoStringFunction(_SphMollweideProj_XY2ThetaPhi);
return ( potma );
}

/* --Methode-- */
template <class T> 
NTupleInterface* NOMAdapter_PixelMap<T>::GetNTupleInterface(bool& adel)
{
adel = true;
return( new NTupInt_PixelMap<T>(mMap) );
}

/* --Methode-- */
template <class T> 
void NOMAdapter_PixelMap<T>::MeanSig(T& gmoy, T& gsig)
{
  gmoy=0.;
  gsig = 0.;
  T valok;
  for(int k=0; k<mMap->NbPixels(); k++) {
    valok = (*mMap)(k);
    gmoy += valok;  gsig += valok*valok;
  }
  gmoy /= (T)mMap->NbPixels();
  gsig = gsig/(T)mMap->NbPixels() - gmoy*gmoy;

}

/* --Methode-- */
template <class T> 
void NOMAdapter_PixelMap<T>::Project_Mol(TMatrix<T> & mtx, T defval)
{
  r_8 xa, yd, teta,phi, facteur;
  int_4 l,c,k;
  int_4 nl = mtx.NRows();
  int_4 nc = mtx.NCols();
  mtx = defval;   // On met tout a defval
//  cout << " NRows= " << nl << "  NCols= " << nc  << endl;
  for(l=0; l<nl; l++) {
    yd = (r_8)(l+0.5)/(r_8)nl-0.5;
    facteur=2.*M_PI/sin(acos((double)yd*2));
    teta = (yd+0.5)*Pi;
    //    teta = (0.5-yd)*M_PI;
    for(c=0; c<nc; c++)  {
      xa = (r_8)(c+0.5)/(r_8)nc-0.5;  
      phi = xa*facteur+M_PI;
      if ( (phi <= 2*M_PI) && (phi >= 0.) ) {
        k = mMap->PixIndexSph(teta, phi);
        mtx(l,c) = (*mMap)(k);
      }
    }
  }
}

static char m_buff[128];
char * _SphMollweideProj_XY2ThetaPhi(P2DArrayAdapter * aa, int ix, int iy)
{
  r_8 xa, yd, teta,phi, facteur;
  int_4 nl = aa->YSize();
  int_4 nc = aa->XSize();
  int l, c;
  aa->IJ2xy(ix, iy, c, l);
  yd = (r_8)(l+0.5)/(r_8)nl-0.5;
  facteur=2.*M_PI/sin(acos((double)yd*2));
  teta = (yd+0.5)*Pi;
  xa = (r_8)(c+0.5)/(r_8)nc-0.5;  
  phi = xa*facteur+M_PI;
  if ( (phi <= 2*M_PI) && (phi >= 0.) ) 
    sprintf(m_buff,"(l,b=%5.1f,%5.1f)", phi*180./M_PI, 90.-teta*180./M_PI);
  else 
    sprintf(m_buff,"(l,b=?,?)");

  return (m_buff);
}

// -------------------------------------------------------------

/* --Methode-- */
template <class T> 
NTupInt_PixelMap<T>::NTupInt_PixelMap(PixelMap<T>* m)
{
mMap = m;
}

/* --Methode-- */
template <class T> 
NTupInt_PixelMap<T>::~NTupInt_PixelMap()
{
}

/* --Methode-- */
template <class T> 
sa_size_t NTupInt_PixelMap<T>::NbLines() const
{
return( mMap->NbPixels() );
}

/* --Methode-- */
template <class T> 
sa_size_t NTupInt_PixelMap<T>::NbColumns() const 
{
return(8);
}

/* --Methode-- */
template <class T> 
r_8* NTupInt_PixelMap<T>::GetLineD(sa_size_t n) const 
{
int i;
if ((n < 0) || (n >= (int)(mMap->NbPixels()) )) 
    for(i=0; i<8; i++)  mRet[i] = 0.;
else { 
  double teta,phi;
  mMap->PixThetaPhi(n, teta, phi);
  mRet[0] = n;  mRet[1] = mMap->PixVal(n);
  mRet[2] = mRet[1];  mRet[3] = 0.;
  mRet[4] = mRet[1];  mRet[5] = 0.;
  mRet[6] = teta;     mRet[7] = phi;
  }
return(mRet);
}

/* --Methode-- */
template <class T> 
string NTupInt_PixelMap<T>::VarList_C(const char* nx) const 
{
string nomx;
if (nx) nomx = nx;
else nomx = "_xh_";
string vardec = "double i,k,val,real,imag,mod,phas,teta,phi; \n";
vardec += "i = " + nomx + "[0]; k = " + nomx + "[0];  val = " + nomx + "[1]; \n";
vardec += "real = " + nomx + "[2];  imag = " + nomx + "[3]; \n";
vardec += "mod  = " + nomx + "[4];  phas = " + nomx + "[5]; \n";
vardec += "teta = " + nomx + "[6];  phi  = " + nomx + "[7]; \n";
return(vardec);
}


DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
r_8* NTupInt_PixelMap< complex<float> >::GetLineD(int n) const 
{
int i;
if ((n < 0) || (n >= (int)(mMap->NbPixels()) )) 
    for(i=0; i<8; i++)  mRet[i] = 0.;
else { 
  double teta,phi;
  mMap->PixThetaPhi(n, teta, phi);
  mRet[0] = n;  
  mRet[2] = mMap->PixVal(n).real();  mRet[3] = mMap->PixVal(n).imag();
  mRet[1] = mRet[4] = sqrt(mRet[2]*mRet[2]+mRet[3]*mRet[3]);
  mRet[5] = atan2(mRet[3], mRet[2]);
  mRet[6] = teta;     mRet[7] = phi;
  }
return(mRet);
}

DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
r_8* NTupInt_PixelMap< complex<double> >::GetLineD(int n) const 
{
int i;
if ((n < 0) || (n >= (int)(mMap->NbPixels()) )) 
    for(i=0; i<8; i++)  mRet[i] = 0.;
else { 
  double teta,phi;
  mMap->PixThetaPhi(n, teta, phi);
  mRet[0] = n;  
  mRet[2] = mMap->PixVal(n).real();  mRet[3] = mMap->PixVal(n).imag();
  mRet[1] = mRet[4] = sqrt(mRet[2]*mRet[2]+mRet[3]*mRet[3]);
  mRet[5] = atan2(mRet[3], mRet[2]);
  mRet[6] = teta;     mRet[7] = phi;
  }
return(mRet);
}


#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template NOMAdapter_PixelMap<int_4>
#pragma define_template NOMAdapter_PixelMap<r_4>
#pragma define_template NOMAdapter_PixelMap<r_8>
#pragma define_template NOMAdapter_PixelMap< complex<float> >
#pragma define_template NOMAdapter_PixelMap< complex<double> >
#pragma define_template NTupInt_PixelMap<int_4>
#pragma define_template NTupInt_PixelMap<r_4>
#pragma define_template NTupInt_PixelMap<r_8>
#pragma define_template NTupInt_PixelMap< complex<float> >
#pragma define_template NTupInt_PixelMap< complex<double> >
#endif

#if defined(ANSI_TEMPLATES) 
template class NOMAdapter_PixelMap<int_4>;
template class NOMAdapter_PixelMap<r_4>;
template class NOMAdapter_PixelMap<r_8>;
template class NOMAdapter_PixelMap< complex<float> >;
template class NOMAdapter_PixelMap< complex<double> >;
template class NTupInt_PixelMap<int_4>;
template class NTupInt_PixelMap<r_4>;
template class NTupInt_PixelMap<r_8>;
template class NTupInt_PixelMap< complex<float> >;
template class NTupInt_PixelMap< complex<double> >;
#endif