lagSphericTransform.h 5.89 KB
Newer Older
1 2
#ifndef LAGSPHERICTRANSFORM_SEEN
#define LAGSPHERICTRANSFORM_SEEN
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
/*
 *  This file is part of LagSHT.
 *
 *  LagSHT is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  LagSHT is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with libsharp; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 */
20

21 22 23 24 25
/*
 *  LagSHT is being developed at the Linear Accelerateur Laboratory (LAL)
 *  91898 ORSAY CEDEX - FRANCE
 *  main author: J.E Campagne
 */
26

27

28
#include "lagsht_spheregeom.h"
29 30 31 32 33 34 35 36 37 38
#include "laguerreTransform.h"
#include "lagsht_exceptions.h"


namespace LagSHT {


////////////////////////////////////////////////////////////////
//// --------- Class Spherical Harmonic Laguerre Transform -- //
////////////////////////////////////////////////////////////////
39 40 41 42 43 44 45 46 47 48 49
/* Example to access the Laguerre Coeff flmn and the intermediate alm(r_k) Coeff.
  for(int n=0;n<Nmax;n++){
    for(int l=0;l<Lmax;l++){
      for (int m=0;m<=l;m++) {
	int id= n*Nalm + l+m*Lmax-m*(m+1)/2; // LagSHT numbering 
	r_8 flmn2 = flmn[id].real()*flmn[id].real() + flmn[id].imag()*flmn[id].imag();
	r_8 almk2 = almk[id].real()*almk[id].real() + almk[id].imag()*almk[id].imag();
      }
    }
  }
 */
50
class LaguerreSphericalTransform {
51 52
 public:

53

54 55
  LaguerreSphericalTransform(string spheregeom, int Lmax, int Mmax, int Nmax, 
			     r_8 R, int Nrings = -1, int Nphi = -1, int alpha=2): 
56
    Lmax_(Lmax), Mmax_(Mmax), R_(R) {
57 58 59
    
    //Factory
    if (spheregeom == "ECP" || spheregeom == "Fejer1") {
60
      sphere_ = new ECPGeometry(Lmax, Mmax, Nrings, Nphi);
61
    } else if (spheregeom == "Gauss") {
62
      sphere_ = new GaussGeometry(Lmax, Mmax, Nrings, Nphi);      
63
    } else if (spheregeom == "Healpix") {
64
      sphere_ = new HealpixGeometry(Lmax, Mmax, Nrings); //here Nsides = Nrings
65 66 67 68 69
    } else {
      cout << "FAILURE: LaguerreSphericalTransform::Ctor wrong geometry type" << endl;
      sphere_ = 0;
      throw LagSHTError("LaguerreSphericalTransform::Ctor wrong geometry type");
    }
70 71 72

    lagTrans_ = new LaguerreTransform(Nmax, R, alpha); 

73 74 75 76 77 78
    
    //do the map pixelization
    if(sphere_) {
      //store localy usefull parameters for the algorithms Analysis/Synthesis
      Npix_ = sphere_->Npix();
      Nalm_ = sphere_->Nalm();
79 80
    } else {
      throw LagSHTError("LaguerreSphericalTransform::Ctor sphere_ ptr NULL");
81 82
    }

83 84 85 86 87 88 89 90 91

    if(lagTrans_) {
      N_    = lagTrans_->GetOrder();
      alpha_= lagTrans_->GetAlpha();
    } else {
      throw LagSHTError("LaguerreSphericalTransform::Ctor lagtrans_ ptr NULL");
    }


92 93 94 95
  }//Ctor

    
  //! Destructor
96
  virtual ~LaguerreSphericalTransform() {
97 98
    if(sphere_) delete sphere_;  sphere_ = 0;
    if(lagTrans_) delete lagTrans_; lagTrans_ = 0;
99
  }
100
  
101 102 103
  //! Geometry helper
  BaseGeometry* GetSphericalGeometry() const { return sphere_;}

104 105 106
  //! Laguerre Transform helper
  LaguerreTransform* GetLagTransform() const { return lagTrans_;}

107 108


109
   //! Getters
110 111 112 113
  int GetOrder()   const { return N_; }
  int GetAlpha()   const { return alpha_; }
  int GetNPixels() const { return Npix_; }
  int GetNAlm()    const { return Nalm_; }
114 115


116 117 118 119 120 121 122 123 124 125 126 127 128 129
  //! Get the f(l,m,n) coefficient in the vector collection
  //make a Flmn Class
  inline complex<r_8> GetFlmn(const vector< complex<r_8> >& flmn, int l , int m, int n) const {
    int id= n*Nalm_ + l +m*Lmax_ -m*(m+1)/2;
    return flmn[id];
 }
  //! Set the f(l,m,n) coefficient in the vector collection
  inline void SetFlmn(vector< complex<r_8> >& flmn, int l , int m, int n, complex<r_8> value) const {
     int id= n*Nalm_ + l +m*Lmax_ -m*(m+1)/2;
     flmn[id] = value;
  }
  


130 131
 //! Synthesis 
 /*! \brief Coeffs -> Pixels with Input/Output using T floating representation
132 133
   \input flmn: 3D complex spherical-laguerre coefficients
   \output fijk: 3D real function values
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
134
   \output almk: 2D alm complex spherical coefficients for each subshell k
135
 */
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
136 137 138 139 140 141 142 143 144 145 146 147 148 149
  void Synthesis(const vector< complex<r_8> >& flmn, 
		 vector<r_8>& fijk, 
		 vector< complex<r_8> >& almk);

  //backward compatibility
  void Synthesis(const vector< complex<r_8> >& flmn, vector<r_8>& fijk) {
    int NtotCoeff = Nalm_*N_;
    vector< complex<r_8> > almk(NtotCoeff);
    Synthesis(flmn,fijk,almk);
  }




150 151 152

  //! Analysis
  /*! \brief Pixels -> Coeffs with Input/Output using T floating representation
153 154
    \input fijk: 3D real function values
    \output flmn: 3D complex spherical-laguerre coefficients
155
  */
156 157 158 159 160 161 162 163 164 165
  void Analysis(const vector<r_8>& fijk, vector< complex<r_8> >& flmn){
    int NtotCoeff = Nalm_*N_;
    vector< complex<r_8> > almk(NtotCoeff);
    Analysis(fijk,flmn,almk);
  }



  //! Analysis
  /*! \brief Pixels -> Coeffs with Input/Output using T floating representation
166
     \input fijk: 3D real function values
167 168 169 170 171 172 173
    \output flmn: 3D complex spherical-laguerre coefficients
    \output almk: 2D alm complex spherical coefficients for each subshell k
  */
  void Analysis(const vector<r_8>& fijk, 
		vector< complex<r_8> >& flmn,
		vector< complex<r_8> >& almk
		);
174 175

protected:
176 177 178

  BaseGeometry* sphere_; //<! the 2D pixelization of the shells 

179
  LaguerreTransform* lagTrans_; //!< the Laguerre Transform
180
  
181 182
  int Npix_; //!< Total number of 2D pixels
  int Nalm_; //!< Total number of Alm coefficients
183 184 185 186
  
  int Lmax_; //!< Spherical harmonic limit such that l:0,...,Lmax-1   m:0,...,l
  int Mmax_; //!< maximum value of m

187 188 189 190

  int N_; //!< Order of the Generalized Lagurerre polynomial
  int alpha_;  //!< second parameter of the generalized Laguerre polynomial 

191 192 193
  r_8 R_; //!< Radial maximal value (not yet used)


194 195 196 197 198 199
}; //class


}//end namespace
#endif