Changeset 2013 in Sophya for trunk

Timestamp:

May 24, 2002, 4:54:53 PM (23 years ago)

Author:

lemeur

Message:

amelioration de la projection de localmap

Location:

trunk

Files:

• SophyaExt/FitsIOServer/fitslocalmap.cc (modified) (5 diffs)
• SophyaLib/SkyMap/fiolocalmap.cc (modified) (5 diffs)
• SophyaLib/SkyMap/localmap.cc (modified) (22 diffs)
• SophyaLib/SkyMap/localmap.h (modified) (6 diffs)

trunk/SophyaExt/FitsIOServer/fitslocalmap.cc

@@ -5 +5 @@
 //    pout LocalMap
 ///////////////////////////////////////////////////////////
-
+ 
 /*!
   \class SOPHYA::FITS_LocalMap
@@ -93 +93 @@
     }
   dobj_->ReSize(nSzX, nSzY); 
-  int_4 localMappingDone = dvl.GetI("LCMP");
-  if (localMappingDone)
-    {
+  //  int_4 localMappingDone = dvl.GetI("LCMP");
+  // if (localMappingDone)
+  //  {
       int_4 x0 = dvl.GetI("X0");
       int_4 y0 = dvl.GetI("Y0");
@@ -105 +105 @@
       double angleY = dvl.GetD("ANGLEY");
       dobj_->SetSize(angleX, angleY);
-    }
+      //  }
   // On lit les DataBlocks;
   dobj_->DataBlock().ReSize(nPix);
@@ -131 +131 @@
   dvl.SetComment("NPIX", "number of pixels in local map");
 
-  if(dobj_->LocalMap_isDone())
-    {
-      dvl["LCMP"] = (int_4) 1;
-      dvl.SetComment("LCMP", "local mapping 1= done, 0= not done");
+  // if(dobj_->LocalMap_isDone())
+  // {
+  //   dvl["LCMP"] = (int_4) 1;
+  //   dvl.SetComment("LCMP", "local mapping 1= done, 0= not done");
       int_4 x0 = 0;
       int_4 y0 = 0;
@@ -159 +159 @@
       dvl["ANGLEY"] = angleY;
       dvl.SetComment("ANGLEY", "sphere angle covered by map y-extension ");
-    }
-  else
-    {
-      dvl["LCMP"] = (int_4) 0;
-    }
+      // }
+  // else
+  //  {
+  //    dvl["LCMP"] = (int_4) 0;
+  //  }
   dvl["Content"]= "LocalMap";
   dvl.SetComment("Content", "name of SOPHYA object");

trunk/SophyaLib/SkyMap/fiolocalmap.cc

@@ -2 +2 @@
 
 #include "fiolocalmap.h"
+#include "fioarr.h"
 #include "pexceptions.h"
 #include "fiondblock.h"
@@ -115 +116 @@
   //  dobj->setSize_y(nSzY);
 
-  int_4 nPix;
-  is.GetI4(nPix);
+  // int_4 nPix;
+  // is.GetI4(nPix);
   //  dobj->setNbPixels(nPix);
   dobj->ReSize(nSzX, nSzY); 
-  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);
+  //  string ss("local mapping is done");
+  // string sso;
+  // is.GetStr(sso);
+  // if(sso == ss)
+  //  {
+  //   cout<<" ReadSelf:: local mapping"<<endl;
 
       double angleX, angleY;
@@ -138 +131 @@
       is.GetR8(angleY);
       dobj->SetSize(angleX, angleY);
-    }
+
+
+      int_4 x0, y0;
+      double theta, phi, angle;
+      is.GetR8(theta);
+      is.GetR8(phi);
+      is.GetI4(x0);
+      is.GetI4(y0);
+      is.GetR8(angle);
+      dobj->SetOrigin(theta, phi, x0, y0, angle);
+
+      //  }
 
 // On lit le DataBlock;
-  FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
+      //  FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
+// On lit la matrice;
+  FIO_TArray<T> fio_nd(&dobj->Matrix());
   fio_nd.Read(is);
 }
@@ -176 +182 @@
   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);
+  //  os.PutI4(nPix);
+
+  // if(dobj->LocalMap_isDone())
+  // {
+  //   string ss("local mapping is done");
+  //    os.PutStr(ss);
 
       double angleX, angleY;
@@ -195 +193 @@
       os.PutR8(angleX);
       os.PutR8(angleY);
-    }
-  else
-    {
-      string ss("no local mapping");
-      os.PutStr(ss);
-    }
+
+
+      int_4 x0, y0;
+      double theta, phi, angle;
+      dobj->Origin(theta, phi, x0, y0, angle);
+      os.PutR8(theta);
+      os.PutR8(phi);
+      os.PutI4(x0);
+      os.PutI4(y0);
+      os.PutR8(angle);
+
+      // }
+// else
+//  {
+//   string ss("no local mapping");
+//   os.PutStr(ss);
+//  }
 
 // On ecrit le dataBlock 
-  FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
+      //  FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
+// On ecrit la matrice 
+  FIO_TArray<T> fio_nd(&dobj->Matrix());
   fio_nd.Write(os);
 }

trunk/SophyaLib/SkyMap/localmap.cc

@@ -7 +7 @@
 #include <string.h>
 #include <iostream.h>
+#include "timing.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)
+  A local map is a 2 dimensional matrix, with ny rows and nx columns.
+  The local map has an origin in (theta0, phi0), mapped to pixel(x0, y0)
  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.
+ Each pixel(ip, it) is defined by its "phi-like" index ip (column 
+ index in the matrix) and its "theta-like" index it (row index in 
+ the matrix). Index ip is associated with x-axis in a map-coordinate 
+ system ; index it is associated with y-axis in the same map-coordinate system.
+  
+    The map is supposed to lie on a plan tangent to the celestial sphere 
+    at a point, with spherical coordinates (theta0, phi0), whose pixel 
+    numbers  are (x0,y0) in the local map. The aperture of the map is defined
+    by two values of angles, angleX and angleY, covered respectively by all 
+    the pixels in x direction and all the pixels in y direction.
+
+    Each pixel has angleX/nx and angleY/ny, as angle extensions. So, in
+    map-coordinate system the pixel (i,j) has following coordinates : 
+
+    x = (i-x0)*angleX/nx
+ 
+    y = (j-y0)*angleY/ny 
     
-    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(...))
+    (angles in radians)
+
+    The projection (method : ProjectionToSphere() )of the map onto a 
+    sphere is made by the following procedure : 
+
+    the sphere is supposed to have radius=1. The map is 
+    considered to be tangent to the sphere, on a point with (theta0, phi0)
+    spherical coodinates. A reference coordinate system (plane-coordinate
+    system) , is chosen in the plane of the map  with reference axes : 
+
+    x-axis : vector tengent to a parallel in (theta0, phi0) on the sphere
+    (components in "3-dim cartesian system : -sin(phi0) ; cos(phi0) ; 0)
+
+    z-axis : vector-radius with 3-dim cartesian cordinates : 
+             sin(theta0)*cos(phi0) ; sin(theta0*sin(phi0) ; cos(theta0)
+
+    y-axis = z-axis^x-axis : tangent to the meridian at (theta0, phi0)
+
+    note that the map-coordinate system may be rotated with respect to
+    plane-coordinate system (parameter "angle" below).
+    
+    the projection of a map pixel is defined as the intersection of 
+    the vector-radius, from sphere center to the pixel defined by
+    his coordinates  in the plane-coordinate system (computed from x,y 
+    above, with eventual rotation), with the sphere.
 */
 template<class T>
