Index: trunk/SophyaLib/Samba/localmap.cc
===================================================================
--- trunk/SophyaLib/Samba/localmap.cc	(revision 604)
+++ trunk/SophyaLib/Samba/localmap.cc	(revision 682)
@@ -6,4 +6,729 @@
 #include <string.h>
 #include <iostream.h>
+
+
+//*****************************************************************************
+//++ 
+// Class	LocalMap
+//
+// include      localmap.h nbmath.h
+//
+//    A local map of a region of the sky, in cartesian coordinates.
+//    It 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(...))
+//    
+//
+
+//
+//-- 
+//++
+//
+// Links	Parents
+//
+//    PixelMap
+//    
+//-- 
+//++
+//
+// Links	Childs
+//
+//     
+//-- 
+//++
+// Titre	Constructors
+//--
+//++
+template<class T>
+LocalMap<T>::LocalMap()
+//
+//     
+//--
+{
+  InitNul();
+  pixels_.Reset();
+}
+
+//++
+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)
+
+//
+//    copy constructor
+//--
+{
+  cout<<" LocalMap:: Appel du constructeur de recopie " << endl;
+
+  if(lm.mInfo_) mInfo_= new DVList(*lm.mInfo_);
+  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_;
+}
+
+//++
+// Titre	Destructor
+//--
+//++
+template<class T>
+LocalMap<T>::~LocalMap()
+//
+//--
+{
+  InitNul();
+}
+
+
+
+//++
+// Titre	Public Methods
+//--
+
+//++
+template<class T>
+void LocalMap<T>::ReSize(int_4 nx, int_4 ny) 
+//
+//    Resize storage area for pixels 
+//--
+{
+  InitNul();
+  nSzX_ = nx;
+  nSzY_ = ny;
+  nPix_= nx*ny;
+  pixels_.ReSize(nPix_);
+  pixels_.Reset();
+}
+
+//++
+template<class T>
+int_4 LocalMap<T>::NbPixels() const
+//
+//    Return number of pixels
+//--
+{
+  return(nPix_);
+}
+
+//++
+template<class T>
+T& LocalMap<T>::PixVal(int_4 k)
+//
+//    Return value of pixel with index k
+//--
+{
+  if((k < 0) || (k >= nPix_)) 
+    {
+      cout << " LocalMap::PIxVal : exceptions  a mettre en place" <<endl;
+      //  THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
+      THROW(rangeCheckErr);
+      //throw "LocalMap::PIxVal Pixel index out of range";
+    }
+  return(pixels_(k));
+}
+
+//++
+
+template<class T>
+T const& LocalMap<T>::PixVal(int_4 k) const
+//
+//   const version of previous method
+//--
+{
+  if((k < 0) || (k >= nPix_)) 
+    {
+      cout << " LocalMap::PIxVal : exceptions  a mettre en place" <<endl;
+      //    THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
+      
+      throw "LocalMap::PIxVal Pixel index out of range"; 
+    } 
+  return *(pixels_.Data()+k);
+}
+
+template<class T>
+bool LocalMap<T>::ContainsSph(double theta, double phi) const
+{
+// $CHECK$  Reza 11/11/99 - A modifier
+return(true);
+}
+
+//++
+template<class T>
+int_4 LocalMap<T>::PixIndexSph(double theta,double phi) const
+//
+//    Return index of the pixel with spherical coordinates (theta,phi)
+//--
+{
+  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);
+}
+
+//++
+template<class T>
+void LocalMap<T>::PixThetaPhi(int_4 k,double& theta,double& phi) const
+//
+//    Return (theta, phi) coordinates of pixel with index k
+//--
+{
+  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);
+}
+
+template <class T>
+T LocalMap<T>::SetPixels(T v)
+{
+pixels_.Reset(v);
+return(v);
+}
+ 
+//++
+template<class T>
+double LocalMap<T>::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.   
+//--
+{
+  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;
+}
+
+//++
+template<class T>
+void LocalMap<T>::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) 
+//--
+{
+  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;
+}
+ 
+//++
+template<class T>
+void LocalMap<T>::SetOrigin(double theta0,double phi0,int_4 x0,int_4 y0,double angle)
+//
+//    set the referential of the map (angles in degrees)
+//--
+{
+  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;
+} 
+
+//++
+template<class T>
+void LocalMap<T>::SetSize(double angleX,double angleY)
+//
+//    angle range of tthe map (angles in degrees)
+//--
+{
+  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;
+}
+
+//++
+template<class T>
+void LocalMap<T>::Project(SphericalMap<T>& sphere) const
+//
+//    Projection to a spherical map   
+//--
+{
+  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$
+}
+
+//++
+template<class T>
+void LocalMap<T>::Getij(int_4 k,int_4& i,int_4& j) const
+//
+//    Return 2 indices corresponding to the pixel number k
+//--
+{
+  i= (k+1)%nSzX_-1;
+  if(i == -1) i= nSzX_-1;
+  j= (k-i+2)/nSzX_;
+}
+
+//++
+template<class T>
+void  LocalMap<T>::ReferenceToUser(double& theta,double& phi) const
+//
+//    Transform a pair of coordinates (theta, phi) given in 
+//    reference coordinates into map coordinates
+//--     
+{
+  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);
+    }
+}
+
+//++
+template<class T>
+void  LocalMap<T>::UserToReference(double& theta,double& phi) const
+//
+//    Transform a pair of coordinates (theta, phi) given in 
+//    map coordinates into reference coordinates
+//--     
+{
+  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);
+    }
+}
+
+//++
+template<class T>
+void LocalMap<T>::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.   
+//--
+{
+  theta= Pi*0.5-atan(2.*y*tgAngleY_/(double)nSzY_);
+  phi= atan2(2.*x*tgAngleX_,(double)nSzX_);
+  if(phi < 0.) phi += DeuxPi;
+}
+
+//++
+template<class T>
+void LocalMap<T>::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.
+//--
+{
+  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;
+}
+//++
+// Titre	class FIO_LocalMap
+//    Delegated objects for persitance management
+//--
+
+//*******************************************************************
+// class FIO_LocalMap<T>
+//  Les objets delegues pour la gestion de persistance 
+//*******************************************************************
+
+//++
+template <class T>
+FIO_LocalMap<T>::FIO_LocalMap()
+//
+//--
+{
+  dobj= new LocalMap<T>;
+  ownobj= true;
+}
+//++
+template <class T>
+FIO_LocalMap<T>::FIO_LocalMap(string const& filename)
+//
+//--
+{
+  dobj= new LocalMap<T>;
+  dobj->DataBlock().SetTemp(true);
+  ownobj= true;
+  Read(filename);
+}
+
+//++
+template <class T>
+FIO_LocalMap<T>::FIO_LocalMap(const LocalMap<T>& obj)
+//
+//--
+{
+  dobj= new LocalMap<T>(obj, true);
+  dobj->DataBlock().SetTemp(true);
+  ownobj= true;
+}
+
+template <class T>
+FIO_LocalMap<T>::FIO_LocalMap(LocalMap<T>* obj)
+{
+  dobj= obj;
+  ownobj= false;
+}
+
+//++
+template <class T>
+FIO_LocalMap<T>::~FIO_LocalMap()
+//
+//--
+{
+  if (ownobj && dobj) delete dobj;
+}
+//++
+template <class T>
+AnyDataObj* FIO_LocalMap<T>::DataObj()
+//
+//--
+{
+  return(dobj);
+}
+
+//++
+template <class T>
+void FIO_LocalMap<T>::ReadSelf(PInPersist& is)
+//
+//--
+{
+
+  if(dobj == NULL) 
+    {
+      dobj= new LocalMap<T>;
+      dobj->DataBlock().SetTemp(true);
+      ownobj= true;
+    }
+
+  // Pour savoir s'il y avait un DVList Info associe
+  char strg[256];
+  is.GetLine(strg, 255);
+  bool hadinfo= false;
+  if(strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0)  hadinfo= true;
+  if(hadinfo) 
+    {    // Lecture eventuelle du DVList Info
+      is >> dobj->Info(); 
+    }
+
+  int_4 nSzX;
+  is.GetI4(nSzX);
+  dobj->setSize_x(nSzX);
+
+  int_4 nSzY;
+  is.GetI4(nSzY);
+  dobj->setSize_y(nSzY);
+
+  int_4 nPix;
+  is.GetI4(nPix);
+  dobj->setNbPixels(nPix);
+
+  string ss("local mapping is done");
+  string sso;
+  is.GetStr(sso);
+  if(sso == ss)
+    {
+      cout<<" ReadSelf:: local mapping"<<endl;
+      int_4 x0, y0;
+      double theta, phi, angle;
+      is.GetI4(x0);
+      is.GetI4(y0);
+      is.GetR8(theta);
+      is.GetR8(phi);
+      is.GetR8(angle);
+      dobj->SetOrigin(theta, phi, x0, y0, angle);
+
+      double angleX, angleY;
+      is.GetR8(angleX);
+      is.GetR8(angleY);
+      dobj->SetSize(angleX, angleY);
+    }
+
+// On lit le DataBlock;
+  is >> dobj->DataBlock();
+}
+
+//++
+template <class T>
+void FIO_LocalMap<T>::WriteSelf(POutPersist& os) const
+//
+//--
+{
+  if(dobj == NULL) 
+    {
+      cout << " FIO_LocalMap::WriteSelf:: dobj= null " << endl;
+      return;
+    }
+
+  char strg[256];
+  int_4 nSzX= dobj->Size_x();
+  int_4 nSzY= dobj->Size_y();
+  int_4 nPix= dobj->NbPixels();
+ 
+  if(dobj->ptrInfo()) 
+    {
+      sprintf(strg,"LocalMap: NPixX=%6d  NPixY=%9d HasInfo",nSzX,nSzY);
+      os.PutLine(strg);
+      os << dobj->Info();
+    }
+  else 
+    { 
+      sprintf(strg,"LocalMap: NPixX=%6d  NPixY=%9d ",nSzX,nSzY);
+      os.PutLine(strg);  
+    }
+
+  os.PutI4(nSzX);
+  os.PutI4(nSzY);
+  os.PutI4(nPix);
+
+  if(dobj->LocalMap_isDone())
+    {
+      string ss("local mapping is done");
+      os.PutStr(ss);
+      int_4 x0, y0;
+      double theta, phi, angle;
+      dobj->Origin(theta, phi, x0, y0, angle);
+      os.PutI4(x0);
+      os.PutI4(y0);
+      os.PutR8(theta);
+      os.PutR8(phi);
+      os.PutR8(angle);
+
+      double angleX, angleY;
+      dobj->Aperture(angleX, angleY);
+      os.PutR8(angleX);
+      os.PutR8(angleY);
+    }
+  else
+    {
+      string ss("no local mapping");
+      os.PutStr(ss);
+    }
+
+// On ecrit le dataBlock 
+  os << dobj->DataBlock();
+}
 
 #ifdef __CXX_PRAGMA_TEMPLATES__
