source: Sophya/trunk/SophyaLib/SkyMap/localmap.cc@ 1827

#include "localmap.h"
#include "smathconst.h"
#include <complex>
#include "piocmplx.h"
#include "fiondblock.h"

#include <string.h>
#include <iostream.h>


/*!
  \class SOPHYA::LocalMap

  A local map has an origin in (theta0, phi0), mapped to pixel(x0, y0)
   (x0, y0 might be outside of this local map)
 default value of (x0, y0) is middle of the map, center of pixel(nx/2, ny/2)
    A local map is a 2 dimensional array, with i as column index and j
    as row index. The map is supposed to lie on a plan tangent to the 
    celestial sphere in a point whose coordinates are (x0,y0) on the local
   map and (theta0, phi0) on the sphere. The range of the map is defined
    by two values of angles covered respectively by all the pixels in 
    x direction and all the pixels in y direction (SetSize()).

    A "reference plane" is considered : this plane is tangent to the 
    celestial sphere in a point with angles theta=Pi/2 and phi=0. This 
    point is the origine of coordinates is of the reference plane. The 
    x-axis is the tangent parallel to the equatorial line and oriented 
    toward the increasing phi's ; the y-axis is parallel to the meridian 
    line and oriented toward the north pole.
    
    Internally, a map is first defined within this reference plane and
    tranported until the point (theta0, phi0) in such a way that both 
    axes are kept parallel to meridian and parallel lines of the sphere.
    The user can define its own map with axes rotated with respect to
    reference axes (this rotation is characterized by angle between 
    the local parallel line and the wanted x-axis-- see method 
    SetOrigin(...))
*/
template<class T>
LocalMap<T>::LocalMap() :  pixels_()
{
  InitNul();
}

template<class T>
LocalMap<T>::LocalMap(int_4 nx, int_4 ny) : nSzX_(nx), nSzY_(ny)
{
  InitNul();
  nPix_= nx*ny;
  pixels_.ReSize(nPix_);
  pixels_.Reset();
}

template<class T>
LocalMap<T>::LocalMap(const LocalMap<T>& lm, bool share)
  : pixels_(lm.pixels_, share)
{

  if(lm.mInfo_) mInfo_= new DVList(*lm.mInfo_);
  recopierVariablesSimples(lm);
}

template<class T>
LocalMap<T>::LocalMap(const LocalMap<T>& lm)
  : pixels_(lm.pixels_)

{

  if(lm.mInfo_) mInfo_= new DVList(*lm.mInfo_);
  recopierVariablesSimples(lm);
}

template<class T>
LocalMap<T>::~LocalMap()
{
  InitNul();
}



/*!  \fn void SOPHYA::LocalMap::ReSize(int_4 nx, int_4 ny) 

  Resize storage area for pixels 
*/
template<class T>
void LocalMap<T>::ReSize(int_4 nx, int_4 ny) 
{
  InitNul();
  nSzX_ = nx;
  nSzY_ = ny;
  nPix_= nx*ny;
  pixels_.ReSize(nPix_);
  pixels_.Reset();
}

////////////////////////// methodes de copie/share

template<class T>
void LocalMap<T>::CloneOrShare(const LocalMap<T>& a)
{
      recopierVariablesSimples(a);     
      pixels_.CloneOrShare(a.pixels_);
      if (mInfo_) {delete mInfo_; mInfo_ = NULL;}
      if (a.mInfo_) mInfo_ = new DVList(*(a.mInfo_));
}
template<class T>
void LocalMap<T>::Share(const LocalMap<T>& a)
{
      recopierVariablesSimples(a);     
      pixels_.Share(a.pixels_);
      if (mInfo_) {delete mInfo_; mInfo_ = NULL;}
      if (a.mInfo_) mInfo_ = new DVList(*(a.mInfo_));
}

template<class T>
LocalMap<T>& LocalMap<T>::CopyElt(const LocalMap<T>& a)
{
  if (NbPixels() < 1) 
    throw RangeCheckError("LocalMap<T>::CopyElt(const LocalMap<T>& )  - Not Allocated Array ! ");
  if (NbPixels() != a.NbPixels()) 
    throw(SzMismatchError("LocalMap<T>::CopyElt(const LocalMap<T>&) SizeMismatch")) ;
      recopierVariablesSimples(a);     
  int k;
  for (k=0; k< nPix_; k++) pixels_(k) = a.pixels_(k);
  return(*this);

}

template<class T>
LocalMap<T>& LocalMap<T>::Set(const LocalMap<T>& a)
   {
  if (this != &a) 
    { 
      if (a.NbPixels() < 1) 
        throw RangeCheckError("LocalMap<T>::Set(a ) - Array a not allocated ! ");
      if (NbPixels() < 1) CloneOrShare(a);
      else CopyElt(a);
      if (mInfo_) delete mInfo_;
      mInfo_ = NULL;
      if (a.mInfo_) mInfo_ = new DVList(*(a.mInfo_));
    }
  return(*this);
   }


/*! \fn int_4 SOPHYA::LocalMap::NbPixels() const

   \return number of pixels 
*/
template<class T>
int_4 LocalMap<T>::NbPixels() const
{
  return(nPix_);
}

/*! \fn T& SOPHYA::LocalMap::PixVal(int_4 k)

   \return value of pixel with index k 
*/
template<class T>
T& LocalMap<T>::PixVal(int_4 k)
{
  if((k < 0) || (k >= nPix_)) 
    {
      cout << " LocalMap::PIxVal : exceptions  a mettre en place" <<endl;
      THROW(rangeCheckErr);
    }
  return(pixels_(k));
}

/*! \fn T const& SOPHYA::LocalMap::PixVal(int_4 k) const
  const version of previous method 
*/
template<class T>
T const& LocalMap<T>::PixVal(int_4 k) const
{
  if((k < 0) || (k >= nPix_)) 
    {
      cout << " LocalMap::PIxVal : exceptions  a mettre en place" <<endl;
      throw "LocalMap::PIxVal Pixel index out of range"; 
    } 
  return *(pixels_.Data()+k);
}

/*! \fn bool SOPHYA::LocalMap::ContainsSph(double theta, double phi) const
 \return true if teta,phi in map  
 <b> Not yet implemented </b>
*/
template<class T>
bool LocalMap<T>::ContainsSph(double theta, double phi) const
{
// $CHECK$  Reza 11/11/99 - A modifier
  throw NotAvailableOperation("LocalMap<T>::ContainsSph() - Not yet implemented !");
  return(true);
}

/*!  \fn int_4 SOPHYA::LocalMap::PixIndexSph(double theta,double phi) const

  \return index of the pixel with spherical coordinates (theta,phi) 
*/
template<class T>
int_4 LocalMap<T>::PixIndexSph(double theta,double phi) const
{
  int_4 i,j;
  if(!(originFlag_) || !(extensFlag_)) 
    {
      cout << " LocalMap: correspondance carte-sphere non etablie" << endl;
      exit(0);
    }

  // theta et phi en coordonnees relatives (on se ramene a une situation par rapport au plan de reference)
  double theta_aux= theta;
  double phi_aux  = phi;
  UserToReference(theta_aux, phi_aux);

  // coordonnees dans le plan local en unites de pixels
  double x,y;
  AngleProjToPix(theta_aux,phi_aux, x, y);

  double xmin= -x0_-0.5;
  double xmax= xmin+nSzX_;
  if((x > xmax) || (x < xmin)) return(-1);
  double xcurrent= xmin;
  for(i = 0; i < nSzX_; i++ ) 
    {
      xcurrent += 1.;
      if( x < xcurrent ) break;
    }
  double ymin= -y0_-0.5;
  double ymax= ymin+nSzY_;
  if((y >  ymax) || (y < ymin)) return(-1);
  double ycurrent= ymin;
  for(j = 0; j < nSzY_; j++ ) 
    {
      ycurrent += 1.;
      if( y < ycurrent ) break;
    }
  return (j*nSzX_+i);
}