@@ -43 +65 @@
 }
 
-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) 
+//! full constructor of the local map
+/*!
+  \param nx : number of pixels in x direction
+  \param ny : number of pixels in y direction
+  \param angleX : total angle aperture in x direction (degrees)
+  \param angleY : total angle aperture in y direction (degrees)
+  \param theta0,phi0 : spherical coordinates of reference point at which 
+                       the map is considered to be tangent to the sphere 
+  \param x0, y0 : coodinates (in pixels) of the reference point just defined
+  \param angle : angle (degrees) of the rotation between x-axis of 
+                 map-coordinate system) and the tangent to parallel on 
+                 the sphere (default : 0.).
+ */
+template<class T>
+LocalMap<T>::LocalMap(int_4 nx, int_4 ny, double angleX,double angleY, double theta0,double phi0,int_4 x0,int_4 y0,double angle)  
 {
   InitNul();
@@ -90 +85 @@
   nSzY_ = ny;
   nPix_= nx*ny;
-  pixels_.ReSize(nPix_);
-  pixels_.Reset();
+  pixels_.ReSize(ny,nx);
+  SetSize(angleX, angleY);
+  SetOrigin(theta0, phi0, x0, y0, angle);
+}
+
+
+//! standard constructor of the local map
+/*!
+  \param nx : number of pixels in x direction
+  \param ny : number of pixels in y direction
+  \param angleX : total angle aperture in x direction (degrees)
+  \param angleY : total angle aperture in y direction (degrees)
+  \param theta0,phi0 : spherical coordinates of reference point at which 
+                       the map is considered to be tangent to the sphere 
+  \param angle : angle (degrees) of the rotation between x-axis of 
+                 map-coordinate system) and the tangent to parallel on 
+                 the sphere (default : 0.).
+ */
+template<class T>
+LocalMap<T>::LocalMap(int_4 nx, int_4 ny, double angleX,double angleY, double theta0,double phi0, double angle) 
+{
+  InitNul();
+  nSzX_ = nx;
+  nSzY_ = ny;
+  nPix_= nx*ny;
+  pixels_.ReSize(ny,nx);
+  cout << " tableau pixels initialise " << endl;
+  SetSize(angleX, angleY);
+  cout << " fin du setsize " << endl;
+  SetOrigin(theta0, phi0, angle);
+  cout << " fin du set oorigine " << endl;
+}
+
+//! copy constructor 
+/*!
+  \param share : if true, share data. If false copy data
+ */
+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);
+}
+
+//! copy constructor 
+/*!
+  \warning datas are \b SHARED with \b lm.
+  \sa TMatrix::TMatrix(const TMatrix<T>&)
+*/
+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() {;}
+
+
+
+/*!  \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) 
+{
+    // angles par pixel,  en radians 
+  deltaPhi_   *=  nSzX_/nx;   
+  deltaTheta_ *=  nSzY_/ny;
+
+  nSzX_ = nx;
+  nSzY_ = ny;
+  nPix_= nx*ny;
+  pixels_.ReSize(ny,nx);
 }
 
@@ -122 +196 @@
       recopierVariablesSimples(a);     
   int k;
-  for (k=0; k< nPix_; k++) pixels_(k) = a.pixels_(k);
+  pixels_ = a.pixels_;
   return(*this);
 
@@ -163 +237 @@
   if((k < 0) || (k >= nPix_)) 
     {
-      cout << " LocalMap::PIxVal : exceptions  a mettre en place" <<endl;
-      THROW(rangeCheckErr);
-    }
-  return(pixels_(k));
+    throw RangeCheckError("LocalMap::PIxVal Pixel index out of range ");
+    }
+  //
+  int_4 i,j;
+  Getij(k,i,j);
+  return(pixels_(j,i));
+
+  //
 }
 
@@ -177 +255 @@
   if((k < 0) || (k >= nPix_)) 
     {
-      cout << " LocalMap::PIxVal : exceptions  a mettre en place" <<endl;
-      throw "LocalMap::PIxVal Pixel index out of range"; 
+    throw RangeCheckError("LocalMap::PIxVal Pixel index out of range ");
     } 
-  return *(pixels_.Data()+k);
+  //
+  int_4 i,j;
+  Getij(k,i,j);
+  return(pixels_(j,i));
 }
 
@@ -197 +277 @@
 /*!  \fn int_4 SOPHYA::LocalMap::PixIndexSph(double theta,double phi) const
 
+   angles in radians
   \return index of the pixel with spherical coordinates (theta,phi) 
 */