@@ -27,727 +752,2 @@
 template class FIO_LocalMap< complex<double> >;
 #endif
-
-//*****************************************************************************
-//++ 
-// Class	LocalMap
-//
-// include      localmap.h nbmath.h
-//
-//    A local map of a region of the sky, in cartesian coordinates.
-//    It 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(...))
-//    
-//
-
-//
-//-- 
-//++
-//
-// Links	Parents
-//
-//    PixelMap
-//    
-//-- 
-//++
-//
-// Links	Childs
-//
-//     
-//-- 
-//++
-// Titre	Constructors
-//--
-//++
-template<class T>
-LocalMap<T>::LocalMap()
-//
-//     
-//--
-{
-  InitNul();
-  pixels_.Reset();
-}
-
-//++
-template<class T>
-LocalMap<T>::LocalMap(int nx, int 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)
-
-//
-//    copy constructor
-//--
-{
-  cout<<" LocalMap:: Appel du constructeur de recopie " << endl;
-
-  if(lm.mInfo_) mInfo_= new DVList(*lm.mInfo_);
-  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_;
-}
-
-//++
-// Titre	Destructor
-//--
-//++
-template<class T>
-LocalMap<T>::~LocalMap()
-//
-//--
-{
-  InitNul();
-}
-
-
-
-//++
-// Titre	Public Methods
-//--
-
-//++
-template<class T>
-void LocalMap<T>::ReSize(int nx, int ny) 
-//
-//    Resize storage area for pixels 
-//--
-{
-  InitNul();
-  nSzX_ = nx;
-  nSzY_ = ny;
-  nPix_= nx*ny;
-  pixels_.ReSize(nPix_);
-  pixels_.Reset();
-}
-
-//++
-template<class T>
-int LocalMap<T>::NbPixels() const
-//
-//    Return number of pixels
-//--
-{
-  return(nPix_);
-}
-
-//++
-template<class T>
-T& LocalMap<T>::PixVal(int k)
-//
-//    Return value of pixel with index k
-//--
-{
-  if((k < 0) || (k >= nPix_)) 
-    {
-      cout << " LocalMap::PIxVal : exceptions  a mettre en place" <<endl;
-      //  THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
-      THROW(rangeCheckErr);
-      //throw "LocalMap::PIxVal Pixel index out of range";
-    }
-  return(pixels_(k));
-}
-
-//++
-
-template<class T>
-T const& LocalMap<T>::PixVal(int k) const
-//
-//   const version of previous method
-//--
-{
-  if((k < 0) || (k >= nPix_)) 
-    {
-      cout << " LocalMap::PIxVal : exceptions  a mettre en place" <<endl;
-      //    THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
-      
-      throw "LocalMap::PIxVal Pixel index out of range"; 
-    } 
-  return *(pixels_.Data()+k);
-}
-
-template<class T>
-bool LocalMap<T>::ContainsSph(double theta, double phi) const
-{
-// $CHECK$  Reza 11/11/99 - A modifier
-return(true);
-}
-
-//++
-template<class T>
-int LocalMap<T>::PixIndexSph(double theta,double phi) const
-//
-//    Return index of the pixel with spherical coordinates (theta,phi)
-//--
-{
-  int 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);
-}
-
-//++
-template<class T>
-void LocalMap<T>::PixThetaPhi(int k,double& theta,double& phi) const
-//
-//    Return (theta, phi) coordinates of pixel with index k
-//--
-{
-  if(!(originFlag_) || !(extensFlag_)) 
-    {
-      cout << " LocalMap: correspondance carte-sphere non etablie" << endl;
-      exit(0);
-    }
-
-  int 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);
-}
-
-template <class T>
-T LocalMap<T>::SetPixels(T v)
-{
-pixels_.Reset(v);
-return(v);
-}
- 
-//++
-template<class T>
-double LocalMap<T>::PixSolAngle(int 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.   
-//--
-{
-  int 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;
-}
-
-//++
-template<class T>
-void LocalMap<T>::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) 
-//--
-{
-  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;
-}
- 
-//++
-template<class T>
-void LocalMap<T>::SetOrigin(double theta0,double phi0,int x0,int y0,double angle)
-//
-//    set the referential of the map (angles in degrees)
-//--
-{
-  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;
-} 
-
-//++
-template<class T>
-void LocalMap<T>::SetSize(double angleX,double angleY)
-//
-//    angle range of tthe map (angles in degrees)
-//--
-{
-  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;
-}
-
-//++
-template<class T>
-void LocalMap<T>::Project(SphericalMap<T>& sphere) const
-//
-//    Projection to a spherical map   
-//--
-{
-  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$
-}
-
-//++
-template<class T>
-void LocalMap<T>::Getij(int k,int& i,int& j) const
-//
-//    Return 2 indices corresponding to the pixel number k
-//--
-{
-  i= (k+1)%nSzX_-1;
-  if(i == -1) i= nSzX_-1;
-  j= (k-i+2)/nSzX_;
-}
-
-//++
-template<class T>
-void  LocalMap<T>::ReferenceToUser(double& theta,double& phi) const
-//
-//    Transform a pair of coordinates (theta, phi) given in 
-//    reference coordinates into map coordinates
-//--     
-{
-  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);
-    }
-}
-
-//++
-template<class T>
-void  LocalMap<T>::UserToReference(double& theta,double& phi) const
-//
-//    Transform a pair of coordinates (theta, phi) given in 
-//    map coordinates into reference coordinates
-//--     
-{
-  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);
-    }
-}
-
-//++
-template<class T>
-void LocalMap<T>::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.   
-//--
-{
-  theta= Pi*0.5-atan(2.*y*tgAngleY_/(double)nSzY_);
-  phi= atan2(2.*x*tgAngleX_,(double)nSzX_);
-  if(phi < 0.) phi += DeuxPi;
-}
-
-//++
-template<class T>
-void LocalMap<T>::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.
-//--
-{
-  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;
-}
-//++
-// Titre	class FIO_LocalMap
-//    Delegated objects for persitance management
-//--
-
-//*******************************************************************
-// class FIO_LocalMap<T>
-//  Les objets delegues pour la gestion de persistance 
-//*******************************************************************
-
-//++
-template <class T>
-FIO_LocalMap<T>::FIO_LocalMap()
-//
-//--
-{
-  dobj= new LocalMap<T>;
-  ownobj= true;
-}
-//++
-template <class T>
-FIO_LocalMap<T>::FIO_LocalMap(string const& filename)
-//
-//--
-{
-  dobj= new LocalMap<T>;
-  dobj->DataBlock().SetTemp(true);
-  ownobj= true;
-  Read(filename);
-}
-
-//++
-template <class T>
-FIO_LocalMap<T>::FIO_LocalMap(const LocalMap<T>& obj)
-//
-//--
-{
-  dobj= new LocalMap<T>(obj, true);
-  dobj->DataBlock().SetTemp(true);
-  ownobj= true;
-}
-
-template <class T>
-FIO_LocalMap<T>::FIO_LocalMap(LocalMap<T>* obj)
-{
-  dobj= obj;
-  ownobj= false;
-}
-
-//++
-template <class T>
-FIO_LocalMap<T>::~FIO_LocalMap()
-//
-//--
-{
-  if (ownobj && dobj) delete dobj;
-}
-//++
-template <class T>
-AnyDataObj* FIO_LocalMap<T>::DataObj()
-//
-//--
-{
-  return(dobj);
-}
-
-//++
-template <class T>
-void FIO_LocalMap<T>::ReadSelf(PInPersist& is)
-//
-//--
-{
-
-  if(dobj == NULL) 
-    {
-      dobj= new LocalMap<T>;
-      dobj->DataBlock().SetTemp(true);
-      ownobj= true;
-    }
-
-  // Pour savoir s'il y avait un DVList Info associe
-  char strg[256];
-  is.GetLine(strg, 255);
-  bool hadinfo= false;
-  if(strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0)  hadinfo= true;
-  if(hadinfo) 
-    {    // Lecture eventuelle du DVList Info
-      is >> dobj->Info(); 
-    }
-
-  int nSzX;
-  is.GetI4(nSzX);
-  dobj->setSize_x(nSzX);
-
-  int nSzY;
-  is.GetI4(nSzY);
-  dobj->setSize_y(nSzY);
-
-  int nPix;
-  is.GetI4(nPix);
-  dobj->setNbPixels(nPix);
-
-  string ss("local mapping is done");
-  string sso;
-  is.GetStr(sso);
-  if(sso == ss)
-    {
-      cout<<" ReadSelf:: local mapping"<<endl;
-      int x0, y0;
-      double theta, phi, angle;
-      is.GetI4(x0);
-      is.GetI4(y0);
-      is.GetR8(theta);
-      is.GetR8(phi);
-      is.GetR8(angle);
-      dobj->SetOrigin(theta, phi, x0, y0, angle);
-
-      double angleX, angleY;
-      is.GetR8(angleX);
-      is.GetR8(angleY);
-      dobj->SetSize(angleX, angleY);
-    }
-
-// On lit le DataBlock;
-  is >> dobj->DataBlock();
-}
-
-//++
-template <class T>
-void FIO_LocalMap<T>::WriteSelf(POutPersist& os) const
-//
-//--
-{
-  if(dobj == NULL) 
-    {
-      cout << " FIO_LocalMap::WriteSelf:: dobj= null " << endl;
-      return;
-    }
-
-  char strg[256];
-  int nSzX= dobj->Size_x();
-  int nSzY= dobj->Size_y();
-  int nPix= dobj->NbPixels();
- 
-  if(dobj->ptrInfo()) 
-    {
-      sprintf(strg,"LocalMap: NPixX=%6d  NPixY=%9d HasInfo",nSzX,nSzY);
-      os.PutLine(strg);
-      os << dobj->Info();
-    }
-  else 
-    { 
-      sprintf(strg,"LocalMap: NPixX=%6d  NPixY=%9d ",nSzX,nSzY);
-      os.PutLine(strg);  
-    }
-
-  os.PutI4(nSzX);
-  os.PutI4(nSzY);
-  os.PutI4(nPix);
-
-  if(dobj->LocalMap_isDone())
-    {
-      string ss("local mapping is done");
-      os.PutStr(ss);
-      int x0, y0;
-      double theta, phi, angle;
-      dobj->Origin(theta, phi, x0, y0, angle);
-      os.PutI4(x0);
-      os.PutI4(y0);
-      os.PutR8(theta);
-      os.PutR8(phi);
-      os.PutR8(angle);
-
-      double angleX, angleY;
-      dobj->Aperture(angleX, angleY);
-      os.PutR8(angleX);
-      os.PutR8(angleY);
-    }
-  else
-    {
-      string ss("no local mapping");
-      os.PutStr(ss);
-    }
-
-// On ecrit le dataBlock 
-  os << dobj->DataBlock();
-}
-
