laguerre2bessel.h 4.94 KB
Newer Older
1 2
#ifndef LAGUERRE2BESSEL_SEEN
#define LAGUERRE2BESSEL_SEEN
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
/*
 *  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
 */

/*
 *  LagSHT is being developed at the Linear Accelerateur Laboratory (LAL)
 *  91898 ORSAY CEDEX - FRANCE
 *  main author: J.E Campagne
 */
26 27 28 29

#include <iostream>
#include <vector>
#include <complex>
30
#include <utility> //move
31 32

#include "lagsht_numbers.h"
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
33 34 35
#include "laguerreTransform.h"
#include "lagsht_spheregeom.h"
#include "lagsht_bessel.h"
36
#include "lagsht_exceptions.h"
37

38 39
#include "quadinteg.h" //JEC 12/11/15

40 41 42 43 44 45 46 47
using namespace std;

namespace LagSHT {

////////////////////////////////////////////////////////////////////////////////
//// --------- Class to transform Spherical Laguerre into Spherical Bessel -- //
////////////////////////////////////////////////////////////////////////////////

48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68
/* //implementation of the Kahan summation algorithm */
/* template <class InIt> */
/* typename std::iterator_traits<InIt>::value_type accumulate(InIt begin, InIt end) { */
/*     typedef typename std::iterator_traits<InIt>::value_type real; */
/*     real sum = real(); */
/*     real running_error = real(); */
/*     real temp; */
/*     real difference; */

/*     for (; begin != end; ++begin) { */
/*         difference = *begin; */
/*         difference -= running_error; */
/*         temp = sum; */
/*         temp += difference; */
/*         running_error = temp; */
/*         running_error -= sum; */
/*         running_error -= difference; */
/*         sum = std::move(temp); */
/*     } */
/*     return sum; */
/* } */
69

70
class Laguerre2Bessel {
71
public:
72

Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
73 74 75
  //! Ctor
  Laguerre2Bessel(BaseGeometry* sphere, 
		  LaguerreTransform* lagTrans,
76 77 78 79 80
		  int Nmax, int Pmax, r_8 Rmax,
		  string clenshawDir,
		  string jlnpDir,
		  bool recomputeJlnp=false
		  );
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
81 82 83 84 85 86 87 88

  //! Destructor
  virtual ~Laguerre2Bessel() {
    if(jnu_) delete jnu_; jnu_ = 0;
  }

  //! Geometry helper
  BaseGeometry* GetSphericalGeometry() const { return sphere_;}
89 90 91

  //! Transform Fourier-Laguerre into Fourier-Bessel
  /*! \brief Fourier-Laguerre Coeffs -> Fourier-Bessel Coeffs
92 93
    \input FLlmn: 3D complex spherical-laguerre coefficients (see LaguerreSphericalTransform)
    \output FBlmp: 3D complex spherical-bessel coefficients
94
   */
95
  void Lag2BessCoeff(const vector< complex<r_8> >& FLlmn, vector< complex<r_8> >& FBlmp);
96

Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
97 98 99 100 101 102 103 104 105 106 107
  
  //! Synthesis
  /*! \brief Fourier-Bessel Coeffs -> Alm(r_k) -> pixel fijk
    \input FBlmp: 3D Fourier Bessel Coefficients
    \output FBalmk: 2D alm(r_k) Harm. Spherical reconstruction of each r_k given by Laguerre radius
    \output fijk: 3D real values of each pixel 
   */
  void Synthesis(const vector< complex<r_8> >& FBlmp, 
		 vector< complex<r_8> >& FBalmk,
		 vector<r_8>& fijk);

108 109 110
  //!Debug
  void Print(ostream& os, int what = -1);

111 112 113 114 115
private:
  
  /*! \brief
      J_{lnp} = J_{ln}(k_{lp}) 
              = tau^(3/2) \int_0^\infty dx x^2 {\cal L}^{(2)}_n (x) j_l(tau k_{lp} x)
116 117
      \input  clenshawDir directory to find the Clenshaw-Curtis quadrature weight & nodes
      \input jlnpDir Directory where to save the result
118
  */
119
  void ComputeJInteg(string clenshawDir, string jlnpDir);
120

121
 
122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
  /*!
    function for numerical integration
   */
  class IntFunc1D : public ClassFunc1D {
  public:
    IntFunc1D(int l, r_8 klp,LaguerreFuncQuad* Ln):
      l_(l),klp_(klp),Ln_(Ln) {}
    virtual double operator()(double x) const {
      return (x*x)*(Bessel::jn(l_,x*klp_))*Ln_->Value(x);
    }
    virtual ~IntFunc1D() {}
  private:
    int l_;
    r_8 klp_;
    LaguerreFuncQuad* Ln_; //no ownership
  };//class IntFunc1D



141
protected:
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
142 143 144 145 146 147 148
  BaseGeometry* sphere_; //<! the 2D pixelization of the shells (not the owner)
  
  LaguerreTransform* lagTrans_; //!< the Laguerre Transform
  
  int Npix_; //!< Total number of 2D pixels
  int Nalm_; //!< Total number of Alm coefficients

149 150 151 152 153 154 155
  int Lmax_; //!< Spherical harmonic limit such that l:0,...,Lmax-1   m:0,...,l
  int Mmax_; //!< maximum value of m (NOT USED)
  int Nmax_; //!< Max. Order of the Generalized Lagurerre polynomial  (alpha=2 by default)

  int Pmax_; //!< Max. order of the Bessel coeff.
  r_8 Rmax_; //!< Radial maximal value  

156
  //JEC 1/2/16  int PPmax_; //!< cut off in the alm computation
157

158
  
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
159
  Bessel* jnu_; //!< Spherical Bessel fnt (owner of the ptr)
160
  vector<r_8> Jlnp_; //!< the  Bessel-Laguerre scalar product
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
161
  
162 163 164 165 166 167

}; //class


}//end namespace
#endif