/*! \fn void SOPHYA::LocalMap::PixThetaPhi(int_4 k,double& theta,double& phi) const

   \return (theta, phi) coordinates of pixel with index k 
*/
template<class T>
void LocalMap<T>::PixThetaPhi(int_4 k,double& theta,double& phi) const
{
  if(!(originFlag_) || !(extensFlag_)) 
    {
      cout << " LocalMap: correspondance carte-sphere non etablie" << endl;
      exit(0);
    }

  int_4 i,j;
  Getij(k,i,j);

  double X= double(i-x0_);
  double Y= double(j-y0_);
  // situation de ce pixel dans le plan de reference
  double x= X*cos_angle_-Y*sin_angle_;
  double y= X*sin_angle_+Y* cos_angle_;
  // projection sur la sphere
  PixProjToAngle(x, y, theta, phi);
  // retour au plan utilisateur
  ReferenceToUser(theta, phi);
}

/*! \fn T SOPHYA::LocalMap::SetPixels(T v)

Set all pixels to value v 
*/
template <class T>
T LocalMap<T>::SetPixels(T v)
{
pixels_.Reset(v);
return(v);
}
 
/*! \fn double SOPHYA::LocalMap::PixSolAngle(int_4 k) const

   Pixel Solid angle  (steradians)

    All the pixels have not necessarly the same size in (theta, phi)
    because of the projection scheme which is not yet fixed. 
*/  
template<class T>
double LocalMap<T>::PixSolAngle(int_4 k) const
{
  int_4 i,j;
  Getij(k,i,j);
  double X= double(i-x0_);
  double Y= double(j-y0_);
  double XR= X+double(i)*0.5;
  double XL= X-double(i)*0.5;
  double YU= Y+double(j)*0.5;
  double YL= Y-double(j)*0.5;

  // situation  dans le plan de reference
  double x0= XL*cos_angle_-YL*sin_angle_;
  double y0= XL*sin_angle_+YL*cos_angle_;
  double xa= XR*cos_angle_-YL*sin_angle_;
  double ya= XR*sin_angle_+YL*cos_angle_;
  double xb= XL*cos_angle_-YU*sin_angle_;
  double yb= XL*sin_angle_+YU*cos_angle_;

  // projection sur la sphere
  double thet0,phi0,theta,phia,thetb,phib;
  PixProjToAngle(x0, y0, thet0, phi0);
  PixProjToAngle(xa, ya, theta, phia);
  PixProjToAngle(xb, yb, thetb, phib);

  // angle solide
  double sol= fabs((xa-x0)*(yb-y0)-(xb-x0)*(ya-y0));
  return sol;
}

/*! \fn void SOPHYA::LocalMap::SetOrigin(double theta0,double phi0,double angle)

   set the referential of the map (angles in degrees)

    (default x0=siz_x/2,  y0=siz_y/2) 
*/
template<class T>
void LocalMap<T>::SetOrigin(double theta0,double phi0,double angle)
{
  theta0_= theta0;
  phi0_  = phi0;
  angle_ = angle;
  x0_= nSzX_/2;
  y0_= nSzY_/2;
  cos_angle_= cos(angle*Pi/180.);
  sin_angle_= sin(angle*Pi/180.);
  originFlag_= true; 
  cout << " LocalMap:: set origin 1 done" << endl;
}
 