@@ -203 +284 @@
 {
   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 csTheta = cos(theta);
+  double snTheta = sin(theta);
+  double csPhi   = cos(phi);
+  double snPhi   = sin(phi); 
+  double csPhiMPhiC = cos (phi - phiC_);
+  double snPhiMPhiC = sin (phi - phiC_);
+  // le point sur la sphere est note M. 
+  // intersection P de OM avec le plan de la carte (tangent en C)
+  double denom = snTheta*snthC_*csPhiMPhiC + csTheta*csthC_;
+  if ( denom == 0.) 
+    {
+      throw PException("LocalMap::PixIndexSph : vector radius parallel to map plane ");
+    }
+  double lambda = 1./denom;
+  // coordonnes du point d'intersection P dans le repere lie au plan
+
+  double XP = lambda*snTheta*snPhiMPhiC;
+  double YP = lambda*snTheta*csthC_*csPhiMPhiC;
+
+  // coordonnees dans le repere lie a la carte (eventuellement tourne par
+  // rapport au precedent.
+
+  double X = XP*cosAngle_ + YP*sinAngle_;
+  double Y = -XP*sinAngle_ + YP*cosAngle_;
+
+  // en unites de pixels 
+
+  X /= deltaPhi_;
+  Y /= deltaTheta_;
+
 
   double xmin= -x0_-0.5;
   double xmax= xmin+nSzX_;
-  if((x > xmax) || (x < xmin)) return(-1);
+  if((X > xmax) || (X < xmin)) return(-1);
   double xcurrent= xmin;
   for(i = 0; i < nSzX_; i++ ) 
     {
       xcurrent += 1.;
-      if( x < xcurrent ) break;
+      if( X < xcurrent ) break;
     }
   double ymin= -y0_-0.5;
   double ymax= ymin+nSzY_;
-  if((y >  ymax) || (y < ymin)) return(-1);
+  if((Y >  ymax) || (Y < ymin)) return(-1);
   double ycurrent= ymin;
   for(j = 0; j < nSzY_; j++ ) 
     {
       ycurrent += 1.;
-      if( y < ycurrent ) break;
+      if( Y < ycurrent ) break;
     }
   return (j*nSzX_+i);
@@ -243 +341 @@
    \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);
