Skip to content

GitLab

Commit 2c360460 authored by Jean-Eric Campagne's avatar Jean-Eric Campagne
Browse files

(JEC) 5/5/15 imporve Doc and clean code

parent cb625b23
Hide whitespace changes
Inline Side-by-side
Showing with 178 additions and 190 deletions
Documentation.md
View file @ 2c360460
......@@ -58,21 +58,21 @@ The analysis process starts by gathering the 3D pixels real values \f$f_{ijk} \e
are the 2D angles defined by a pixelization scheme.
For instance, one can use the ECP Fejer's sum rule with \f$N_\theta = 2L-1\f$ colatitude rings of \f$N_\phi \geq 2L-1\f$ iso-latitude pixels.
Then, one proceeds as followed
+ for the \f$k\f$th-radial shell, one performs a Spherical Harmonic decomposition using for instance the `map2alm` routine of `libsharp` to get the set of complex coefficients \f$f_{lmk}\f$ with \f$l\in\{0,\dots,L-1\}\f$ and \f$m \in\{0,\dots,l\}\f$;
+ for each \f$(l,m)\f$, one uses the set of values \f$f_{lmk}\f$ with \f$k\in\{0,\dots,N-1\}\f$ to determine the \f$f_{lmn}\f$ complex coefficients thanks to the following simplified relation
+ for the \f$k\f$th-radial shell, one performs a Spherical Harmonic decomposition using for instance the `map2alm` routine of `libsharp` to get the set of complex coefficients \f$a_{lmk}=a_{lm}(r_k)\f$ with \f$l\in\{0,\dots,L-1\}\f$ and \f$m \in\{0,\dots,l\}\f$ (the notation remind the underlying 2D Spherical Harmonic transform and the usual way to note the resulting coefficients \f$a_{lm}\f$);
+ for each \f$(l,m)\f$, one uses the set of values \f$a_{lmk}\f$ with \f$k\in\{0,\dots,N-1\}\f$ to determine the \f$f_{lmn}\f$ complex coefficients thanks to the following simplified relation
\f[
f_{lmn} = \sum_{k=0}^{N-1}w_k\ f_{lmk}\ {\cal L}_n(r_k) \quad\quad n\in\{0,\dots,N-1\}
f_{lmn} = \sum_{k=0}^{N-1}w_k\ a_{lmk}\ {\cal L}_n(r_k) \quad\quad n\in\{0,\dots,N-1\}
\f]
+ **Synthesis**:
The synthesis also proceeds in two steps as the above detailed analysis. The algorithm receives as input the \f$f_{lmn}\f$ complex coefficients
+ to first compute for each \f$(l,m)\f$ couple of indices the intermediate \f$f_{lmk}\f$ complex values obtained from the Inverse Laguerre Transform:
+ to first compute for each \f$(l,m)\f$ couple of indices the intermediate \f$a_{lmk}\f$ complex values obtained from the Inverse Laguerre Transform:
\f[
f_{lmk} = \sum_{n=0}^{N-1} f_{lmn}\ {\cal L}_n^{(\alpha)}(r_k)
a_{lmk} = \sum_{n=0}^{N-1} f_{lmn}\ {\cal L}_n^{(\alpha)}(r_k)
\f]
+ then, thanks to the Inverse Spherical Harmonic Transform, `alm2map` routine of `libsharp`, ones determines the real \f$f_{ijk}\f$ pixel values.
We take advantage from the Matrix Multiplication writing of the \f$f_{lmn} \leftrightarrow a_{lmk}\f$ passage to use efficient algorithm or even more efficient the BLAS-like libraries (OpenBLAS for Linux and the native Accelerate framework on Mac OS X).
The main classes in brief
------------------------
......
doxydoc
View file @ 2c360460
......@@ -1552,7 +1552,7 @@ EXTRA_SEARCH_MAPPINGS =
# If the GENERATE_LATEX tag is set to YES doxygen will generate LaTeX output.
# The default value is: YES.
GENERATE_LATEX = NO
GENERATE_LATEX = YES
# The LATEX_OUTPUT tag is used to specify where the LaTeX docs will be put. If a
# relative path is entered the value of OUTPUT_DIRECTORY will be put in front of
......
fig-schema.png
View replaced file @ cb625b23 View file @ 2c360460
fig-schema.png

21.9 KB | W: | H:

fig-schema.png

21.8 KB | W: | H:

fig-schema.png
fig-schema.png
fig-schema.png
fig-schema.png
  • 2-up
  • Swipe
  • Onion skin
lagsht_spheregeom.cc
View file @ 2c360460
......@@ -5,222 +5,210 @@
namespace LagSHT {
//----------------- Base Class of Geometry ----------------------
//------------------------- 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
void BaseGeometry::SetCommonMapFeatures() {
//sharp initialization
sharpjob_.set_triangular_alm_info(L_-1,M_-1);
Npix_ = get_npix(sharpjob_.get_geom_info());
Nalm_ = get_nalms(sharpjob_.get_alm_info());
int itheta=int(Nrings_*theta/M_PI);
int iphi = int(Nphi_+ Nphi_*(phi-phi0)/(2*M_PI)+0.5)%Nphi_;
SetPixelPositions();
cout << "SetThetaPhiMap: Geom, Nring, Nphi, Npix, Nalm: " << geomType_ << ", "
<< Nrings_ << ", "
<< Nphi_ << ", "
<< Npix_ << ", "
<< Nalm_
<< endl;
}//SetThetaPhiMap
return stride_phi*iphi+stride_theta*itheta;
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
}//PixIndexSph
void ECPGeometry::PixThetaPhi(int k, r_8& theta, r_8& phi) const {
//------------------------- ECP/Fejer1 Geom ---------------------
//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 ECPGeometry::PixIndexSph(r_8 theta, r_8 phi) const {
int iphi=k%Nphi_;
int itheta=k/Nphi_;
phi=phi0 + iphi*2*M_PI/Nphi_;
theta=(itheta+0.5)*M_PI/Nrings_;
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 {
}//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_);
//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
//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;
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_;
//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;
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() {
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);
sharpjob_.set_triangular_alm_info(L_-1,M_-1);
}//PixIndexSph
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 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 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));
}
}
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
sort(vecPix_.begin(),vecPix_.end()); //See if it is necessary...
sort(thetaPix_.begin(), thetaPix_.end()); //mandatory for PixIndexSph
}//SetPixelPositions
}// Fin de namespace
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment