(JEC) 24/11/15 add more detils on Fourier-Bessel part

doc/Documentation.md

@@ -3,7 +3,7 @@ LagSHT                        {#mainpage}
 
 Introduction
 ------------
-**LagSHT** is an acronym for *Laguerre Spherical Harmonic Transform*. It is also the **namespace** use throughout the C++ code to isolate the classes and the fucntions. 
+**LagSHT** is an acronym for *Laguerre Spherical Harmonic Transform*, although it includes since the version v1.7 a *Bessel Spherical Harmonic Transform* deduced from the *Laguerre Transform*. **LagSHT** is also the **namespace** use throughout the C++ code to isolate the classes and the fucntions as well as the root for the header files to be included in user applications.
 
 The original work by **Boris Leistedt** and **Jason D. McEwen** ([see the ArXiv paper](http://arxiv.org/pdf/1205.0792v2.pdf)) as well as their [FLAG software](http://www.jasonmcewen.org/codes/flag/) have been the starting point of the present library development although we differ from the FLAG 1.0b1 version by some aspects those the most important are: 
 
@@ -34,7 +34,7 @@ $$\begin{eqnarray}
 
 with $ L_n^{(\alpha)}(r)$ the generalized Laguerre polynomials (notice that by default $ \alpha = 2$). The $\tau$-scaling parameter is a free parameter but can be nicely set to $\tau = R_{max}/r_{N-1}$ with $R_{max}$ the limit in radial direction of the function $f(r,\Omega)$, and $r_{N-1}$ the largest node of the Gauss-Laguerre quadrature (see below).
 
-For bnad-limited functions the discret sums are truncated with $ l <L$ and $ n < N$ (notice the L band-limit definition is not the same used for instance in `libsharp`). Focusing on the radial part only, one gets
+For band-limited functions the discret sums are truncated with $ l <L$ and $ n < N$ (notice the L band-limit definition is not the same used for instance in `libsharp`). Focusing on the radial part only, one gets
 $$
 f(r) = \sum_{n=0}^{N-1} f_n\ {\cal K}_n(r)
 $$
@@ -74,16 +74,16 @@ We take advantage from the Matrix Multiplication writing of the $f_{lmn} \leftri
 In addition to the **Fourier-Laguerre** transform, the computation of the **Fourier-Bessel** is performed thanks to the link between the two sets of coefficients $f_{lmp}^{FB}$ and $f_{lmn}^{FL}$: 
 $$\begin{eqnarray}
 f_{lmp}^{FB} &=& \sum_n f_{lmn}^{FL} J_{ln}(k_{lp})\\
-J_{ln}(k_{lp}) &=& \sqrt{\frac{2}{\pi}}\int_0^\infty j_l(k_{lp}r){\cal K}_n(r,\tau) r^2 dr
+J_{lnp}\equiv J_{ln}(k_{lp}) &=& \sqrt{\frac{2}{\pi}}\int_0^\infty j_l(k_{lp}r){\cal K}_n(r,\tau) r^2 dr
 \end{eqnarray}$$
-where the $k_{lp}$ are related to the zeros of the $j_l$ Spherical Bessel functions and the maximal radius as $k_{lp} = q_{lp}/R_{max}$. From the set of $f_{lmp}^{FB}$ coefficients then one can synthesis the $f(r,\Omega)$ function using
+where the $k_{lp}$ are related to the zeros of the $j_l$ Spherical Bessel functions and the maximal radius as $k_{lp} = q_{lp}/R_{max}$. The $J_{lnp}$ are computed using a combnaison of a Discrete Fourier Bessel Transform (B. Leistedt et al., A&A 540, A92 (2012)) and a 39-pts Clenshaw-Curtis quadrature. As the computation is CPU intensive the set of $J_{lnp}$ are loaded from previous computation of the same set of input parameters $(L_{max}, N_{max}, P_{max})$ (ie. the $R_{max}$ parameter through the $\tau$ parameter is a scaling factor). The input/output $J_{lnp}$ are stored in the **data** directory. From the set of $f_{lmp}^{FB}$ coefficients then one can synthesis the $f(r,\Omega)$ function using
 $$
 \begin{equation}
-f(r,\Omega) = \frac{\sqrt{2 \pi}}{R_{max}^3} \sum_{lm} \left( \sum_p^\infty f^{FB}_{lmp} \ \frac{j_l(k_{lp} r)}{j_{l+1}^2(k_{lp}R_{max})} \right) Y_{lm}(\Omega)
+f(r,\Omega) = \frac{\sqrt{2 \pi}}{R_{max}^3} \sum_{lm} \left( \sum_{p=0}^{P-1} f^{FB}_{lmp} \ \frac{j_l(k_{lp} r)}{j_{l+1}^2(k_{lp}R_{max})} \right) Y_{lm}(\Omega)
 \label{eq-fr-discret-FB2}
 \end{equation}
 $$
-Notice in practice the infinite sum over the $p$ indice, i.e there is an infinity of zeros of the $j_l$ function, is truncated and this leads to numerical inaccuracy. 
+The truncation of the sum leads to numerical inaccuracy which indicates that a user tuning is necessary to get reliable results especially for the innermost shell with $P_{max}$ much higher than $N_{max}$. 
 The main classes in brief
 ------------------------
 + The LagSHT::LaguerreFuncQuad class deals with the computation of the weights and nodes of the Gauss-Laguerre quadrature with special care of the numerical validity of the computation. 
@@ -98,4 +98,4 @@ with $r_k$ the $k$-th roots of the Laguerre polynomial of order $N$.
 
 + The LagSHT::LaguerreSphericalTransform class performs the combinaison of the Laguerre transform and the Spherical Harmonic transform. 
 
-+ The LagSHT::Laguerre2Bessel class deals with the Fourier-Bessel transform based on the Fourier-Laguerre coefficients. To compute the $J_{ln}(k_{lp})$ integral I use specific quadrature and integration strategy located in the file `quadinteg.h`.
++ The LagSHT::Laguerre2Bessel class deals with the Fourier-Bessel transform based on the Fourier-Laguerre coefficients. To compute the $J_{ln}(k_{lp})$ integral I use specific quadrature and integration strategy located in the file `quadinteg.h` as well as a Discrete Fourier-Bessel Transform.
\ No newline at end of file