📣 An issue occured with the embedded container registry on October 25 2021, between 10:30 and 12:10 (UTC+2). Any persisting issues should be reported to CC-IN2P3 Support. 🐛

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

Documentation style according to GitLab std

parent 6a4661dc
......@@ -5,7 +5,7 @@ Introduction
------------
**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:
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:
* we use the `libsharp` library ([see the ArXiv paper by M. Reinecke and D. S. Seljebotn](http://arxiv.org/pdf/1303.4945v2.pdf)) to deal with the Spherical Harmonic Transform part and use possibly different pixelization schemes (ECP/Fejer1 sum rule, Gauss, Healpix). Notice that user pixelization which are available in libsharp are not yet interfaced in LagSHT.
* we have also extended the numerical viability of the Laguerre Transform using simple 64-bits floating-point representation thanks to an on fly re-scaling technics to compute the Laguerre polynomial values. The same technics was used in the `libsharp` library to compute the Wigner d-matrix coefficients. This is particularly important to compute the weights and nodes of the Laguerre Quadrature that we can manage easily with N ~ 20,000 (we do not study exactly up to which N the computation fails).
......@@ -16,40 +16,45 @@ Mathematics in summary
Any real 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; \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; \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 K}_n$:
$$\begin{eqnarray}
```math
\begin{array}{rcl}
f(r,\Omega) &=& \sum_{n=0}^\infty \sum_{l=0}^{\infty}\sum_{m=-l}^{l}\ f_{lmn}\ K_{lmn}(r,\Omega; \tau)\\
f_{lmn} &=& \int_{\mathrm{B}^3} dr d\Omega\ r^2\ f(r,\Omega)\ K^\ast_{lmn}(r,\Omega; \tau)
\end{array}
```
based on the orthogonality of the set of functions
```math
K_{lmn}(r,\Omega; \tau) \equiv Y_{l,m}(\Omega)\times \mathcal{K}_n(r,\tau)
````
using the spherical harmonic functions $Y_{l,m}(\Omega)$ and the genralized Laguerre functions $\mathcal{K}_n$:
```math
\begin{array}{rcl}
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 K}_n(r,\tau) &=& \tau^{-3/2} \sqrt{\frac{n!}{(n+\alpha)!}} e^{-r/2\tau} \left(\frac{r}{\tau}\right)^{\frac{\alpha}{2}-1}
L_n^{(\alpha)}(r/\tau) = \tau^{-3/2} {\cal K}_n(r/\tau,1)
\end{eqnarray}$$
\mathcal{K}_n(r,\tau) &=& \tau^{-3/2} \sqrt{\frac{n!}{(n+\alpha)!}} e^{-r/2\tau} \left(\frac{r}{\tau}\right)^{\frac{\alpha}{2}-1}
L_n^{(\alpha)}(r/\tau) = \tau^{-3/2} \mathcal{K}_n(r/\tau,1)
\end{array}
```
with $ L_n^{(\alpha)}(r)$ the generalized Laguerre polynomials (notice that by default $ \alpha = 2$ but one is free to use real numbers). 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 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
$$\begin{eqnarray}
f(R_i) &=& f(\tau r_i) = \sum_{n=0}^{N-1} f_n\ {\cal K}_n(r_i,\tau) \\
```math
\begin{array}{rcl}
f(R_i) &=& f(\tau r_i) = \sum_{n=0}^{N-1} f_n\ \mathcal{ K}_n(r_i,\tau) \\
&=& \tau^{-3/2} \sum_{n=0}^{N-1} f_n \sqrt{\frac{n!}{(n+\alpha)!}} e^{-r_i/2} L^{(\alpha)}_n(r_i) r_i^{\frac{\alpha}{2}-1}
\end{eqnarray}$$
\end{array}
```
The function is evaluated at $R_i$ values which are scaled values of the roots of the generalized Laguerre polynomial of order $N$ noted $r_i$ with $i$ in $[0,N-1]$. To get the $f_n$ coefficients one uses Gauss-Laguerre quadrature to compute
$$
f_n = \int_{R^+} dR\ R^2\ f(R)\ {\cal K}_n(R,\tau)
f_n = \int_{R^+} dR\ R^2\ f(R)\ \mathcal{ K}_n(R,\tau)
$$
It yields
$$
f_n = \tau^{3/2} \sqrt{\frac{n!}{(n+\alpha)!}} \sum_{i=0}^{N-1} w_i\ f(\tau r_i)\ e^{-r_i/2} L^{(\alpha)}_n(r_i) r_i^{1-\frac{\alpha}{2}}\quad\quad n\in \{ 0,\dots, N-1\}
f_n = \tau^{3/2} \sqrt{\frac{n!}{(n+\alpha)!}} \sum_{i=0}^{N-1} w_i\ f(\tau r_i)\ e^{-r_i/2} L^{(\alpha)}_n(r_i) r_i^{1-\frac{\alpha}{2}}\quad\quad n \in \{0, ..., N-1\}
$$
with the weight *(scaled with $e^{r_i}$ for numerical purpose as it is included in the computation by recurence of Laguerre polynomials)*
$$
w_i = \frac{(N + \alpha)!}{N!(N+1)^2}\frac{r_i e^{r_i}}{\left[L^{(\alpha)}_{N+1}(r_i)\right]^2}
......@@ -61,54 +66,53 @@ The combined transformation (spherical and Laguerre) can be summarized by the fo
![img def](./fig-schema.png "Schema of global Analysis and Synthesis")
+ **Analysis**:
The analysis process starts by gathering the 3D pixels real values $f_{ijk} \equiv f(\theta_i, \phi_j, r_k)$ where $r_k$ ($k\in\{0,\dots,N-1\}$) are the $N$ Laguerre zeros of $L_N(x)$ and $(\theta_i, \phi_j)$ are the 2D angles defined by a pixelization scheme. For instance, one can use the ECP Fejer's sum rule with $N_\theta = 2L-1$ colatitude rings of $N_\phi \geq 2L-1$ iso-latitude pixels.
The analysis process starts by gathering the 3D pixels real values $f_{ijk} \equiv f(\theta_i, \phi_j, r_k)$ where $r_k$ ($k\in\{0,...,N-1\}$) are the $N$ Laguerre zeros of $L_N(x)$ and $(\theta_i, \phi_j)$ are the 2D angles defined by a pixelization scheme. For instance, one can use the ECP Fejer's sum rule with $N_\theta = 2L-1$ colatitude rings of $N_\phi \geq 2L-1$ iso-latitude pixels.
Then, one proceeds as followed
+ 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
+ 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,...,L-1\}$ and $m \in\{0,...,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,...,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(r_k) \quad\quad n\in\{0,\dots,N-1\}
f_{lmn} = \sum_{k=0}^{N-1}w_k\ a_{lmk}\ \mathcal{ K}_n(r_k) \quad\quad n\in\{0,...,N-1\}
$$
*(here to simplify the notation I have introduced the ${\cal L}$ function but it is derived from the Gauss-Quadrature expression given above)*
*(here to simplify the notation I have introduced the $\mathcal{ K}$ function but it is derived from the Gauss-Quadrature expression given above)*
+ **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 K}_n(r_k)
a_{lmk} = \sum_{n=0}^{N-1} f_{lmn}\ \mathcal{ 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).
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).
### Spin > 0
Since **v2.0** it is introduced the spin-weighted decomposition of complex functions. In practice it is restricted to spin up to 2. The Spherical decomposition uses the *gradient* **E** and *curl* **B** coeficients in place of the above $a_{lmk}$ and the output of the Analysis are the $E_{lmn}$ and $B_{lmn}$ complex coefficients. For praticle usage of the ligsharp library C++ interface one should feed sperately the real and imaginary parts of the complex pixelized function.
### Fourier-Bessel (spin = 0)
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}
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}$:
```math
\begin{array}{rcl}
f_{lmp}^{FB} &=& \sum_n f_{lmn}^{FL} J_{ln}(k_{lp})\\
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}$.
J_{lnp}\equiv J_{ln}(k_{lp}) &=& \sqrt{\frac{2}{\pi}}\int_0^\infty j_l(k_{lp}r)\mathcal{ K}_n(r,\tau) r^2 dr
\end{array}
```
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}$.
####For Version < 2.1
##### For Version < 2.1
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=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}
$$
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 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}$.
#### Version 2.1
##### Version 2.1
The computation of the $J_{lnp}$ is quite CPU intensive and so I have developped a new method suggested by Jason Mc Ewen. In brief the $J_{lnp}$ is an integral of the product of two oscillatory functions. An idea is to perform a Chebyshev Decomposition of both functions in a basis of order N, then project the product in the basis of order 2*N and use the Clenshaw-Curtis weights to compute the integral. The gain is that the Chebyshev decomposition and its invserse can be performed by Fast Discrete Cosine Transform as well as the computation of the Clenshaw-curtis weights. To minimize the number of FFT planning which increases a little the overhead, one uses estimations of the number of zeros and extrema of both Bessel function and Laguerre polynomials, to get the order N of the Chebyshev decomposition.
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.
+ 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.
+ The LagSHT::LaguerreTransform class implement the radial Laguerre Transformation part in the forward (Analysis) and backward (Synthesis) direction. The class also provides access to the radial pixelization (methods **R2Index** and **Index2R**) as the k-th shell has the radius (in physical units) $R_k$:
$$
......@@ -118,6 +122,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.
+ 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` as well as a Discrete Fourier-Bessel Transform.
\ No newline at end of file
+ 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 preview for this file type
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