lagsht_spheregeom.cc

#include "lagsht_spheregeom.h" 
#include "lagsht_exceptions.h"
#include "lagsht_healpixhelper.h"

namespace LagSHT {

//----------------- Base Class of Geometry ----------------------

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());

  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


//------------------------- 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



}// Fin de namespace