/*!  \fn void SOPHYA::LocalMap::SetOrigin(double theta0,double phi0,int_4 x0,int_4 y0,double angle)

  set the referential of the map (angles in degrees) 
*/
template<class T>
void LocalMap<T>::SetOrigin(double theta0,double phi0,int_4 x0,int_4 y0,double angle)
{
  theta0_= theta0;
  phi0_  = phi0;
  angle_ = angle;
  x0_= x0;
  y0_= y0;
  cos_angle_= cos(angle*Pi/180.);
  sin_angle_= sin(angle*Pi/180.);
  originFlag_= true;
  cout << " LocalMap:: set origin 2 done" << endl;
} 

/*! \fn void SOPHYA::LocalMap::SetSize(double angleX,double angleY)

   angle range of tthe map (angles in degrees) 
*/
template<class T>
void LocalMap<T>::SetSize(double angleX,double angleY)
{
  angleX_= angleX;
  angleY_= angleY;

  // tangente de la moitie de l'ouverture angulaire totale
  tgAngleX_= tan(0.5*angleX_*Pi/180.);
  tgAngleY_= tan(0.5*angleY_*Pi/180.);

  extensFlag_= true;
  cout << " LocalMap:: set extension done" << endl;
}

/*! \fn void SOPHYA::LocalMap::Project(SphericalMap<T>& sphere) const

Projection to a spherical map 
*/
template<class T>
void LocalMap<T>::Project(SphericalMap<T>& sphere) const
{
  for(int m = 0; m < nPix_; m++) 
    {
      double theta,phi;
      PixThetaPhi(m,theta,phi);
      sphere(theta,phi)= pixels_(m);
      // cout << "theta " << theta << " phi " << phi << " valeur " << sphere(theta,phi)<< endl;
    }
}
// Titre        Private Methods
//++
template<class T>
void LocalMap<T>::InitNul()
//
//    set some attributes to zero
//--
{ 
  originFlag_= false;
  extensFlag_= false;
  cos_angle_= 1.0;
  sin_angle_= 0.0;
//  pixels_.Reset(); Pas de reset par InitNul (en cas de share) - Reza 20/11/99 $CHECK$
}


/*!   \fn void SOPHYA::LocalMap::Getij(int_4 k,int_4& i,int_4& j) const

 \return 2 indices corresponding to the pixel number k 
*/
template<class T>
void LocalMap<T>::Getij(int_4 k,int_4& i,int_4& j) const
{
  i= (k+1)%nSzX_-1;
  if(i == -1) i= nSzX_-1;
  j= (k-i+2)/nSzX_;
}

/*! \fn void  SOPHYA::LocalMap::ReferenceToUser(double& theta,double& phi) const

   Transform a pair of coordinates (theta, phi) given in 
    reference coordinates into map coordinates
*/
template<class T>
void  LocalMap<T>::ReferenceToUser(double& theta,double& phi) const
{
  if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi) 
    {
      throw "LocalMap::ReferenceToUser (theta,phi) out of range"; 
    } 

  theta= theta0_*Pi/180.+theta-Pi*0.5;
  if(theta < 0.) 
    {
      theta= -theta;
      phi += Pi;
    } 
  else 
    {
      if(theta > Pi) 
        {
          theta= 2.*Pi-theta;
          phi += Pi;
        } 
    }

  phi= phi0_*Pi/180.+phi;
  while(phi >= 2.*Pi) phi-= 2.*Pi;

  if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi) 
    {
      cout <<  " LocalMap::ReferenceToUser : erreur bizarre dans le transfert a la carte utilisateur " << endl;
      cout << " theta= " << theta << " phi= " << phi << endl;
      exit(0);
    }
}

