Skip to content

GitLab

Commit 0ceb5d88 authored by Jean-Eric Campagne's avatar Jean-Eric Campagne
Browse files

(JEC) 7/10/15 Fourier-Bessel coef from Fourier-Laguerre ones.

parent 9ea9156e
Hide whitespace changes
Inline Side-by-side
Showing with 235 additions and 110 deletions
Makefile
View file @ 0ceb5d88
......@@ -111,17 +111,18 @@ CXXSHOBJ = laguerreBuilder.o \
#C++ common Headers
CXXHDR = $(SRC)lagsht_exceptions.h \
$(SRC)lagsht_numbers.h \
$(SRC)lagsht_utils.h \
$(SRC)lagsht_geom.h \
$(SRC)lagsht_spheregeom.h \
$(SRC)lagsht_healpixhelper.h \
$(SRC)laguerreBuilder.h \
$(SRC)laguerreTransform.h \
$(SRC)lagSphericTransform.h \
$(SRC)walltimer.h \
$(SRC)lagsht_bessel.h \
$(SRC)laguerre2bessel.h
$(SRC)lagsht_numbers.h \
$(SRC)lagsht_utils.h \
$(SRC)lagsht_geom.h \
$(SRC)lagsht_spheregeom.h \
$(SRC)lagsht_healpixhelper.h \
$(SRC)laguerreBuilder.h \
$(SRC)laguerreTransform.h \
$(SRC)lagSphericTransform.h \
$(SRC)walltimer.h \
$(SRC)lagsht_bessel.h \
$(SRC)laguerre2bessel.h \
$(SRC)quadinteg.h
CPPFLAGS += $(BLASYES) $(SHARPINC) $(BLASINC)
LDFLAGS += $(SHARPLIBN) $(BLASLIBN) -lm
......
src/lagsht_testsuite.cc
View file @ 0ceb5d88
......@@ -295,8 +295,8 @@ void TestLaguerre2Bessel() {
cout << " ______________ TestLaguerre2Bessel TEST ___________ " << endl;
cout << "tstack 1 Start" <<endl;
tstack_push("data input");
cout << "tstack 1 Start" <<endl;
tstack_push("data input");
int state = 1234567 + 8912 ; //random seed
......@@ -313,7 +313,6 @@ void TestLaguerre2Bessel() {
BaseGeometry* sphere = sphlagtrans.GetSphericalGeometry();
LaguerreTransform* lagTrans = sphlagtrans.GetLagTransform();
int Nalm = sphere->Nalm();
int Nshell = sphlagtrans.GetOrder();
int Ntot = Nshell * Nalm; //total number of Coefficients of the Spherical Laguerre transform
......@@ -350,18 +349,18 @@ void TestLaguerre2Bessel() {
cout << "tstack 2 Start" <<endl;
tstack_push("processing");
int Pmax = Nmax; //To see
int Pmax = 32*Nmax; //To see
cout << "tstack 2a Start" <<endl;
tstack_push("Init Laguerre 2 Bessel");
Laguerre2Bessel lag2bess(sphere, lagTrans,
Nmax,Pmax,Rmax);
Laguerre2Bessel lag2bess(sphere,lagTrans,Nmax,Pmax,Rmax);
cout << "tstack 2a End" <<endl;
tstack_pop("Init Laguerre 2 Bessel");
cout << "tstack 2b Start" <<endl;
tstack_push("Compute Fourier-Bessel coeffs.");
vector< complex<r_8> > FBlmp(Nalm*Pmax);
lag2bess.Transform(flmnOrig,FBlmp);
lag2bess.Lag2BessCoeff(flmnOrig,FBlmp);
cout << "tstack 2b End" <<endl;
tstack_pop("Compute Fourier-Bessel coeffs.");
......@@ -369,13 +368,12 @@ void TestLaguerre2Bessel() {
cout << "tstack 2c Start" <<endl;
tstack_push("Compute Alm(rk) from FB coeffs.");
vector< complex<r_8> > FBalmk;
vector<r_8> fFBijk;
vector<r_8> fFBijk; //Not yet computed JEC 22/9/15
lag2bess.Synthesis(FBlmp,FBalmk,fFBijk);
cout << "tstack 2c End" <<endl;
tstack_pop("Compute Alm(rk) from FB coeffs.");
//A faire le calcul des Almk via LaguerreSphericalTransform::Synthesis
cout << "tstack 2d Start" <<endl;
tstack_push("Compute Alm(rk) from FB coeffs.");
vector< complex<r_8> > FLalmk;
......@@ -384,8 +382,8 @@ void TestLaguerre2Bessel() {
cout << "tstack 2d End" <<endl;
tstack_pop("Compute Alm(rk) from FB coeffs.");
cout << "tstack 2 End" <<endl;
tstack_pop("processing");
cout << "tstack 2 End" <<endl;
tstack_pop("processing");
// //FLaguerre Coeff.
// for(int n=0;n<Nmax;n++){
......@@ -400,21 +398,24 @@ void TestLaguerre2Bessel() {
// }
// }
// //FBessel Coeff.
// for(int p=0;p<Pmax;p++){
// for(int l=0;l<Lmax;l++){
// for (int m=0;m<=l;m++) {
// int id= p*Nalm + l+m*Lmax-m*(m+1)/2; // LagSHT numbering
// cout << "FB("
// <<l<<","
// <<m<<","
// <<p<<"): " << setprecision(12) << FBlmp[id] << endl;
// }
// }
// }
// // //FBessel Coeff.
// // for(int p=0;p<Pmax;p++){
// // for(int l=0;l<Lmax;l++){
// // for (int m=0;m<=l;m++) {
// // int id= p*Nalm + l+m*Lmax-m*(m+1)/2; // LagSHT numbering
// // cout << "FB("
// // <<l<<","
// // <<m<<","
// // <<p<<"): " << setprecision(12) << FBlmp[id] << endl;
// // }
// // }
// // }
cout << "Dump FL or FB reconstructed Alm(r_k)" <<endl;
std::ofstream ofs;
ofs.open ("Almi.txt", std::ofstream::out);
for(int n=0;n<Nmax;n++){
for(int l=0;l<Lmax;l++){
for (int m=0;m<=l;m++) {
......@@ -424,10 +425,16 @@ void TestLaguerre2Bessel() {
<<m<<","
<<n<< "): FL : " << setprecision(12) << FLalmk[id]
<< " FB : " << FBalmk[id] << endl;
ofs <<l<<" "
<<m<<" "
<<n<< " " << setprecision(12) << FLalmk[id].real() << " " << FLalmk[id].imag()
<< " " << FBalmk[id].real() << " " << FBalmk[id].imag() << endl;
}
}
}
ofs.close();
......@@ -597,8 +604,8 @@ int main(int narg, char *arg[]) {
{
tstack_push("TestLaguerre2Bessel");
TestLaguerre2Bessel();
tstack_pop("TestLaguerre2Bessel");
tstack_report("TestLaguerre2Bessel");
tstack_pop("TestLaguerre2Bessel");
tstack_report("TestLaguerre2Bessel");
}
break;
......
src/laguerre2bessel.cc
View file @ 0ceb5d88
......@@ -5,8 +5,15 @@
#include "walltimer.h" //timing
namespace LagSHT {
/* 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
*/
Laguerre2Bessel::Laguerre2Bessel(BaseGeometry* sphere,
LaguerreTransform* lagTrans,
int Nmax, int Pmax, r_8 Rmax
......@@ -24,60 +31,115 @@ namespace LagSHT {
throw LagSHTError("LaguerreSphericalTransform::Ctor sphere_ ptr NULL");
}
/*!compute J_{lpn} by Gauss-Laguerre quadrature.
J_{lpn} = J_{ln}(k_{lp}) = \int_0^\infty dr r^2 {\cal L}^{(2)}_n (r) j_l(k_{lp}r)
J_{lpn} approximated by Sum_{k=0}^K0 w_k {\cal L}^{(2)}(r_k) j_l(k_{lp} r_k)
Use laguerreTransform class even if it computes more coeff as needed.
*/
ComputeJInteg();
K0_ = 4*Nmax; //Heuristic Quadrature Order; should be > Nmax
LaguerreTransform lagT(K0_,Rmax);
int stride = Lmax_*Pmax; // nbre of (l,p) couples
vector<r_8> glnodes = lagT.GetNodes(); //r_k of Gauss-Laguerre Quadrature
jnu_ = new Bessel(Lmax_,Pmax); //compute qlp: l:0...,Lmax-1 et p:0,...,Pmax-1
vector<r_8> jlpk(stride*K0_); //values of j_l(k_{lp} r_k)
int lp=0;
for (int lc = 0; lc < Lmax_; lc++) {
for (int pc = 0; pc < Pmax; pc++) {
r_8 klp = (*jnu_)(lc,pc)/Rmax;
for (int k=0; k<K0_ ; k++) {
jlpk[lp+k*stride] = Bessel::jn(lc, klp * glnodes[k]);
}//end k
lp++;
}//end pc
}//end lc
//J_{lpn} approximated by Sum_{k=0}^K0 w_k {\cal L}^{(2)}(r_k) j_l(k_{lp} r_k)
// vector<r_8> Jlpn; //will be resize to stride*K0 elements while we need n:0...,Nmax-1
lagT.MultiAnalysisR(jlpk,Jlpn_,stride);
}//Ctor
/*! 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)) );
}
// cout << "Rmax, Lmax, Nmax, Pmax: "
// << Rmax_ << ", "
// << Lmax_ << ", "
// << Nmax_ << ", "
// << Pmax_ << endl;
jnu_ = new Bessel(Lmax_,Pmax_);
int Jstride = Lmax_*Pmax_; // nbre of (l,p) couples
Jlnp_.resize(Nmax_*Jstride);
vector<r_8> vTest = lagTrans_->GetNodes();
r_8 integLowerBound = 0.;
r_8 integUpperBound = 2*Rmax_/tau; //for the integral
r_8 tol = 1e-6;
tstack_push("J Integ compute");
for(int p=0; p<Pmax_; p++){
for(int l=0; l<Lmax_; l++){
int ploff7 = p + l*Pmax_;
r_8 klp = (*jnu_)(l,p)/Rmax_*tau; //tau-rescaling of the original klp definition
for (int n=0; n<Nmax_; n++){
LaguerreFuncQuad* Ln = new LaguerreFuncQuad(n);
IntFunc1D f(l,klp,Ln);
theQuad.SetFuncBounded(f,integLowerBound,integUpperBound);
Quadrature<r_8,Integrator_t>::values_t loc_val = theQuad.LocalAdapStrat(tol);
Jlnp_[ploff7 + n*Jstride] = loc_val.first*facts[n];
delete Ln;
}//n
}//l
}//p
tstack_pop("J Integ compute");
/*
std::ofstream ofs;
ofs.open ("jlnp.txt", std::ofstream::out);
Print(ofs,1);
ofs.close();
*/
}//ComputeJlpnInteg
/*!
FBlmp = FB_{lm}(k_{lp}) = Sqrt[2/Pi] Sum_n FL_{lmn} J_{lpn}
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
*/
void Laguerre2Bessel::Transform(const vector< complex<r_8> >& FLlmn,
vector< complex<r_8> >& FBlmp){
void Laguerre2Bessel::Lag2BessCoeff(const vector< complex<r_8> >& FLlmn,
vector< complex<r_8> >& FBlmp){
int Jstride = Lmax_*Pmax_;
r_8 fact = sqrt(2.0/M_PI);
//try here a pedestrian computation
for(int p=0; p<Pmax_; p++){
for(int l=0; l<Lmax_; l++){
int ploff7 = p + l*Pmax_;
for (int m=0; m<=l; m++) {
int almoff7 = l+m*Lmax_-m*(m+1)/2;
int idBess = p*Nalm_ + almoff7;
FBlmp[idBess] = 0.0;
for (int n=0; n<Nmax_; n++){
int idLag = n*Nalm_ + almoff7;
FBlmp[idBess] += FLlmn[idLag]*Jlpn_[ploff7 + n*Jstride];
FBlmp[idBess] += FLlmn[idLag]*Jlnp_[ploff7 + n*Jstride];
}//n
FBlmp[idBess] *= fact;
}//m
}//l
}//p
}//Transform
......@@ -88,6 +150,8 @@ void Laguerre2Bessel::Synthesis(const vector< complex<r_8> >& FBlmp,
int NtotCoeffFB = Nalm_*Pmax_;
if(FBlmp.size() != (size_t)NtotCoeffFB) throw LagSHTError("Laguerre2Bessel::Synthesis size error");
cout << "Dump Laguerre2Bessel::Synthesis START "<<endl;
// int Ntot3DPix = Npix_*Nmax_;
// fijk.resize(Ntot3DPix); //Not used for the moment
......@@ -96,14 +160,21 @@ void Laguerre2Bessel::Synthesis(const vector< complex<r_8> >& FBlmp,
FBalmk.resize(NtotCoeffFL);
/*
FBalmk = alm^{FB}(r_k) = Sum_{p=0}^{Pmax-1} kapa_{lp} FBlmp k_{lp} j_l(k_{lp} r_k)
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)
with
kapa_{lp} = Sqrt(2pi) Rmax^{-3}/j_{l+1}^2(q_{lp})
k_{lp} = q_{lp}/Rmax
R_k = tau r_k
r_k: k-th nodes of Laguerre of order Nmax as for the Spherical-Laguerre Transform.
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
Rdeminder: Fourier-Laguerre case where
alm^{FL}(r_k) = Sum_{n=0}^{Nmax-1} FLlmn LaguerreFunc_n(r_k)
alm^{FL}(R_k) = Sum_{n=0}^{Nmax-1} FLlmn LaguerreFunc_n(r_k)
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}.)
......@@ -111,24 +182,33 @@ void Laguerre2Bessel::Synthesis(const vector< complex<r_8> >& FBlmp,
*/
vector<r_8> vrk = lagTrans_->GetNodes();
//OK cout << "Verif: tau " << lagTrans_->GetTau() << ", Rmax/r[N-1]: " << Rmax_/vrk[Nmax_-1] << endl;
r_8 rNm1 = vrk[Nmax_-1];
//pedestrian computation
r_8 fact = sqrt(2.*M_PI)/(Rmax_*Rmax_*Rmax_);
for(int k=0;k<Nmax_;k++){
r_8 rk = vrk[k];
r_8 rk_red = vrk[k]/rNm1;
for(int l=0;l<Lmax_;l++){
for(int m=0;m<=l;m++){ //TODO use Mmax also to limit "m"
int lmOff7 = l+m*Lmax_-m*(m+1)/2;
int id0 = k*Nalm_ + lmOff7;
FBalmk[id0]=0.;
FBalmk[id0]=0.;
for(int p=0;p<Pmax_;p++){
r_8 qlp = (*jnu_)(l,p);
r_8 klp = qlp/Rmax_;
r_8 jlp1 = Bessel::jn(l+1,qlp);
r_8 kapalp = fact/(jlp1*jlp1);
int id1 = p*Nalm_ + lmOff7;
FBalmk[id0] += kapalp*klp*Bessel::jn(l,klp*rk)*FBlmp[id1];
FBalmk[id0] += kapalp*Bessel::jn(l,qlp*rk_red)*FBlmp[id1];
}//loop on p
}//loop on m
}//loop on l
......@@ -138,17 +218,16 @@ void Laguerre2Bessel::Synthesis(const vector< complex<r_8> >& FBlmp,
void Laguerre2Bessel::Print(ostream& os, int what){
//debug of the Bessel-Laguerre scalar product
if(what>0){
//Debug the Jlnp coefficients
int stride = Lmax_*Pmax_;
int lp=0;
for (int lc = 0; lc < Lmax_; lc++) {
for (int pc = 0; pc < Pmax_; pc++) {
for (int pc = 0; pc < Pmax_; pc++) {
for (int lc = 0; lc < Lmax_; lc++) {
int ploff7 = pc + lc*Pmax_;
for (int nc = 0; nc < Nmax_; nc++) {
os << lc<<" "<<pc<<" "<<nc<<" " << setprecision(12)<< Jlpn_[lp + stride*nc] << endl;
os << lc<<" "<<nc<<" "<<pc<<" " << setprecision(12)<< Jlnp_[ploff7 + stride*nc] << endl;
}//loop nc
lp++;
}//loop pc
}//loop lc
}//what
......
src/laguerre2bessel.h
View file @ 0ceb5d88
......@@ -10,6 +10,8 @@
#include "lagsht_spheregeom.h"
#include "lagsht_bessel.h"
#include "quadinteg.h"
using namespace std;
namespace LagSHT {
......@@ -19,7 +21,7 @@ namespace LagSHT {
////////////////////////////////////////////////////////////////////////////////
class Laguerre2Bessel {
public:
public:
//! Ctor
Laguerre2Bessel(BaseGeometry* sphere,
......@@ -36,12 +38,10 @@ class Laguerre2Bessel {
//! Transform Fourier-Laguerre into Fourier-Bessel
/*! \brief Fourier-Laguerre Coeffs -> Fourier-Bessel Coeffs
\input flmn: 3D complex spherical-laguerre coefficients (see LaguerreSphericalTransform)
\input Nalm: number of (l,m) couples to access the data
\output flmp: 3D complex spherical-bessel coefficients
\input FLlmn: 3D complex spherical-laguerre coefficients (see LaguerreSphericalTransform)
\output FBlmp: 3D complex spherical-bessel coefficients
*/
void Transform(const vector< complex<r_8> >& FLlmn,
vector< complex<r_8> >& FBlmp);
void Lag2BessCoeff(const vector< complex<r_8> >& FLlmn, vector< complex<r_8> >& FBlmp);
//! Synthesis
......@@ -57,8 +57,29 @@ class Laguerre2Bessel {
//!Debug
void Print(ostream& os, int what = -1);
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)
*/
void ComputeJInteg();
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
protected:
protected:
BaseGeometry* sphere_; //<! the 2D pixelization of the shells (not the owner)
LaguerreTransform* lagTrans_; //!< the Laguerre Transform
......@@ -74,10 +95,10 @@ class Laguerre2Bessel {
r_8 Rmax_; //!< Radial maximal value
int K0_; //!< Order of the LaguerreTransform to compute the Bessel-Laguerre scalar products Jlpn
//JEC 28/9/15 int K0_; //!< Order of the LaguerreTransform to compute the Bessel-Laguerre scalar products Jlpn
Bessel* jnu_; //!< Spherical Bessel fnt (owner of the ptr)
vector<r_8> Jlpn_; //!< the Bessel-Laguerre scalar product
vector<r_8> Jlnp_; //!< the Bessel-Laguerre scalar product
}; //class
......
src/laguerreBuilder.h
View file @ 0ceb5d88
......@@ -32,17 +32,17 @@ public:
fbig(std::ldexp(1.,large_exponent2)),
xi(N+2), cf(maxscale+1-minscale) {
using namespace std;
xi[0]=0.;
for (size_t i=1; i<xi.size(); ++i)
xi[i]=1./i;
for (size_t i=0; i<cf.size(); ++i)
cf[i] = ldexp(1.,(int(i)+minscale)*large_exponent2);
}
//! Destructor
virtual ~LaguerreFuncQuad() {}
virtual ~LaguerreFuncQuad() {
}
//! Getters
int GetOrder() const { return N_; }
......
src/laguerreTransform.cc
View file @ 0ceb5d88
......@@ -21,7 +21,7 @@ using namespace std;
namespace LagSHT {
LaguerreTransform::LaguerreTransform(int N, r_8 R, int alpha) :
LaguerreTransform::LaguerreTransform(int N, r_8 R, int alpha, bool R2tau) :
BaseLaguerreTransform(N,alpha), R_(R) {
#if DEBUG >= 1
......@@ -32,18 +32,32 @@ LaguerreTransform::LaguerreTransform(int N, r_8 R, int alpha) :
LaguerreFuncQuad lag(N_,alpha_);
lag.QuadWeightNodes(nodes_,weights_);
int alphaFact = 1;
for(int i=1;i<=alpha_; i++) alphaFact *= i;
//JEC 30/9/15 using alphaFact & alphaFact_ with different types is missleading
int tmp = 1;
for(int i=1;i<=alpha_; i++) tmp *= i;
alphaFact_ = sqrt((r_8)alphaFact); //sqrt(alpha!)
alphaFact_ = sqrt((r_8)tmp); //sqrt(alpha!)
//JEC 22/9/15
if(R2tau) {
tau_ = R_/nodes_[N_-1];
r_8 pow3inv2 = 3./2.;
tau3sqrt_ = pow(tau_,pow3inv2);
tau3sqrtinv_ = 1./tau3sqrt_;
} else { //this case may be used when Rmax is not bounded
tau_ = 1.;
tau3sqrt_ = 1.;
tau3sqrtinv_ = 1.;
R_ = nodes_[N_-1];
}
}//Ctor
int LaguerreTransform::R2Index(r_8 r) const {
using namespace std;
r_8 rscaled = r*nodes_[N_-1]/R_;
//JEC 22/9/15 r_8 rscaled = r*nodes_[N_-1]/R_;
r_8 rscaled = r/tau_;
vector<r_8> dist(N_);
transform(nodes_.begin(),nodes_.end(),dist.begin(), distance(rscaled));
return min_element(dist.begin(),dist.end()) - dist.begin();
......@@ -62,7 +76,7 @@ void LaguerreTransform::MultiAnalysis(const vector< complex<r_8> >& fi,
#if DEBUG >= 2
std::copy(fi.begin(), fi.end(), std::ostream_iterator<char>(std::cout, " "));
#endif
r_8 invalphaFact = 1.0/alphaFact_; //1/alpha!
r_8 invalphaFact = 1.0/alphaFact_; //1/sqrt(alpha)!
vector<r_8> facts(N_); //sqrt(n!/(n+alpha)!) en iteratif
facts[0] = invalphaFact;
......@@ -76,8 +90,7 @@ void LaguerreTransform::MultiAnalysis(const vector< complex<r_8> >& fi,
r_8 rk = nodes_[k];
lag.Values(rk,LnAll);
for (int n = 0; n<N_; n++ ){
// LnkMtx[n*N_+k] = LnAll[n]*facts[n];
LnkMtx[n+N_*k] = LnAll[n]*facts[n]*weights_[k];
LnkMtx[n+N_*k] = LnAll[n]*facts[n]*weights_[k]*tau3sqrt_; //JEC 22/9/15 add tau-fact
}
}
......@@ -132,8 +145,7 @@ void LaguerreTransform::MultiAnalysisR(const vector<r_8>& fi,
r_8 rk = nodes_[k];
lag.Values(rk,LnAll);
for (int n = 0; n<N_; n++ ){
// LnkMtx[n*N_+k] = LnAll[n]*facts[n];
LnkMtx[n+N_*k] = LnAll[n]*facts[n]*weights_[k];
LnkMtx[n+N_*k] = LnAll[n]*facts[n]*weights_[k]*tau3sqrt_; //JEC 22/9/15 add tau-fact
}
}
......@@ -191,7 +203,7 @@ void LaguerreTransform::MultiSynthesis(const vector< complex<r_8> >& fn,
r_8 rk = nodes_[k];
lag.Values(rk,LnAll);
for (int n = 0; n<N_; n++ ){
LnkMtx[n*N_+k] = LnAll[n]*facts[n];
LnkMtx[n*N_+k] = LnAll[n]*facts[n]*tau3sqrtinv_; //JEC 22/9/15 add tau-fact
}
}
......
src/laguerreTransform.h
View file @ 0ceb5d88
......@@ -33,20 +33,22 @@ protected:
class LaguerreTransform : public BaseLaguerreTransform {
public:
//! Creator
LaguerreTransform(int N, r_8 R, int alpha=2);
LaguerreTransform(int N, r_8 R, int alpha=2, bool R2tau=true);
//! Destructor
virtual ~LaguerreTransform() {}
//! Getters
vector<r_8> GetNodes() const { return nodes_; }
vector<r_8> GetWeights() const { return weights_; }
r_8 GetTau() const { return tau_; }
//! Find the closest nodes: no boundchecking
int R2Index(r_8 r) const;
//! return the nodes position: no boundchecking
r_8 Index2R(int k) const {
return nodes_[k]*R_/nodes_[N_-1];
// return nodes_[k]*R_/nodes_[N_-1];