-}
+  PixThetaPhi(i,j, theta, phi);
+
+}
+
+/*! \fn void SOPHYA::LocalMap::PixThetaPhi(int_4 ip,int_4 it, double& theta,double& phi) const
+
+
+   \return (theta, phi) coordinates of pixel of map with indices (ip,it) corresponding to x and y directions
+
+ \param ip : phi-like index
+ \param it : theta-like index
+*/
+
+template<class T>
+void LocalMap<T>::PixThetaPhi(int_4 ip,int_4 it, double& theta,double& phi) const
+{
+  double XP, YP, ZP;
+  PixToSphereC(ip,it,XP,YP,ZP);
+
+  theta = acos(ZP);
+  phi = atan2(YP, XP);
+  while (phi < 0.)  phi += 2.*Pi;
+
+}
+
+
+
 
 /*! \fn T SOPHYA::LocalMap::SetPixels(T v)
@@ -273 +385 @@
 T LocalMap<T>::SetPixels(T v)
 {
-pixels_.Reset(v);
+pixels_ = v;
 return(v);
 }
@@ -281 +393 @@
    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. 
+   All the pixels are considered to have  the same size in (theta, phi).
+      
 */  
 template<class T>
 double LocalMap<T>::PixSolAngle(int_4 k) const
 {
+  double XP, YP, ZP;
   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);
+  PixToSphereC(i,j,XP,YP,ZP);
+
+  double theta = acos(ZP);
+
+  
 
   // angle solide
-  double sol= fabs((xa-x0)*(yb-y0)-(xb-x0)*(ya-y0));
+
+  double sol= 2.* deltaPhi_*sin(theta)*sin(0.5*deltaTheta_);
   return sol;
 }
 
+/*! \fn void SOPHYA::LocalMap::PixToSphereC(int_4 ip, int_4 it, double& XP, double& YP, double& ZP)
+
+projection of a pixel of map, onto the unity sphere ; result in cartesian coordinates.
+*/
+
+template<class T>
+void LocalMap<T>::PixToSphereC(int_4 ip, int_4 it, double& XP, double& YP, double& ZP) const
+{
+  double X= double(ip-x0_)* deltaPhi_;
+  double Y= double(it-y0_)*deltaTheta_; 
+
+    // situation dans le plan de reference tangent
+  double dx= X*cosAngle_-Y*sinAngle_;
+  double dy= X*sinAngle_+Y*cosAngle_;
+
+  XP = XC_ - dx*snphC_ - dy*csthC_*csphC_;
+  YP = YC_ + dx*csphC_ - dy*csthC_*snphC_;
+  ZP = ZC_ + dy*snthC_;
+
+  // on renormalise pour eviter les probleme de precision lors
+  // de la prise ulterieure de lignes trigonometriques
+  double norme = sqrt(XP*XP + YP*YP + ZP*ZP);
+  XP /= norme;
+  YP /= norme;
+  ZP /= norme;
+}
+
 /*! \fn void SOPHYA::LocalMap::SetOrigin(double theta0,double phi0,double angle)
 
@@ -324 +450 @@
 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;
+  SetOrigin(theta0, phi0,  nSzX_/2,  nSzY_/2, angle);
 }
  
@@ -342 +460 @@
 void LocalMap<T>::SetOrigin(double theta0,double phi0,int_4 x0,int_4 y0,double angle)
 {
-  theta0_= theta0;
-  phi0_  = phi0;
-  angle_ = angle;
+  // if (originFlag_)
+  //   {
+  //     throw PException("LocalMap::SetOrigin : redefining origin is not allowed at the moment ");
+  //  }
+
+  if ( theta0 < 0. || theta0 > 180. || phi0 < 0. || phi0 > 360. || angle < -90. || angle > 90.)
+    {
+      throw "LocalMap::SetOrigin  angle out of range"; 
+    }
+  SetCoorC(theta0,phi0);
+  angleDegres_ = angle; 
+  // en radians 
+  angle_ = angle*Pi/180.;
   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;
+  cosAngle_= cos(angle_);
+  sinAngle_= sin(angle_);
 } 
 
+template<class T>
+void LocalMap<T>::SetCoorC(double theta0, double phi0)
+{
+
+  thetaDegresC_ = theta0;
+  phiDegresC_ = phi0;
+
+  // passage en radians
+  thetaC_= theta0*Pi/180.;
+  phiC_  = phi0*Pi/180.;
+  csthC_ = cos(thetaC_);
+  snthC_ = sin(thetaC_);
+  csphC_ = cos(phiC_);
+  snphC_ = sin(phiC_);
+  XC_ = snthC_*csphC_;
+  YC_ = snthC_*snphC_;
+  ZC_ = csthC_;
+}
+
+
+// angles en RADIANS
+template<class T>
+TMatrix<double> LocalMap<T>::CalculMatricePassage()
+{
+  cout << " calcul matrice de passage " << endl;
+  TMatrix<double> passage(3,3);
+  double cos_th_axeZ;
+  double sin_th_axeZ;
+  double cos_phi_axeZ;
+  double sin_phi_axeZ;
+
+  double cth,sth, cdeltaPhi, phi, cphi, sphi;
+
+  if ( snthC_ <= 0.) // carte centree au pole 
+    {
+      cos_th_axeZ = 1.;
+      sin_th_axeZ = 0.;
+      cos_phi_axeZ = 1.;
+      sin_phi_axeZ = 0.;
+      
+    }
+  else
+    {
+      cth = cosAngle_*snthC_;
+      double arg = 1.-cth*cth;
+      if (arg <= 0. ) // carte centree sur l'equateur, sans rotation
+        {
+          cos_th_axeZ = 1.;
+          sin_th_axeZ = 0.;
+          cos_phi_axeZ = csphC_;
+          sin_phi_axeZ = snphC_;
+        }
+      else
+        {
+          sth = sqrt(arg);
+          cdeltaPhi = -csthC_*cth/(snthC_*sth);
+          if (cdeltaPhi < -1. ) cdeltaPhi = -1.;
+          if (cdeltaPhi > 1. ) cdeltaPhi = 1.;
+          phi = phiC_ + acos( cdeltaPhi );
+
+          cos_th_axeZ = cth;
+          sin_th_axeZ = sth;
+          cos_phi_axeZ = cos(phi);
+          sin_phi_axeZ = sin(phi);
+        }
+    }
+  passage(0,0) = snthC_*csphC_;
+  passage(1,0) = snthC_*snphC_;
+  passage(2,0) = csthC_;
+
+  passage(0,2) = sin_th_axeZ*cos_phi_axeZ;
+  passage(1,2) = sin_th_axeZ*sin_phi_axeZ;
+  passage(2,2) = cos_th_axeZ;
+
+  passage(0,1) = passage(1,2)*passage(2,0) - passage(2,2)*passage(1,0);
+  passage(1,1) = passage(2,2)*passage(0,0) - passage(0,2)*passage(2,0); 
+  passage(2,1) = passage(0,2)*passage(1,0) - passage(1,2)*passage(0,0);
+
+  //  passage.Print(cout);
+  // cout << " fin calcul passage " << endl;
+
+  return passage;
+
+
+}
+
 /*! \fn void SOPHYA::LocalMap::SetSize(double angleX,double angleY)
 
@@ -360 +572 @@
 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
+  if ( angleX <= 0. || angleX > 180. || angleY <= 0. || angleY > 180.)
+    {
+      throw "LocalMap::SetSize extension angle out of range"; 
+    }
+
+  angleDegresX_ = angleX;
+  angleDegresY_ = angleY;
+  // angles par pixel,  en radians 
+  cout << " taille x " << nSzX_ <<"  taille y " << nSzY_ << endl;
+  deltaPhi_   = angleX*Pi/(180.*nSzX_);
+  deltaTheta_ = angleY*Pi/(180.*nSzY_);
+
+
+}
+
+
+
+/*! \fn void SOPHYA::LocalMap::ProjectionToSphere(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;
-    }
-}
+void LocalMap<T>::ProjectionToSphere(SphericalMap<T>& sphere) const
+{
+   double theta, phi;
+  int it, ip;
+   for (it = 0; it < nSzY_; it++)
+    {
+      for (ip = 0; ip < nSzX_; ip++) 
+        {
+          PixThetaPhi(ip,it, theta, phi);
+          sphere(theta,phi)= pixels_(it, ip);
+        }
+    }
+}
+
+
+
 // Titre        Private Methods
 //++
@@ -391 +615 @@
 void LocalMap<T>::InitNul()
 //
-//    set some attributes to zero
+//    set 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$
+  nSzX_ = 0;
+  nSzY_ = 0;
+  angleDegresX_ = 0.;
+  angleDegresY_ = 0.;
+  thetaDegresC_ = 0.;
+  phiDegresC_   = 0.;
+  x0_ = 0;
+  y0_ = 0;
+  angleDegres_ = 0.;
+
+
+  nPix_ = 0;
+
+  thetaC_ = 0.;
+  phiC_   = 0.;
+  csthC_  = 1.;
+  snthC_  = 0.;
+  csphC_  = 1.;
+  snphC_  = 0.;
+  XC_     = 0.;
+  YC_     = 0.;
+  ZC_     = 1.;
+
+  angle_      = 0.;
+  cosAngle_   = 1.;
+  sinAngle_   = 0.;
+  deltaPhi_   = 0.;
+  deltaTheta_ =0.;
 }
 
@@ -414 +661 @@
 }
 
-/*! \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>
@@ -527 +670 @@
 {
   os<<" SzX= "<<nSzX_<<", SzY= "<<nSzY_<<", NPix= "<<nPix_<<endl;
-  if(LocalMap_isDone())
-    {
-      os<<" theta0= "<<theta0_<<", phi0= "<<phi0_<<", angle= "<<angle_<<endl;
+      os<<" theta0= "<<thetaC_*180./Pi <<", phi0= "<<phiC_*180./Pi <<", angle= "<<angle_*180./Pi<<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<<" cos= "<< cosAngle_ <<", & sin= "<<sinAngle_<<endl;
+      os<<"  deltaPhi_= "<<  deltaPhi_ <<", deltaTheta = "<< deltaTheta_ <<endl;
+      os <<endl;
 
   os << " contenu de pixels : ";
@@ -540 +680 @@
     {
       if(i%5 == 0) os << endl;
-      os <<  pixels_(i) <<", ";
+  int_4 ik,jk;
+  Getij(i,ik,jk);
+      os <<  pixels_(jk,ik) <<", ";
     }
   os << endl;
@@ -549 +691 @@
 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_;
+
+  nSzX_ = lm.nSzX_;
+  nSzY_ = lm.nSzY_;
+  nPix_ = lm.nPix_;
+  x0_   = lm.x0_;
+  y0_   = lm.y0_;
+
+  thetaC_ = lm.thetaC_;
+  phiC_   = lm.phiC_;
+  csthC_  = lm.csthC_;
+  snthC_  = lm.snthC_;
+  csphC_  = lm.csphC_;
+  snphC_  = lm.snphC_;
+  XC_     = lm.XC_;
+  YC_     = lm.YC_;
+  ZC_     = lm.ZC_;
+
   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_;
+  cosAngle_  = lm.cosAngle_;
+  sinAngle_  = lm.sinAngle_;
+  deltaPhi_  = lm. deltaPhi_;
+  deltaTheta_ = lm.deltaTheta_;
 }
 
@@ -654 +802 @@
     throw(SzMismatchError("LocalMap<T>::MulElt(const LocalMap<T>&) SizeMismatch")) ;
     }
-  pixels_ *= a.pixels_;
+  //  pixels_ *= a.pixels_;
+  pixels_.DataBlock() *= a.pixels_.DataBlock();
   return (*this);
 }
@@ -666 +815 @@
     throw(SzMismatchError("LocalMap<T>::DivElt(const LocalMap<T>&) SizeMismatch")) ;
     }
-  pixels_ /= a.pixels_;
+  //  pixels_ /= a.pixels_;
+  pixels_.DataBlock() /= a.pixels_.DataBlock();
   return (*this);
 }
@@ -673 +823 @@
 
 
-
-//++
-// Titre        class FIO_LocalMap
-//    Delegated objects for persitance management
-//--
 
 

trunk/SophyaLib/SkyMap/localmap.h

@@ -19 +19 @@
 
 
+template <class T>
+class FIO_LocalMap;
+
+template<class T>
+class FITS_LocalMap;
+
 /* Class LocalMap */
 