/*!  \fn void  SOPHYA::LocalMap::UserToReference(double& theta,double& phi) const

  Transform a pair of coordinates (theta, phi) given in 
   map coordinates into reference coordinates
*/
template<class T>
void  LocalMap<T>::UserToReference(double& theta,double& phi) const
{
  if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi) 
    {
      cout<<" LocalMap::UserToReference: exceptions a mettre en place" <<endl;
      // THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
      throw "LocalMap::UserToReference (theta,phi) out of range"; 
    }

  double phi1= phi-phi0_*Pi/180.;
  if(phi1 < 0.) phi1+= 2.*Pi;

  double theta1= theta-theta0_*Pi/180.+Pi*0.5;
  if(theta1 < 0.) 
    {
      theta= -theta1;
      phi1+= Pi;
    } 
  else
    { 
      if(theta1 > Pi) 
        {
          theta= 2.*Pi-theta1;
          phi1+= Pi;
        }
    }

  while(phi1 >= 2.*Pi) phi1-= 2.*Pi;
  phi= phi1;
  if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi) 
    {
      cout <<  " LocalMap::UserToReference : erreur bizarre dans le transfert a la carte de reference " << endl;
      cout << " theta= " << theta << " phi= " << phi << endl;
      exit(0);
    }
}

/*! \fn void SOPHYA::LocalMap::PixProjToAngle(double x,double y,double& theta,double& phi) const

  Given coordinates in pixel units in the REFERENCE PLANE, return
    (theta, phi) in "absolute" referential theta=pi/2 ,phi=0.   
*/
template<class T>
void LocalMap<T>::PixProjToAngle(double x,double y,double& theta,double& phi) const
{
  theta= Pi*0.5-atan(2.*y*tgAngleY_/(double)nSzY_);
  phi= atan2(2.*x*tgAngleX_,(double)nSzX_);
  if(phi < 0.) phi += DeuxPi;
}

/*! \fn void SOPHYA::LocalMap::AngleProjToPix(double theta,double phi,double& x,double& y) const

   Given coordinates  (theta, phi) in "absolute" referential 
    theta=pi/2 ,phi=0  return pixel indices  (i,j) in the REFERENCE PLANE.
*/
template<class T>
void LocalMap<T>::AngleProjToPix(double theta,double phi,double& x,double& y) const
{
  if(phi >= Pi) phi-= DeuxPi;
  //  y=0.5*mSzY_*cot(theta)/tgAngleY_;  $CHECK-REZA-04/99$
  y= 0.5*nSzY_/tan(theta)/tgAngleY_;  // ? cot = 1/tan ?  
  x= 0.5*nSzX_*tan(phi)/tgAngleX_;
}

template<class T>
void LocalMap<T>::print(ostream& os) const
{
  os<<" SzX= "<<nSzX_<<", SzY= "<<nSzY_<<", NPix= "<<nPix_<<endl;
  if(LocalMap_isDone())
    {
      os<<" theta0= "<<theta0_<<", phi0= "<<phi0_<<", angle= "<<angle_<<endl;
      os<<" x0= "<<x0_<<", y0= "<<y0_<<endl;
      os<<" cos= "<<cos_angle_<<", & sin= "<<sin_angle_<<endl;
      os<<" angleX= "<<angleX_<<", angleY= "<<angleY_<<endl;
      os<<" tg(angleX)= "<<tgAngleX_<<", tg(angleY)= "<<tgAngleY_<<endl;
    }

  os << " contenu de pixels : ";
  for(int i=0; i < nPix_; i++)
    {
      if(i%5 == 0) os << endl;
      os <<  pixels_(i) <<", ";
    }
  os << endl;
}


template<class T>
void LocalMap<T>::recopierVariablesSimples(const LocalMap<T>& lm)
{
  nSzX_= lm.nSzX_;
  nSzY_= lm.nSzY_;
  nPix_= lm.nPix_;
  originFlag_= lm.originFlag_;
  extensFlag_= lm.extensFlag_;
  x0_= lm.x0_;
  y0_= lm.y0_;
  theta0_= lm.theta0_;
  phi0_= lm.phi0_;
  angle_= lm.angle_;
  cos_angle_= lm.cos_angle_;
  sin_angle_= lm.sin_angle_;
  angleX_= lm.angleX_;
  angleY_= lm.angleY_;
  tgAngleX_= lm.tgAngleX_;
  tgAngleY_= lm.tgAngleY_;
}

