(JEC) 7/10/15 complete the documentation for tau-factor and Fourier-Bessel transform.

doc/Documentation.md

@@ -15,28 +15,28 @@ Mathematics in summary
 Any function $f(r,\Omega)$ square-integrable on $\mathrm{B}^3 = \mathrm{R}^+ \times [0,\pi] \times [0,2\pi)$ can be decomposed as:
 
 $$\begin{eqnarray}
-f(r,\Omega) &=& \sum_{n=0}^\infty \sum_{l=0}^{\infty}\sum_{m=-l}^{l}\ f_{lmn}\ K_{lmn}(r,\Omega) \label{eq:FLagfulla}\\
-f_{lmn} &=& \int_{\mathrm{B}^3} dr d\Omega\ r^2\ f(r,\Omega)\ K^\ast_{lmn}(r,\Omega) 
+f(r,\Omega) &=& \sum_{n=0}^\infty \sum_{l=0}^{\infty}\sum_{m=-l}^{l}\ f_{lmn}\ K_{lmn}(r,\Omega; \tau) \label{eq:FLagfulla}\\
+f_{lmn} &=& \int_{\mathrm{B}^3} dr d\Omega\ r^2\ f(r,\Omega)\ K^\ast_{lmn}(r,\Omega; \tau) 
 	\label{eq:FLagfullb}
 \end{eqnarray}$$
 
 based on the orthogonality of the set of functions 
 $$\begin{equation}
-K_{lmn}(r,\Omega) \equiv  Y_{l,m}(\Omega)\times {\cal L}_n(r)
+K_{lmn}(r,\Omega; \tau) \equiv  Y_{l,m}(\Omega)\times {\cal K}_n(r,\tau)
 \label{eq:FLagOrthoFunc}
 \end{equation}$$
-using the spherical harmonic functions $Y_{l,m}$ and the genralized Laguerre functions ${\cal L}_n$ 
+using the spherical harmonic functions $Y_{l,m}$ and the genralized Laguerre functions ${\cal K}_n$: 
 
 $$\begin{eqnarray}
  Y_{l,m}(\Omega) &=& \lambda_{lm}(\theta) e^{i m\phi} = \sqrt{\frac{2 l +1}{4 \pi}\frac{(l-m)!}{(l+m)!}}P_l^m({\cos \theta}) e^{i m\phi} \\
- {\cal L}_n(r) &=& \sqrt{\frac{n!}{(n+\alpha)!}} e^{-r/2} L_n^{(\alpha)}(r)
+ {\cal K}_n(r,\tau) &=& \sqrt{\frac{n!}{(n+\alpha)!}} \frac{e^{-r/2\tau} L_n^{(\alpha)}(r/\tau)}{\tau^{3/2}}
 \end{eqnarray}$$
 
-with $ L_n^{(\alpha)}(r)$ the generalized Laguerre polynomials (notice that by default $ \alpha = 2$). 
+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
 $$
-f(r) = \sum_{n=0}^{N-1} f_n\ {\cal L}_n^{(\alpha)}(r)
+f(r) = \sum_{n=0}^{N-1} f_n\ {\cal K}_n(r)
 $$
 and the r-values are the roots of the generalized Laguerre polynomial of order $N$ noted $r_i$ with $i$ in $[0,N-1]$. From the orthogonality of the Laguerre functions, it yields the Gauss-Laguerre quadrature
 $$
@@ -58,19 +58,32 @@ The analysis process starts by gathering the 3D pixels real values $f_{ijk} \equ
 	+ for the $k$th-radial shell, one performs a Spherical Harmonic decomposition using for instance the `map2alm` routine of `libsharp` to get the set of complex coefficients $a_{lmk}=a_{lm}(r_k)$ with $l\in\{0,\dots,L-1\}$ and $m \in\{0,\dots,l\}$ (the notation remind the underlying 2D Spherical Harmonic transform and the usual way to note the resulting coefficients $a_{lm}$);
 	+ for each $(l,m)$, one uses the set of values $a_{lmk}$ with $k\in\{0,\dots,N-1\}$ to determine the $f_{lmn}$ complex coefficients thanks to the following simplified relation
 $$
-	f_{lmn} = \sum_{k=0}^{N-1}w_k\ a_{lmk}\ {\cal L}_n^{(\alpha)}(r_k) \quad\quad n\in\{0,\dots,N-1\}
+	f_{lmn} = \sum_{k=0}^{N-1}w_k\ a_{lmk}\ {\cal K}_n(r_k) \quad\quad n\in\{0,\dots,N-1\}
 $$
 
 + **Synthesis**:  
 The synthesis also proceeds in two steps as the above detailed analysis. The algorithm receives as input the $f_{lmn}$ complex coefficients  
 	+ to first compute for each $(l,m)$ couple of indices the intermediate $a_{lmk}$ complex values obtained from the Inverse Laguerre Transform:
 $$
-a_{lmk} = \sum_{n=0}^{N-1} f_{lmn}\ {\cal L}_n^{(\alpha)}(r_k)
+a_{lmk} = \sum_{n=0}^{N-1} f_{lmn}\ {\cal K}_n(r_k)
 $$
  	+ then, thanks to the Inverse Spherical Harmonic Transform, `alm2map` routine of `libsharp`, ones determines the real $f_{ijk}$ pixel values.
 
 We take advantage from the Matrix Multiplication writing of the $f_{lmn} \leftrightarrow a_{lmk}$ passage to use efficient algorithm or even more efficiently the BLAS-like libraries (OpenBLAS for Linux and the native Accelerate framework on Mac OS X). 
 
+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
+\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
+$$
+\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)
+\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 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. 
@@ -83,4 +96,6 @@ with $r_k$ the $k$-th roots of the Laguerre polynomial of order $N$.
 
 + The *lagsht_spheregeom.h* file includes the different 2D-Sphere pixelization schemes available derived from the LagSHT::BaseGeometry class. This class provides access to the pixelization throw the methods **PixIndexSph** and **PixThetaPhi**. The first method provides the pixel index at a given $(\theta, \phi)$ position, while the second method does the reverse by providing the $(\theta, \phi)$ position of a given pixel index.
 
-+ The LagSHT::LaguerreSphericalTransform class performs the combinaison of the Laguerre transform and the Spherical Harmonic transform. 
\ No newline at end of file
++ 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`.
\ No newline at end of file