@@ -26 +32 @@
 {
 
+ // friend declaration for classes which handle persistence and FITS IO
+  friend class FIO_LocalMap<T>; 
+  friend class FITS_LocalMap<T>; 
+
 public:
 
 LocalMap();
-LocalMap(int_4 nx, int_4 ny);
+LocalMap(int_4 nx, int_4 ny, double angleX,double angleY, double theta0,double phi0,int_4 x0,int_4 y0,double angle=0.);
+LocalMap(int_4 nx, int_4 ny, double angleX,double angleY, double theta0,double phi0, double angle=0.);
 LocalMap(const LocalMap<T>& lm, bool share);
 LocalMap(const LocalMap<T>& lm);
@@ -69 +80 @@
 inline virtual char* TypeOfMap() const {return "LOCAL";};
  
-virtual void SetOrigin(double theta=90.,double phi=0.,double angle=0.);
-virtual void SetOrigin(double theta,double phi,int_4 x0,int_4 y0,double angle=0.);
-
-virtual void SetSize(double angleX,double angleY);
 
 /*! Check to see if the local mapping is done */
-inline bool LocalMap_isDone() const {return(originFlag_ && extensFlag_);};
-
-virtual void Project(SphericalMap<T>& sphere) const;
+  //inline bool LocalMap_isDone() const {return(originFlag_ && extensFlag_);};
+
+void PixThetaPhi(int_4 ip,int_4 it, double& theta,double& phi) const;
+
+
+void ProjectionToSphere(SphericalMap<T>&) const;
   
 /* There should be a more complex algorithm somewhere to combine *several* local maps to a full sphere.
@@ -92 +102 @@
 inline int_4 SizeY() const {return nSzY_;}
 
-inline void Origin(double& theta,double& phi,int_4& x0,int_4& y0,double& angle) const {theta= theta0_; phi= phi0_; x0= x0_; y0= y0_;angle= angle_;}
-
-inline void Aperture(double& anglex,double& angley) const {anglex= angleX_; angley= angleY_;}
+inline void Origin(double& theta,double& phi,int_4& x0,int_4& y0,double& angle) const {theta= thetaDegresC_; phi= phiDegresC_; x0= x0_; y0= y0_;angle= angleDegres_;}
+
+/*! total aperture in theta and phi, in degrees ( from SetSize() ) */
+inline void Aperture(double& anglex,double& angley) const {anglex= angleDegresX_; angley= angleDegresY_;}
 
 
 /*  Acces to the DataBlock  */
-inline       NDataBlock<T>& DataBlock()       {return pixels_;}
-inline const NDataBlock<T>& DataBlock() const {return pixels_;}
+inline       NDataBlock<T>& DataBlock()       {return pixels_.DataBlock();}
+inline const NDataBlock<T>& DataBlock() const {return pixels_.DataBlock();}
+
+/*  Acces to the matrix  */
+//! access to matrix 
+/*!
+  \warning : a pixel is defined by the pair of a phi-like index 
+   (x-axis) and a theta-like index (y-axis). The phi-like index denotes
+   the column number in the matrix ; the theta-like index denotes the
+   row number.
+*/
+inline       TMatrix<T>& Matrix()       {return pixels_;}
+inline const TMatrix<T>& Matrix() const {return pixels_;}
 
 /* impression */ 
@@ -156 +178 @@
 
 void InitNul();
+void SetOrigin(double theta=90.,double phi=0.,double angle=0.);
+void SetOrigin(double theta,double phi,int_4 x0,int_4 y0,double angle=0.);
+void SetSize(double angleX,double angleY);
+void SetCoorC(double theta0, double phi0);
+TMatrix<double> CalculMatricePassage();
 void Getij(int_4 k,int_4& i,int_4& j) const;
-void ReferenceToUser(double& theta,double& phi) const; 
-void UserToReference(double& theta,double& phi) const;
-void PixProjToAngle(double x,double y,double& theta,double& phi) const; 
-void AngleProjToPix(double theta,double phi,double& x,double& y) const; 
+  void PixToSphereC(int_4 ip, int_4 it, double& XP, double& YP, double& ZP) const;
 
 void recopierVariablesSimples(const LocalMap<T>& lm);
@@ -167 +191 @@
 // ---------- Variables internes ----------------------------
 
-int_4 nSzX_;
-int_4 nSzY_;
-int_4 nPix_;
-bool originFlag_;
-bool extensFlag_;
-int_4 x0_;
-int_4 y0_;
-double theta0_;
-double phi0_;
-double angle_;
-double cos_angle_;
-double sin_angle_;
-double angleX_;
-double angleY_;
-double tgAngleX_;
-double tgAngleY_;
-NDataBlock<T> pixels_;
+
+  // variables suffisantes pour reconstruire l'objet
+
+  int_4 nSzX_;
+  int_4 nSzY_;
+  double angleDegresX_;
+  double angleDegresY_;
+  double thetaDegresC_;
+  double phiDegresC_;
+  int_4 x0_;
+  int_4 y0_;
+  double angleDegres_;
+
+  //  NDataBlock<T> pixels_;
+  TMatrix<T> pixels_;
+
+
+  // variables derivees (redondantes, precalculees pour ameliorer 
+  // les performances)
+
+  int_4 nPix_;
+
+
+  double thetaC_;
+  double phiC_;
+  double csthC_;
+  double snthC_;
+  double csphC_;
+  double snphC_;
+  double XC_;
+  double YC_;
+  double ZC_;
+
+
+  double angle_;
+  double cosAngle_;
+  double sinAngle_;
+  double deltaPhi_;
+  double deltaTheta_;
+
+
+
 };