Skip to content

GitLab

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

Start Bessel introduction

parent 6735ca6f
Hide whitespace changes
Inline Side-by-side
Showing with 312 additions and 28 deletions
src/lagSphericTransform.cc
View file @ d620dee7
......@@ -15,8 +15,11 @@
namespace LagSHT {
// void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& flmn,
// vector<r_8>& fijk) {
void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& flmn,
vector<r_8>& fijk) {
vector<r_8>& fijk,
vector< complex<r_8> >& flmk) {
int Ntot3DPix = Npix_*N_;
int NtotCoeff = Nalm_*N_;
if(flmn.size() != (size_t)NtotCoeff) throw LagSHTError("SphLagTransform::Synthesis size error");
......@@ -24,7 +27,9 @@ void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& flmn,
fijk.resize(Ntot3DPix);
//Perform the Inverse Laguerre Transform first
vector< complex<r_8> > flmk(NtotCoeff);
// vector< complex<r_8> > flmk(NtotCoeff);
flmk.resize(NtotCoeff);
#if DEBUG >= 1
cout << "SphLagTransform::Synthesis type <r_8> I/O .... Go..." << endl;
......
src/lagSphericTransform.h
View file @ d620dee7
......@@ -117,8 +117,22 @@ class LaguerreSphericalTransform {
/*! \brief Coeffs -> Pixels with Input/Output using T floating representation
\input flmn: 3D complex spherical-laguerre coefficients
\output fijk: 3D real function values
\output almk: 2D alm complex spherical coefficients for each subshell k
*/
void Synthesis(const vector< complex<r_8> >& flmn, vector<r_8>& fijk);
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);
}
//! Analysis
/*! \brief Pixels -> Coeffs with Input/Output using T floating representation
......
src/lagsht_spheregeom.h
View file @ d620dee7
......@@ -51,6 +51,10 @@ class BaseGeometry {
int Npix() const { return Npix_; }
//! Total number of Alm coefficients
int Nalm() const { return Nalm_; }
//! Lmax
int Lmax() const { return L_; }
//! Mmax
int Mmax() const { return M_; }
//! Type of geomery
string Geometry() const { return geomType_; }
......
src/lagsht_testsuite.cc
View file @ d620dee7
......@@ -15,6 +15,8 @@ using namespace std;
#include "lagSphericTransform.h"
#include "laguerre2bessel.h"
using namespace LagSHT;
#define DEBUG 1
......@@ -190,15 +192,6 @@ void MultiSphericalLaguerreTransform() {
nphi
);
// LaguerreSphericalTransform sphlagtrans(Lmax,-1,N,R);
// sphlagtrans.SetThetaPhiMap(geometry,ntheta,nphi);
// int Nalm = sphlagtrans.Nalm();
// int Ntot = N * Nalm; //total number of 3D-pixels & number of Coefficients of the Spherical Laguerre transform
// int Npix = sphlagtrans.Npix();
// int NpTot= N * Npix;
int Nalm = sphlagtrans.GetSphericalGeometry()->Nalm();
int Nshell = sphlagtrans.GetOrder();
......@@ -288,6 +281,158 @@ void MultiSphericalLaguerreTransform() {
}//MultiSphericalLaguerreTransform
//-------------------------------------------------------------
void TestLaguerre2Bessel() {
//Todo introduce Pmax !
int Nmax = param.N;
int Lmax = param.Lmax;
r_8 Rmax = param.R;
string geometry = param.geometry;
int ntheta = param.ntheta;
int nphi = param.nphi;
cout << " ______________ TestLaguerre2Bessel TEST ___________ " << endl;
cout << "tstack 1 Start" <<endl;
tstack_push("data input");
int state = 1234567 + 8912 ; //random seed
LaguerreSphericalTransform sphlagtrans(geometry,
Lmax,
-1,
Nmax,
Rmax,
ntheta,
nphi
);
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
int Npix = sphere->Npix();
int NpTot= Nshell * Npix; //totoal number of 3D-pixels
cout << "Verif: Npix, Nptot, Nalm, Nshell, Ntot: "
<< Npix << " "
<< NpTot << " "
<< Nalm << " "
<< Nshell << " "
<< Ntot << endl;
vector< complex<r_8> > flmnOrig(Ntot);
int id=-1;
for(int n=0;n<Nmax;n++){ //on each shell
for(int m=0;m<Lmax;m++){ //Warning the filling of alm is adapted for libsharp memory
for(int l=m;l<Lmax;l++){
id++;
r_8 rv = drand(-1,1,&state);
r_8 iv = (m==0) ? 0.0 : drand(-1,1,&state);
flmnOrig[id] = complex<r_8>(rv,iv);
}//end l-loop
}//end m-loop
}//end loop on shell
#if DEBUG >=2
std::copy(flmnOrig.begin(), flmnOrig.end(), std::ostream_iterator<char>(std::cout, " "));
#endif
cout << "tstack 1 End" <<endl;
tstack_pop("data input");
cout << "tstack 2 Start" <<endl;
tstack_push("processing");
int Pmax = Nmax; //To see
cout << "tstack 2a Start" <<endl;
tstack_push("Init Laguerre 2 Bessel");
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);
cout << "tstack 2b End" <<endl;
tstack_pop("Compute Fourier-Bessel coeffs.");
cout << "tstack 2c Start" <<endl;
tstack_push("Compute Alm(rk) from FB coeffs.");
vector< complex<r_8> > FBalmk;
vector<r_8> fFBijk;
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;
vector<r_8> fFLijk;
sphlagtrans.Synthesis(flmnOrig,fFLijk,FLalmk);
cout << "tstack 2d End" <<endl;
tstack_pop("Compute Alm(rk) from FB coeffs.");
cout << "tstack 2 End" <<endl;
tstack_pop("processing");
// //FLaguerre 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
// cout << "FL("
// <<l<<","
// <<m<<","
// <<n<<"): " << setprecision(12) << flmnOrig[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;
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
cout << "Almk ("
<<l<<","
<<m<<","
<<n<< "): FL : " << setprecision(12) << FLalmk[id]
<< " FB : " << FBalmk[id] << endl;
}
}
}
}//TestLaguerre2Bessel
//-------------------------------------------------------------
void TestPixelization() {
......@@ -397,6 +542,10 @@ int main(int narg, char *arg[]) {
nphi = atoi(arg[ka+1]);
ka+=2;
}
else if (strcmp(arg[ka],"-r")==0) {
R = atof(arg[ka+1]);
ka+=2;
}
else ka++;
}//eo while
......@@ -411,8 +560,6 @@ int main(int narg, char *arg[]) {
int rc=0;
try {
switch(test) {
case 0:
TestBasic();
......@@ -446,6 +593,15 @@ int main(int narg, char *arg[]) {
}
break;
case 5:
{
tstack_push("TestLaguerre2Bessel");
TestLaguerre2Bessel();
tstack_pop("TestLaguerre2Bessel");
tstack_report("TestLaguerre2Bessel");
}
break;
default:
throw LagSHTError("EROOR Test type unkwown");
}//end of switch
......
src/laguerre2bessel.cc
View file @ d620dee7
#include "laguerre2bessel.h"
#include "lagsht_bessel.h"
#include "laguerreTransform.h"
#include "lagsht_exceptions.h"
namespace LagSHT {
Laguerre2Bessel::Laguerre2Bessel(int Lmax, int Mmax, int Nmax, int Pmax, r_8 Rmax):
Lmax_(Lmax), Mmax_(Mmax), Nmax_(Nmax), Pmax_(Pmax), Rmax_(Rmax) {
Laguerre2Bessel::Laguerre2Bessel(BaseGeometry* sphere,
LaguerreTransform* lagTrans,
int Nmax, int Pmax, r_8 Rmax
):
sphere_(sphere),lagTrans_(lagTrans),
Nmax_(Nmax), Pmax_(Pmax), Rmax_(Rmax) {
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");
}
/*!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)
......@@ -15,17 +32,17 @@ Laguerre2Bessel::Laguerre2Bessel(int Lmax, int Mmax, int Nmax, int Pmax, r_8 Rma
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> nodes = lagT.GetNodes(); //r_k of Gauss-Laguerre Quadrature
Bessel jnu(Lmax,Pmax); //compute qlp: l:0...,Lmax-1 et p:0,...,Pmax-1
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 lc = 0; lc < Lmax_; lc++) {
for (int pc = 0; pc < Pmax; pc++) {
r_8 klp = jnu(lc,pc)/Rmax;
r_8 klp = (*jnu_)(lc,pc)/Rmax;
for (int k=0; k<K0_ ; k++) {
jlpk[lp+k*stride] = Bessel::jn(lc, klp * nodes[k]);
jlpk[lp+k*stride] = Bessel::jn(lc, klp * glnodes[k]);
}//end k
lp++;
}//end pc
......@@ -40,7 +57,6 @@ Laguerre2Bessel::Laguerre2Bessel(int Lmax, int Mmax, int Nmax, int Pmax, r_8 Rma
FBlmp = FB_{lm}(k_{lp}) = Sqrt[2/Pi] Sum_n FL_{lmn} J_{lpn}
*/
void Laguerre2Bessel::Transform(const vector< complex<r_8> >& FLlmn,
int Nalm,
vector< complex<r_8> >& FBlmp){
int Jstride = Lmax_*Pmax_;
......@@ -51,10 +67,10 @@ void Laguerre2Bessel::Transform(const vector< complex<r_8> >& FLlmn,
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;
int idBess = p*Nalm_ + almoff7;
FBlmp[idBess] = 0.0;
for (int n=0; n<Nmax_; n++){
int idLag = n*Nalm + almoff7;
int idLag = n*Nalm_ + almoff7;
FBlmp[idBess] += FLlmn[idLag]*Jlpn_[ploff7 + n*Jstride];
}//n
FBlmp[idBess] *= fact;
......@@ -65,6 +81,62 @@ void Laguerre2Bessel::Transform(const vector< complex<r_8> >& FLlmn,
}//Transform
void Laguerre2Bessel::Synthesis(const vector< complex<r_8> >& FBlmp,
vector< complex<r_8> >& FBalmk,
vector<r_8>& fijk //not computed for the moment
) {
int NtotCoeffFB = Nalm_*Pmax_;
if(FBlmp.size() != (size_t)NtotCoeffFB) throw LagSHTError("Laguerre2Bessel::Synthesis size error");
// int Ntot3DPix = Npix_*Nmax_;
// fijk.resize(Ntot3DPix); //Not used for the moment
int NtotCoeffFL = Nalm_*Nmax_;
FBalmk.resize(NtotCoeffFL);
/*
FBalmk = alm^{FB}(r_k) = Sum_{p=0}^{Pmax-1} kapa_{lp} FBlmp k_{lp} 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: k-th nodes of Laguerre of order Nmax as for the Spherical-Laguerre Transform.
Rdeminder: Fourier-Laguerre case where
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}.)
*/
vector<r_8> vrk = lagTrans_->GetNodes();
//pedestrian computation
r_8 fact = sqrt(2.*M_PI)/(Rmax_*Rmax_*Rmax_);
for(int k=0;k<Nmax_;k++){
r_8 rk = vrk[k];
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.;
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];
}//loop on p
}//loop on m
}//loop on l
}//loop on k
}//Synthesis
void Laguerre2Bessel::Print(ostream& os, int what){
//debug of the Bessel-Laguerre scalar product
if(what>0){
......
src/laguerre2bessel.h
View file @ d620dee7
......@@ -6,6 +6,9 @@
#include <complex>
#include "lagsht_numbers.h"
#include "laguerreTransform.h"
#include "lagsht_spheregeom.h"
#include "lagsht_bessel.h"
using namespace std;
......@@ -18,8 +21,18 @@ namespace LagSHT {
class Laguerre2Bessel {
public:
//!Ctor: TODO see if Mmax is finally used or not
Laguerre2Bessel(int Lmax, int Mmax, int Nmax, int Pmax, r_8 Rmax);
//! Ctor
Laguerre2Bessel(BaseGeometry* sphere,
LaguerreTransform* lagTrans,
int Nmax, int Pmax, r_8 Rmax);
//! Destructor
virtual ~Laguerre2Bessel() {
if(jnu_) delete jnu_; jnu_ = 0;
}
//! Geometry helper
BaseGeometry* GetSphericalGeometry() const { return sphere_;}
//! Transform Fourier-Laguerre into Fourier-Bessel
/*! \brief Fourier-Laguerre Coeffs -> Fourier-Bessel Coeffs
......@@ -28,14 +41,31 @@ class Laguerre2Bessel {
\output flmp: 3D complex spherical-bessel coefficients
*/
void Transform(const vector< complex<r_8> >& FLlmn,
int Nalm,
vector< complex<r_8> >& FBlmp);
//! 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);
//!Debug
void Print(ostream& os, int what = -1);
protected:
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
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)
......@@ -45,7 +75,10 @@ class Laguerre2Bessel {
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
}; //class
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment