(JEC) 9/4/15 helper to get the pixels (theta, phi) positions ordered

README.md

@@ -10,7 +10,8 @@ git clone git@gitlab.in2p3.fr:campagne/LagSHT.git
 # 9/4/15
 * new set of files for a Stand-a-lone version wo `Sophya`
 * extension of the computation of the weights (Laguerre function) up to N ~ 5500 by transfert of the exponential scale a the very beginning (see laguerreBuilder.cc)
-*  Start documentation (Doxygen)
+* setup some helper to get the geometry (theta, phi) positions of geometry 
+* start documentation (Doxygen)
 
 `> doxygen doxydoc`  
 

lagSphericTransform.cc

 #include <stdlib.h>
 #include <iostream>
+#include <algorithm> //copy, sort
+#include <iterator> // for ostream_iterator
 
 #include "lagSphericTransform.h"
 #include "lagsht_exceptions.h"
@@ -39,8 +41,40 @@ void BaseLagSphTransform::SetThetaPhiMap(string choice){
        << Npix_ << ", "
        << Nalm_ 
        << endl;
+  
 }//SetThetaPhiMap
 
+BaseLagSphTransform::geom_t BaseLagSphTransform::GetPixelPositions() {
+  geom_t pixels(Npix_);
+  sharp_geom_info *ginfo; 
+  ginfo = sharpjob_.get_geom_info();
+  for (int i=0; i<ginfo->npairs; ++i) {
+      double theta1 = ginfo->pair[i].r1.theta;
+      // 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 << "pixels ( "<< theta1 << ", " << phi1 << ")" << endl;
+	pixels.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;
+	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 << "pixels ( "<< theta2 << ", " << phi2 << ")" << endl;
+	  pixels.push_back(make_pair(theta2,phi2));
+	}
+      }
+  }
+
+  sort(pixels.begin(),pixels.end());
+  
+  return pixels;
+}//GetPixelPositions
+
 
 
 void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& flmn, 

lagSphericTransform.h

 #ifndef LAGSPHERICTRANSFORM_SEEN
 #define LAGSPHERICTRANSFORM_SEEN
+
+#include <utility>      // std::pair, std::make_pair
+using namespace std;
+
 //Sharp
 #include "sharp_lowlevel.h"
 #include "sharp_geomhelpers.h"
@@ -21,7 +25,7 @@ namespace LagSHT {
 class BaseLagSphTransform {
  public:
 
-  typedef pair<r_8,r_8> pixpos_t;
+  typedef pair<r_8,r_8>    pixpos_t;
   typedef vector<pixpos_t> geom_t;
 
   //! Creator
@@ -45,14 +49,12 @@ class BaseLagSphTransform {
   int Nalm()   { return Nalm_; }
   
   //! List of pixel positions
-  geom_t GetPixelPositions() {
-    throw LagSHTError("BaseLagSphTransform::GetPixelPositions not yet implemented");
-  }
-  
-  
+  geom_t GetPixelPositions();
+    
   //! Choose the 2D map representation & initialization of libsharp
   void SetThetaPhiMap(string choice ="");
 
+
 protected:
 
   //from sharp test suite