(JEC) 25/4/15 the 2D-spherical-geometry management files

lagsht_spheregeom.cc

+
+#include "lagsht_spheregeom.h" 
+#include "lagsht_exceptions.h"
+#include "lagsht_healpixhelper.h"
+
+namespace LagSHT {
+
+
+//------------------------- ECP/Fejer1 Geom ---------------------
+
+  int ECPGeometry::PixIndexSph(r_8 theta, r_8 phi) const {
+
+    int stride_phi=1;         //this migth be also a memeber data
+    int stride_theta=Nphi_;  //"
+    r_8 phi0 = 0.; //for the time beeing
+
+    int itheta=int(Nrings_*theta/M_PI);
+    int iphi = int(Nphi_+ Nphi_*(phi-phi0)/(2*M_PI)+0.5)%Nphi_;
+
+    return stride_phi*iphi+stride_theta*itheta;
+
+  }//PixIndexSph
+
+  void ECPGeometry::PixThetaPhi(int k, r_8& theta, r_8& phi) const {
+
+    //JEC if is not the case change the code       int stride_phi=1;         //this migth be also a memeber data
+    //JEC if is not the case change the code      int stride_theta=Nphi_;  //"
+    r_8 phi0 = 0.;
+
+    int iphi=k%Nphi_;
+    int itheta=k/Nphi_;
+    phi=phi0 + iphi*2*M_PI/Nphi_;
+    theta=(itheta+0.5)*M_PI/Nrings_;
+    
+  }//PixThetaPhi
+
+  void ECPGeometry::SetMap(int nrings, int nphi) {
+
+      geomType_ = "ECP";
+
+      if(geomDone_){
+	cout << "-------------------------------------------------" <<endl;
+	cout << "---------------  New ECP Geometry  ------------- " <<endl;
+	cout << "-------------------------------------------------" <<endl;
+	geomDone_ = false;
+      }
+      
+      //dans sharp_testsuite.c Fejer's 1st rule
+      //nrings = 2Lmax+1
+      //npixper ring = 2Mmax+1
+      
+      Nrings_ = (nrings == -1) ? 2*L_-1 : nrings;
+      Nphi_   = (nphi == -1) ?  2*M_-1 : nphi;
+      
+      sharpjob_.set_ECP_geometry(Nrings_,Nphi_);
+      
+      
+      //do the common features
+      SetCommonMapFeatures();
+
+      geomDone_=true;
+
+  }//SetMap
+  
+//------------------------- Gauss Geom ---------------------
+
+int GaussGeometry::PixIndexSph(r_8 theta, r_8 phi) const {
+  int stride_phi=1;         //this migth be also a memeber data
+    int stride_theta=Nphi_;  //"
+    r_8 phi0 = 0.; //for the time beeing
+
+    int itheta = std::min_element(thetaPix_.begin(),thetaPix_.end(), closer_to(theta)) - thetaPix_.begin();
+    int iphi = int(Nphi_+ Nphi_*(phi-phi0)/(2*M_PI)+0.5)%Nphi_;
+    
+    return stride_phi*iphi+stride_theta*itheta;
+
+}//PixIndexSph
+
+void GaussGeometry::PixThetaPhi(int k, r_8& theta, r_8& phi) const {
+  //Nb this is the same code as ECPGeometry
+//JEC if is not the case change the code       int stride_phi=1;         //this migth be also a memeber data
+//JEC if is not the case change the code      int stride_theta=Nphi_;  //"
+      r_8 phi0 = 0.;
+
+      int iphi=k%Nphi_;
+      int itheta=k/Nphi_;
+      phi=phi0 + iphi*2*M_PI/Nphi_;
+      theta=(itheta+0.5)*M_PI/Nrings_;
+
+}//PixThetaPhi
+
+
+void GaussGeometry::SetMap(int nrings, int nphi) {
+      geomType_ = "Gauss";
+
+      if(geomDone_){
+	cout << "-------------------------------------------------" <<endl;
+	cout << "---------------  New Gauss Geometry  ------------- " <<endl;
+	cout << "-------------------------------------------------" <<endl;
+	geomDone_ = false;
+      }
+
+      //dans sharp_testsuite.c Gauss
+      //nrings = Lmax+1
+      //npixper ring = 2Mmax+1
+
+      Nrings_ = (nrings == -1) ? L_ : nrings;
+      Nphi_   = (nphi == -1) ?  2*M_-1 : nphi;
+
+      sharpjob_.set_Gauss_geometry(Nrings_,Nphi_);
+
+      //do the common features
+      SetCommonMapFeatures();
+
+      geomDone_=true;
+
+
+}//SetMap
+
+
+
+//------------------------- HEALPIX Geom ---------------------
+
+  int HealpixGeometry::PixIndexSph(r_8 theta, r_8 phi) const {
+
+    if(theta<=0 || theta>=M_PI)  LagSHTError("BaseLagSphTransform::PixIndexSph: Healpix theta out of range");
+    return HEALPIX::ang2pix_ring_z_phi(Nrings_,cos(theta),phi);
+
+  }//PixIndexSph
+
+  void HealpixGeometry::PixThetaPhi(int k, r_8& theta, r_8& phi) const {
+
+    r_8 z;
+    HEALPIX::pix2ang_ring_z_phi (Nrings_,k,&z,&phi);
+    theta = acos(z);
+
+  }//PixThetaPhi
+
+  void HealpixGeometry::SetMap(int nsides, int DUMMY) { //second argument not used
+
+      geomType_ = "Healpix";
+
+      if(geomDone_){
+	cout << "-------------------------------------------------" <<endl;
+	cout << "---------------  New Healpix Geometry  --------- " <<endl;
+	cout << "-------------------------------------------------" <<endl;
+	geomDone_ = false;
+      }
+      
+      //dans sharp_testsuite.c Healpix
+      //nside = Lmax/2;
+      Nrings_ = (nsides == -1 ) ? (L_ -1)/2: nsides; // this is the nside argument in fact
+      Nphi_ = -1; // is not used
+      
+      sharpjob_.set_Healpix_geometry(Nrings_);
+
+      //do the common features
+      SetCommonMapFeatures();
+
+      geomDone_=true;
+
+  }//SetMap
+
+  //---------------------------------------
+
+  void BaseGeometry::SetCommonMapFeatures() {
+  
+      
+  
+  sharpjob_.set_triangular_alm_info(L_-1,M_-1); 
+  
+  Npix_ = get_npix(sharpjob_.get_geom_info());    
+  Nalm_ = get_nalms(sharpjob_.get_alm_info());
+
+  //fill the ANN-kdTree to prepare the filling of the 2D spheres
+  //  geom_t vecPix;
+  SetPixelPositions();
+
+    
+  cout << "SetThetaPhiMap: Geom, Nring, Nphi, Npix, Nalm: " << geomType_ << ", "
+       << Nrings_ << ", " 
+       << Nphi_ << ", "
+       << Npix_ << ", "
+       << Nalm_ 
+       << endl;
+
+}//SetThetaPhiMap
+
+
+
+  //----------------------------------------------------------
+
+void BaseGeometry::SetPixelPositions() {
+  const sharp_geom_info *ginfo; 
+  ginfo = sharpjob_.get_geom_info();
+  for (int i=0; i<ginfo->npairs; ++i) {
+      double theta1 = ginfo->pair[i].r1.theta;
+      thetaPix_.push_back(theta1);
+      // and so on for nph, phi0
+      int nph1 = ginfo->pair[i].r1.nph;
+      double phi0 = ginfo->pair[i].r1.phi0;
+      for(int j=0; j<nph1; j++) {
+	double phi1 = phi0 + j*(2.0*M_PI/nph1);
+	//	cout << "vecPix_s ( "<< theta1 << ", " << phi1 << ")" << endl;
+	vecPix_.push_back(make_pair(theta1,phi1));
+      }
+      if (ginfo->pair[i].r2.nph>0) {// the second ring in this pair exists
+	double theta2 = ginfo->pair[i].r2.theta;
+	thetaPix_.push_back(theta2);
+	int nph2 = ginfo->pair[i].r2.nph;
+	double phi0 = ginfo->pair[i].r2.phi0;
+	for(int j=0; j<nph2; j++) {
+	  double phi2 = phi0 + j*(2.0*M_PI/nph2);
+	  //	  cout << "vecPix_s ( "<< theta2 << ", " << phi2 << ")" << endl;
+	  vecPix_.push_back(make_pair(theta2,phi2));
+	}
+      }
+  }
+
+  sort(vecPix_.begin(),vecPix_.end()); //See if it is necessary...
+  sort(thetaPix_.begin(), thetaPix_.end()); //mandatory for PixIndexSph
+
+}//SetPixelPositions
+
+
+}// Fin de namespace

lagsht_spheregeom.h

+#ifndef LAGSHTSPHEREGEOM_SEEN
+#define LAGSHTSPHEREGEOM_SEEN
+
+#include <cmath>
+#include <string>
+#include <iostream>
+#include <algorithm>
+#include <utility>      // std::pair, std::make_pair
+using namespace std;
+
+#include "lagsht_geom.h"
+//Sharp
+#include "sharp_lowlevel.h"
+#include "sharp_geomhelpers.h"
+#include "sharp_almhelpers.h"
+#include "sharp_cxx.h"
+
+namespace LagSHT {
+
+class BaseGeometry {
+ public:
+
+  /*! Creator
+    \input Lmax : maximal l value: l:0,..., Lmax-1
+    \input Mmax : maximal m value by default it is Lmax
+    \input Nrings : in general the total number of iso-lattitude rings (except for Healpix)
+    \input Nphi   : in general the total number of phi equidistant pixel in a ring
+   */
+  BaseGeometry(int Lmax, int Mmax, int Nrings, int Nphi):
+    L_(Lmax), M_(Mmax), 
+    Nrings_(Nrings), Nphi_(Nphi), 
+    Npix_(0), Nalm_(0), geomType_(""), geomDone_(false) {
+
+     if (Mmax<0) M_=Lmax;
+              
+  }
+
+
+  //! Destructor
+  virtual ~BaseGeometry() {}
+
+  //! Number of colatitude rings of the 2D-sphere representation
+  int Nrings() const { return Nrings_; }
+  //! Side parameter for Healpix
+  int Nsides() const {  return Nrings_; } //this is a convention here
+  //! Number of colatitude rings of the 2D-sphere representation
+  int NTheta() const { return Nrings(); }
+  //! Number of phi values (pixels) per colatitude rings
+  int NPhi()   const { return Nphi_; }
+  //! Total number of 2D pixels
+  int Npix()   const { return Npix_; }
+  //! Total number of Alm coefficients
+  int Nalm()   const { return Nalm_; }
+  //! Type of geomery
+  string Geometry() const { return geomType_; }
+  
+  //! Sorted List of pixel center (theta, phi) positions
+  geom_t GetPixelPositions() const { return vecPix_; }
+
+   
+  /*! Get Pixel index in the 2D map
+    \input theta in [0, pi]
+    \input phi   in [0, 2pi]
+  */
+  virtual int PixIndexSph(r_8 theta, r_8 phi) const = 0;
+
+  /*! Get Theta Phi position of a given pixel index
+    \input k in [0, Npix-1] no boundchckeing done
+   */
+  virtual void PixThetaPhi(int k, r_8& theta, r_8& phi) const = 0;
+
+
+  /*! Choose the 2D map representation & initialization of libsharp
+    \input ECP/Fejer1, Gauss, and Healpix
+    \input nrings
+    \input nphi
+   */
+  virtual void SetMap(int nrings=-1, int nphi=-1) = 0;
+
+  
+  virtual void alm2map(const r_8 *alm, r_8 *map, bool add=false){
+    sharpjob_.alm2map(alm,map,add);
+  }
+  virtual void map2alm(const r_8 *map, r_8 *alm, bool add=false){
+    sharpjob_.map2alm(map,alm,add);
+  }
+
+protected:
+
+  //! Set the vector of pixel positions
+  void SetPixelPositions();
+
+  //! Manage the commonality of the different geometries 
+  void SetCommonMapFeatures();
+
+  //! Get the total number of pixels (from sharp test suite)
+  static ptrdiff_t get_npix(const sharp_geom_info *ginfo) {
+    ptrdiff_t res=0;
+    for (int i=0; i<ginfo->npairs; ++i){
+      res += ginfo->pair[i].r1.nph;
+      if (ginfo->pair[i].r2.nph>0) res += ginfo->pair[i].r2.nph;
+    }
+    return res;
+  }
+  
+  //from sharp test suite
+  static ptrdiff_t get_nalms(const sharp_alm_info *ainfo) {
+    ptrdiff_t res=0;
+    for (int i=0; i<ainfo->nm; ++i)
+      res += ainfo->lmax-ainfo->mval[i]+1;
+    return res;
+  }
+
+ 
+
+protected:
+  int L_; //!< Spherical harmonic limit such that l:0,...,L-1   m:0,...,l
+  int M_; //!< maximum value of m
+//   int N_; //!< Radial band limit such that n:0,...,N-1
+  int Nrings_; //!< Number of colatitude rings of the 2D-sphere representation
+  int Nphi_; //!< Number of phi values (pixels) per colatitude rings
+
+//   int alpha_; //!< parameter of the generalized Laguerre Polynomials
+
+  int Npix_; //!< Total number of 2D pixels
+  int Nalm_; //!< Total number of Alm coefficients
+
+  geom_t vecPix_; //!< sorted vector of pixel positions
+  vector<r_8> thetaPix_; //!< sorted vector of theta rings
+  string geomType_; //!<< Type of pixelization
+  bool geomDone_; //!< Set to true once the geometry is properly set (See SetThetaPhiMap(string) )
+
+  sharp_cxxjob<r_8> sharpjob_; //!< libsharp job // voir si r_8 -> r_16 extend the validity of the computation
+
+
+}; //BaseGeometry
+
+  //-----------------------------------------------------
+
+class ECPGeometry : public BaseGeometry {
+  
+ public:
+  //! Constructor
+  ECPGeometry(int Lmax, int Mmax, int Nrings, int Nphi):
+    BaseGeometry(Lmax, Mmax, Nrings, Nphi) {}
+
+  //! Destructor 
+  virtual ~ECPGeometry() {}
+
+  /*! Get Pixel index in the 2D map
+    \input theta in [0, pi]
+    \input phi   in [0, 2pi]
+  */
+  virtual int PixIndexSph(r_8 theta, r_8 phi) const;
+
+  /*! Get Theta Phi position of a given pixel index
+    \input k in [0, Npix-1] no boundchckeing done
+   */
+  virtual void PixThetaPhi(int k, r_8& theta, r_8& phi) const;
+
+  /*! Choose the 2D map representation & initialization of libsharp
+    \input ECP/Fejer1
+    \input nrings
+    \input nphi
+   */
+  virtual void SetMap(int nrings=-1, int nphi=-1);
+
+
+}; //ECPGeometry
+  //-----------------------------------------------------
+
+//alias of different geometry
+#define Fejer1Geometry ECPGeometry
+
+  //-----------------------------------------------------
+
+class GaussGeometry : public BaseGeometry {
+  
+ public:
+  //! Constructor
+  GaussGeometry(int Lmax, int Mmax, int Nrings, int Nphi):
+    BaseGeometry(Lmax, Mmax, Nrings, Nphi) {}
+
+  //! Destructor 
+  virtual ~GaussGeometry() {}
+
+  /*! Get Pixel index in the 2D map
+    \input theta in [0, pi]
+    \input phi   in [0, 2pi]
+  */
+  virtual int PixIndexSph(r_8 theta, r_8 phi) const;
+
+  /*! Get Theta Phi position of a given pixel index
+    \input k in [0, Npix-1] no boundchckeing done
+   */
+  virtual void PixThetaPhi(int k, r_8& theta, r_8& phi) const;
+
+
+  /*! Choose the 2D map representation & initialization of libsharp
+    \input Gauss
+    \input nrings
+    \input nphi
+   */
+  virtual void SetMap(int nrings=-1, int nphi=-1);
+
+protected:
+
+  /*! Utility function object to perform PixIndexSph theta search
+   */
+  struct closer_to {
+    r_8 target_;
+    explicit closer_to( r_8 target ) : target_( target ){}
+
+    bool operator ()( r_8 left, r_8 right ) const
+    {
+        return std::abs( target_ - left ) < std::abs( target_ - right );
+    }
+  };
+ 
+
+}; //GaussGeometry
+  //-----------------------------------------------------
+class HealpixGeometry : public BaseGeometry {
+  
+ public:
+  /*! Constructor
+    \input Lmax 
+    \input Mmax 
+    \input Nsides (= Nrings parameter of BaseGeometry) 
+   */
+  HealpixGeometry(int Lmax, int Mmax, int Nsides):
+    BaseGeometry(Lmax, Mmax, Nsides,-1) {}
+
+  //! Destructor 
+  virtual ~HealpixGeometry() {}
+
+  /*! Get Pixel index in the 2D map
+    \input theta in [0, pi]
+    \input phi   in [0, 2pi]
+  */
+  virtual int PixIndexSph(r_8 theta, r_8 phi) const;
+
+  /*! Get Theta Phi position of a given pixel index
+    \input k in [0, Npix-1] no boundchckeing done
+   */
+  virtual void PixThetaPhi(int k, r_8& theta, r_8& phi) const;
+
+
+  /*! Choose the 2D map representation & initialization of libsharp
+    \input Healpix
+    \input nsides
+   */
+  virtual void SetMap(int nsides=-1, int DUMMY=-1); //2nd arg not used but here for compatibility
+
+
+}; //HealpixGeometry
+
+
+}//end namespace
+#endif