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

(JEC) 4/515 add Documentation.md which gives some hints on the underlying math...

(JEC) 4/515 add Documentation.md which gives some hints on the underlying math and some explanations on the algo and main classes
parent 81b33f0c
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.
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).
Mathematics in summary
----------------
Any function \f$f(r,\Omega)\f$ square-integrable on \f$\mathrm{B}^3 = \mathrm{R}^+ \times [0,\pi] \times [0,2\pi)\f$ can be decomposed as:
\f{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)
\label{eq:FLagfullb}
\f}
based on the orthogonality of the set of functions
\f{equation}
K_{lmn}(r,\Omega) \equiv Y_{l,m}(\Omega)\times {\cal L}_n(r)
\label{eq:FLagOrthoFunc}
\f}
using the spherical harmonic functions \f$Y_{l,m}\f$ and the genralized Laguerre functions \f${\cal L}_n\f$
\f{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)
\f}
with \f$ L_n^{(\alpha)}(r)\f$ the generalized Laguerre polynomials (notice that by default \f$ \alpha = 2\f$).
For bnad-limited functions the discret sums are truncated with \f$ l <L\f$ and \f$ n < N\f$ (notice the L band-limit definition is not the same used for instance in `libsharp`). Focusing on the radial part only, one gets
\f[
f(r) = \sum_{n=0}^{N-1} f_n\ {\cal L}_n^{(\alpha)}(r)
\f]
and the r-values are the roots of the generalized Laguerre polynomial of order \f$N\f$ noted \f$r_i\f$ with i in [0,N-1]. From the orthogonality of the Laguerre functions, it yields the Gauss-Laguerre quadrature
\f[
f_n = \sum_{i=0}^{N-1} w_i\ f(r_i)\ {\cal L}^{(\alpha)}_n (r_i) \quad\quad n\in \{ 0,\dots, N-1\}
\f]
with the weight
\f[
w_i = \frac{(N + \alpha)!}{N!(N+1)^2}\frac{r_i e^{r_i}}{\left[L^{(\alpha)}_{N+1}(r_i)\right]^2}
\f]
The Spherical Harmonic part is done by `libsharp` and the reader is invited to have a look at the reference mentionned in the introduction.
The combined transformation (spherical and Laguerre) can be summarized by the following schema:
![img def](../fig-schema.png "Schema of global Analysis and Synthesis")
+ **Analysis**:
The analysis process starts by gathering the 3D pixels real values \f$f_{ijk} \equiv f(\theta_i, \phi_j, r_k)\f$
where \f$r_k\f$ (\f$k\in\{0,\dots,N-1\}\f$) are the \f$N\f$ Laguerre zeros of \f$L_N(x)\f$ and \f$(\theta_i, \phi_j)\f$
are the 2D angles defined by a pixelization scheme.
For instance, one can use the ECP Fejer's sum rule with \f$N_\theta = 2L-1\f$ colatitude rings of \f$N_\phi \geq 2L-1\f$ iso-latitude pixels.
Then, one proceeds as followed
+ for the \f$k\f$th-radial shell, one performs a Spherical Harmonic decomposition using for instance the `map2alm` routine of `libsharp` to get the set of complex coefficients \f$f_{lmk}\f$ with \f$l\in\{0,\dots,L-1\}\f$ and \f$m \in\{0,\dots,l\}\f$;
+ for each \f$(l,m)\f$, one uses the set of values \f$f_{lmk}\f$ with \f$k\in\{0,\dots,N-1\}\f$ to determine the \f$f_{lmn}\f$ complex coefficients thanks to the following simplified relation
\f[
f_{lmn} = \sum_{k=0}^{N-1}w_k\ f_{lmk}\ {\cal L}_n(r_k) \quad\quad n\in\{0,\dots,N-1\}
\f]
+ **Synthesis**:
The synthesis also proceeds in two steps as the above detailed analysis. The algorithm receives as input the \f$f_{lmn}\f$ complex coefficients
+ to first compute for each \f$(l,m)\f$ couple of indices the intermediate \f$f_{lmk}\f$ complex values obtained from the Inverse Laguerre Transform:
\f[
f_{lmk} = \sum_{n=0}^{N-1} f_{lmn}\ {\cal L}_n^{(\alpha)}(r_k)
\f]
+ then, thanks to the Inverse Spherical Harmonic Transform, `alm2map` routine of `libsharp`, ones determines the real \f$f_{ijk}\f$ pixel values.
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::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) \f$R_k\f$:
\f[
R_k = R_{max} \times \frac{r_k}{r_{N-1}}
\f]
with \f$r_k\f$ 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 \f$(\theta, \phi)\f$ position, while the second method does the reverse by providing the \f$(\theta, \phi)\f$ 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
......@@ -774,7 +774,7 @@ INPUT_ENCODING = UTF-8
# *.md, *.mm, *.dox, *.py, *.f90, *.f, *.for, *.tcl, *.vhd, *.vhdl, *.ucf,
# *.qsf, *.as and *.js.
FILE_PATTERNS = *.h *.cc
FILE_PATTERNS = *.h *.cc *.md
# The RECURSIVE tag can be used to specify whether or not subdirectories should
# be searched for input files as well.
......@@ -789,7 +789,7 @@ RECURSIVE = NO
# Note that relative paths are relative to the directory from which doxygen is
# run.
EXCLUDE =
EXCLUDE = README.md History.md
# The EXCLUDE_SYMLINKS tag can be used to select whether or not files or
# directories that are symbolic links (a Unix file system feature) are excluded
......@@ -1079,7 +1079,7 @@ HTML_STYLESHEET =
# see the documentation.
# This tag requires that the tag GENERATE_HTML is set to YES.
HTML_EXTRA_STYLESHEET =
HTML_EXTRA_STYLESHEET = doxygen_html_style.css
# The HTML_EXTRA_FILES tag can be used to specify one or more extra images or
# other source files which should be copied to the HTML output directory. Note
......@@ -1397,7 +1397,7 @@ EXT_LINKS_IN_WINDOW = NO
# Minimum value: 8, maximum value: 50, default value: 10.
# This tag requires that the tag GENERATE_HTML is set to YES.
FORMULA_FONTSIZE = 10
FORMULA_FONTSIZE = 14
# Use the FORMULA_TRANPARENT tag to determine whether or not the images
# generated for formulas are transparent PNGs. Transparent PNGs are not
......
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