laguerre2bessel.cc 6.63 KB
Newer Older
1 2 3
#include "laguerre2bessel.h"
#include "lagsht_bessel.h"
#include "laguerreTransform.h"
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
4 5 6
#include "lagsht_exceptions.h"


7

8 9
#include "walltimer.h"    //timing

10 11
#include <omp.h> //TEST 

12 13
namespace LagSHT {

14 15 16 17 18
  /* Version 0: test de calcul des J_{lnp} via une quadrature de Laguerre 
     mais ca n'a pas l'air de converger (28 Sept 15)
     Version 1: on fait une integration a la Clenshaw-Curtis avec startegie Localeo ou Globale
     recursive
  */
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
19 20 21 22 23 24
  Laguerre2Bessel::Laguerre2Bessel(BaseGeometry* sphere,
				   LaguerreTransform* lagTrans,
				   int Nmax, int Pmax, r_8 Rmax
				   ):
    sphere_(sphere),lagTrans_(lagTrans),
    Nmax_(Nmax), Pmax_(Pmax), Rmax_(Rmax) {
25

Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
26 27 28 29 30 31 32 33 34 35
    if(sphere_) {
      //store localy usefull parameters for the algorithms Analysis/Synthesis
      Npix_ = sphere_->Npix();
      Nalm_ = sphere_->Nalm();
      Lmax_ = sphere_->Lmax();
      Mmax_ = sphere_->Mmax();
    } else {
      throw LagSHTError("LaguerreSphericalTransform::Ctor sphere_ ptr NULL");
    }
 
36
    ComputeJInteg();
37 38 39

}//Ctor

40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
  /*! compute J_{lnp} = J_{ln}(k_{lp} via Clenshaw-Curtis quadrature and local integration strategy
   */
  void Laguerre2Bessel::ComputeJInteg(){

    typedef  ClenshawCurtisQuadrature<double> Integrator_t;
    Integrator_t theQuad(20,"./data/ClenshawCurtisRuleData-20.txt",false); //false=do not recompute the weights and nodes of the quadrature

    r_8 tau = lagTrans_->GetTau();
    r_8 tau3sqrt = pow(tau,(r_8)(3./2.));
    //    cout << "ComputeJlpnInteg: tau " << tau << ", tau^(3/2): " << tau3sqrt << endl;


    //normalisation factor 
    int alpha = lagTrans_->GetAlpha();
    r_8 alphaFact = 1; // alpha!
    for(int i=1;i<=alpha; i++) alphaFact *= i;
    alphaFact = sqrt((r_8)alphaFact); //sqrt(alpha!)
    r_8 invalphaFact = 1.0/alphaFact; // 1/sqrt(alpha)!                   

    vector<r_8> facts(Nmax_); //sqrt(2/pi)*tau^(3/2)*sqrt(n!/(n+alpha)!) en iteratif
    facts[0] = invalphaFact*tau3sqrt*sqrt(2.0/M_PI); //JEC 7/10/15 add sqrt(2/pi)
    for(int n=1;n<Nmax_;n++) {
      facts[n] = facts[n-1]*sqrt(((r_8)n)/((r_8)(n+alpha)) ); 
    }


    jnu_ = new Bessel(Lmax_,Pmax_);
    int Jstride = Lmax_*Pmax_; // nbre of (l,p) couples
    Jlnp_.resize(Nmax_*Jstride);

    
    r_8 integLowerBound = 0.;
72
    r_8 integUpperBound = 10*Rmax_/tau; //for the integral 
73 74 75 76
    r_8 tol = 1e-6;

    tstack_push("J Integ compute");

77
    Quadrature<r_8,Integrator_t>::values_t loc_val;
78 79
    for(int p=0; p<Pmax_; p++){
      for(int l=0; l<Lmax_; l++){
80

81 82
	int ploff7 = p + l*Pmax_;
 	r_8 klp = (*jnu_)(l,p)/Rmax_*tau; //tau-rescaling of the original klp definition
83

84
	for (int n=0; n<Nmax_; n++){
85

86
	  LaguerreFuncQuad* Ln = new LaguerreFuncQuad(n);
87

88 89
	  IntFunc1D f(l,klp,Ln);
	  theQuad.SetFuncBounded(f,integLowerBound,integUpperBound);
90
	  loc_val = theQuad.LocalAdapStrat(tol);
91 92 93 94 95 96 97
	  Jlnp_[ploff7 + n*Jstride] = loc_val.first*facts[n];

	  delete Ln;
	}//n
      }//l
    }//p

98
    tstack_pop("J Integ compute");
99 100 101

  }//ComputeJlpnInteg

102
/*!
103 104 105 106 107
  FBlmp = FB_{lm}(k_{lp}) = Sqrt[2/Pi] Sum_n FL_{lmn} J_{ln}(k_{lp})
        = Sqrt[2/Pi] Sum_n FL_{lmn} J_{lnp}
  This is the Fourrier-Bessel coefficients related to McEwen & Leistedh definition
  In case one wants to use Lanusse & Starck definition then
  FBlmp_Lanusse = k_{lp} * FBlmp_McEwen
108
*/
109 110 111
void Laguerre2Bessel::Lag2BessCoeff(const vector< complex<r_8> >& FLlmn,
				    vector< complex<r_8> >& FBlmp){

112 113
  
  int Jstride = Lmax_*Pmax_;
114

115 116 117
  //try here a pedestrian computation
  for(int p=0; p<Pmax_; p++){
    for(int l=0; l<Lmax_; l++){
118

119
      int ploff7 = p + l*Pmax_;
120

121 122
      for (int m=0; m<=l; m++) {
	int almoff7 = l+m*Lmax_-m*(m+1)/2;
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
123
	int idBess = p*Nalm_ + almoff7;
124
	FBlmp[idBess] = 0.0;
125

126
	for (int n=0; n<Nmax_; n++){
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
127
	  int idLag =  n*Nalm_ + almoff7;
128
	  FBlmp[idBess] += FLlmn[idLag]*Jlnp_[ploff7 + n*Jstride];	  
129
	}//n
130

131 132 133
      }//m
    }//l
  }//p
134

135 136 137
}//Transform


Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
138 139
void Laguerre2Bessel::Synthesis(const vector< complex<r_8> >& FBlmp, 
				vector< complex<r_8> >& FBalmk,
140
				vector<r_8>& fijk
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
141 142 143 144
				) {

  int NtotCoeffFB = Nalm_*Pmax_; 
  if(FBlmp.size() != (size_t)NtotCoeffFB) throw LagSHTError("Laguerre2Bessel::Synthesis size error");
145 146

  cout << "Dump  Laguerre2Bessel::Synthesis START "<<endl;
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
147 148 149 150 151 152
    

  int NtotCoeffFL = Nalm_*Nmax_; 
  FBalmk.resize(NtotCoeffFL); 
  
  /*
153 154
    This is the McEwen & Leistedh definition of the Fourier-Bessel Transform
    FBalmk = alm^{FB}(R_k) = Sum_{p=0}^{Pmax-1} Kapa_{lp} FBlmp  j_l(k_{lp} R_k)
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
155 156 157
    with
    kapa_{lp} = Sqrt(2pi) Rmax^{-3}/j_{l+1}^2(q_{lp})
    k_{lp} = q_{lp}/Rmax
158
    R_k = tau r_k
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
159
    r_k: k-th nodes of Laguerre of order Nmax as for the Spherical-Laguerre Transform.
160 161 162 163 164
    tau: Rmax/r_{N-1} from Laguerre-Transform

    23/9/15 in fact Leistedt & McEwen  et al in IEEE Vol 60  Num 12 (Dec 2012) do not 
    use the same convention on Bessel Transform as Lanusse et al. A&A 578, A10 (2015) so
    in the Sum_{p=0}^{Pmax-1}... the k_{lp} should be squared
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
165 166

    Rdeminder: Fourier-Laguerre case where
167
    alm^{FL}(R_k) = Sum_{n=0}^{Nmax-1} FLlmn LaguerreFunc_n(r_k)
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
168 169 170 171 172 173 174
    
    in the Fourier-Bessel case the element of the sum, except the FBlmp, are depending
    on "l" also throw kapa_{lp}, k_{lp} and j_l(k_{lp}.) 
   
   */
  
  vector<r_8> vrk = lagTrans_->GetNodes();
175 176 177
  //OK  cout << "Verif: tau " << lagTrans_->GetTau() << ", Rmax/r[N-1]: " << Rmax_/vrk[Nmax_-1] << endl;
  r_8 rNm1 = vrk[Nmax_-1];

Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
178 179 180 181
  //pedestrian computation
  r_8 fact = sqrt(2.*M_PI)/(Rmax_*Rmax_*Rmax_);

  for(int k=0;k<Nmax_;k++){
182 183 184

    r_8 rk_red = vrk[k]/rNm1;

Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
185 186
    for(int l=0;l<Lmax_;l++){
      for(int m=0;m<=l;m++){ //TODO use Mmax also to limit "m"
187

Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
188 189
	int lmOff7 = l+m*Lmax_-m*(m+1)/2;
	int id0 = k*Nalm_ + lmOff7;
190 191 192

 	FBalmk[id0]=0.;

Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
193
	for(int p=0;p<Pmax_;p++){
194

Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
195 196 197 198
	  r_8 qlp = (*jnu_)(l,p);
	  r_8 jlp1 = Bessel::jn(l+1,qlp);
	  r_8 kapalp = fact/(jlp1*jlp1);
	  int id1 = p*Nalm_ + lmOff7;
199 200 201

	  FBalmk[id0] += kapalp*Bessel::jn(l,qlp*rk_red)*FBlmp[id1];

Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
202 203 204 205 206
	}//loop on p
      }//loop on m
    }//loop on l
  }//loop on k

207 208 209 210 211 212 213 214 215 216 217 218 219 220 221

  int Ntot3DPix = Npix_*Nmax_;
  fijk.resize(Ntot3DPix); //Not used for the moment

  int off7n, off7k;
  //perform the synthesis per shell
  for(int n =0; n<Nmax_; n++){
    off7n = n*Nalm_;
    off7k = n*Npix_; 

    sphere_->alm2map(reinterpret_cast<const r_8*>(&(FBalmk.data()[off7n])),&fijk.data()[off7k],false); //the false is to forbidd accumulation

  }//loop on shell
  

Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
222 223 224
}//Synthesis


225 226 227
void Laguerre2Bessel::Print(ostream& os, int what){
  if(what>0){

228
    //Debug the Jlnp coefficients
229
    int stride = Lmax_*Pmax_;
230 231 232
    for (int pc = 0; pc < Pmax_; pc++) {
      for (int lc = 0; lc < Lmax_; lc++) {
	int ploff7 = pc + lc*Pmax_;
233
	for (int nc = 0; nc < Nmax_; nc++) {
234
	  os << lc<<" "<<nc<<" "<<pc<<" " << setprecision(12)<< Jlnp_[ploff7 + stride*nc] << endl;
235 236 237 238 239 240 241
	}//loop nc
      }//loop pc
    }//loop lc
  }//what
}//Print

}//namespace