//   ...... Operations de calcul  ......


//! Fill a LocalMap with a constant value \b a
template <class T>
LocalMap<T>& LocalMap<T>::SetT(T a) 
{
  if (NbPixels() < 1) 
    throw RangeCheckError("LocalMap<T>::SetT(T )  - LocalMap not dimensionned ! ");
  pixels_ = a;
  return (*this);
}

/*! Add a constant value \b x to a LocalMap */
template <class T>
LocalMap<T>& LocalMap<T>::Add(T a)
 {
  if (NbPixels()< 1) 
    throw RangeCheckError("LocalMap<T>::Add(T )  - LocalMap not dimensionned ! ");
  pixels_ += a; 
  return (*this);
}

/*! Substract a constant value \b a to a LocalMap */
template <class T>
LocalMap<T>& LocalMap<T>::Sub(T a,bool fginv) 
{
  if (NbPixels()< 1) 
    throw RangeCheckError("LocalMap<T>::Sub(T )  - LocalMap not dimensionned ! ");
  pixels_.Sub(a,fginv); 
  return (*this);
} 

/*! multiply a LocalMap by a constant value \b a */
template <class T>
LocalMap<T>& LocalMap<T>::Mul(T a) 
{
  if (NbPixels()< 1) 
    throw RangeCheckError("LocalMap<T>::Mul(T )  - LocalMap not dimensionned ! ");
  pixels_ *= a; 
  return (*this);
}

/*! divide a LocalMap by a constant value \b a */
template <class T>
LocalMap<T>& LocalMap<T>::Div(T a) 
{
  if (NbPixels()< 1) 
    throw RangeCheckError("LocalMap<T>::Div(T )  - LocalMap not dimensionned ! ");
  pixels_ /= a; 
  return (*this);
} 

//  >>>> Operations avec 2nd membre de type LocalMap
//! Add two LocalMap

template <class T>
LocalMap<T>& LocalMap<T>::AddElt(const LocalMap<T>& a)
{
  if (nSzX_ != a.nSzX_ || nSzY_ != a.nSzY_)
    {
    throw(SzMismatchError("LocalMap<T>::AddElt(const LocalMap<T>&) SizeMismatch")) ;
    }
  pixels_ += a.pixels_;
  return (*this);
}

//! Substract two LocalMap
template <class T>
LocalMap<T>& LocalMap<T>::SubElt(const LocalMap<T>& a)
{
  if (nSzX_ != a.nSzX_ || nSzY_ != a.nSzY_)
    {
    throw(SzMismatchError("LocalMap<T>::SubElt(const LocalMap<T>&) SizeMismatch")) ;
    }
  pixels_ -= a.pixels_;
  return (*this);
}

//! Multiply two LocalMap (elements by elements)
template <class T>
LocalMap<T>& LocalMap<T>::MulElt(const LocalMap<T>& a)
{
  if (nSzX_ != a.nSzX_ || nSzY_ != a.nSzY_)
    {
    throw(SzMismatchError("LocalMap<T>::MulElt(const LocalMap<T>&) SizeMismatch")) ;
    }
  pixels_ *= a.pixels_;
  return (*this);
}

//! Divide two LocalMaps (elements by elements) - No protection for divide by 0
template <class T>
LocalMap<T>& LocalMap<T>::DivElt(const LocalMap<T>& a)
{
  if (nSzX_ != a.nSzX_ || nSzY_ != a.nSzY_)
    {
    throw(SzMismatchError("LocalMap<T>::DivElt(const LocalMap<T>&) SizeMismatch")) ;
    }
  pixels_ /= a.pixels_;
  return (*this);
}





//++
// Titre        class FIO_LocalMap
//    Delegated objects for persitance management
//--


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