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

Merge branch 'master' of gitlab.in2p3.fr:campagne/LagSHT

parents 930cd612 c87f4555
# Flag preprocesseur
CPPFLAGS = -DDarwin -I/opt/local/include
# Def compilateur C et flags
CC = clang
# Replace [CNFPHFLF] by specific compilation flags below (done by configure)
CFLAGS = -fno-common -O
# With optimization flags
# CFLAGS = -fno-common -g -fastf -mtune=G5 -fPIC
# Def compilateur C++ et flags
CXX = clang++
#JEC add -fopenmp -O3 -Wa,-q
CXXFLAGS = -fno-common -g -O3 -fopenmp -Wall -std=c++11 -ffast-math -m64
# flags specifiques pour templates repository...
CXXTEMPFLG =
# Compilo fortran
FC = gfortran
FFLAGS =
# Pour fabriquer les lib .a
AR = libtool
ARFLAGS = -static -o
# Ordres de compilation
CCOMPILE = $(CC) $(CPPFLAGS) $(CFLAGS) -c
CXXCOMPILE = $(CXX) $(CPPFLAGS) $(CXXFLAGS) -c
FCCOMPILE = $(FC) $(FFLAGS) -c
# Extension de nom pour les librairies partagees
SLEXT = dylib
# Fabrication de librairies partagees
CMDSHLCXX = $(CXX) $(CXXFLAGS) $(CXXTEMPFLG) -dynamiclib
LIBFORT =
# Commande de link pour creation d'executables
CXXLINK = $(CXX) $(CXXFLAGS) -bind_at_load
# Commande de link pour creation de module pouvant etre charge dynamiquement
CXXDLL = $(CXX) $(CXXFLAGS) -bundle
......@@ -11,7 +11,7 @@ CFLAGS = -fno-common -g -O -fPIC
# Def compilateur C++ et flags
CXX = g++
#JEC add -fopenmp -O3 -Wa,-q
CXXFLAGS = -fno-common -g -O3 -fPIC -fopenmp -Wa,-q -Wall -std=c++11
CXXFLAGS = -fno-common -g -O3 -fPIC -fopenmp -Wa,-q -Wall -Wpedantic -Wno-deprecated-declarations -std=c++11 -march=native -ffast-math -m64
# flags specifiques pour templates repository...
CXXTEMPFLG =
# Compilo fortran
......
......@@ -18,6 +18,7 @@ ifeq ($(MACH), Linux)
BOOSTINC = -I$(BOOSTDIR)
else
include Darwin_g++_omp_make.inc
# include Darwin_clang_omp_make.inc
# ===== nNative BLAS (Darwin)
BLASYES = -DCBLAS -flax-vector-conversions
BLASLIBN = -framework Accelerate
......@@ -39,6 +40,12 @@ SHARPLIB = $(SHARPDIR)/lib
SHARPINC = -I$(SHARPDIR)/include
SHARPLIBN = -L$(SHARPLIB) -lsharp -lc_utils -lfftpack
# ===== FFTW3 ======
FFTWDIR = /opt/local
FFTWLIB = $(FFTWDIR)/lib
FFTWINC = -I$(FFTWDIR)/include
FFTWLIBN = -L$(FFTWLIB) -lfftw3 -lfftw3_threads
SRC = ./src/
LIB = ./lib/
OBJ = ./objs/
......@@ -99,8 +106,8 @@ tidy :
.PHONY: clean
clean :
rm -f $(OBJ)/*
rm -f $(EXE)/*
rm -f $(LIB)/*
rm -rf $(EXE)/*
rm -rf $(LIB)/*
#C++ common Objects
CXXOBJ = $(OBJ)laguerreBuilder.o \
......@@ -111,6 +118,7 @@ CXXOBJ = $(OBJ)laguerreBuilder.o \
$(OBJ)lagsht_bessel.o \
$(OBJ)laguerre2bessel.o
CXXSHOBJ = laguerreBuilder.o \
laguerreTransform.o \
lagsht_spheregeom.o \
......@@ -133,10 +141,11 @@ CXXHDR = $(SRC)lagsht_exceptions.h \
$(SRC)walltimer.h \
$(SRC)lagsht_bessel.h \
$(SRC)laguerre2bessel.h \
$(SRC)lagsht_func.h \
$(SRC)quadinteg.h
CPPFLAGS += $(BLASYES) $(SHARPINC) $(BLASINC) $(BOOSTINC)
LDFLAGS += $(SHARPLIBN) $(BLASLIBN) -lm
CPPFLAGS += $(BLASYES) $(SHARPINC) $(BLASINC) $(BOOSTINC) $(FFTWINC)
LDFLAGS += $(SHARPLIBN) $(BLASLIBN) $(FFTWLIBN) -lm
#C++ rule for compiling
......@@ -144,6 +153,7 @@ $(OBJ)%.o: $(SRC)%.cc $(CXXHDR)
echo "compile... $<"
$(CXXCOMPILE) $< -o $@
######################
.PHONY: sharelib
sharelib : $(CXXOBJ)
......
Laguerre Spherical Harmonic Transform
======================================
Version 1.x (last update of this file 7/10/15)
Version 2.x (last update of this file 6/04/16)
LagSHT performs Laguerre Spheraical Harmonic transform over different set of 2D-sphere mappings, choice of Generzalized Laguerre parameter ($\alpha$) and spin (eg. 0 for scalar, 2 for polarization). The Bessel coefficients are also computed from Laguerre ones although this part is more experimental.
Main Author: J.E Campagne - Lab. Accelerateur Lineaire - CNRS/IN2P3 & PARIS-SUD University - France
......@@ -18,6 +20,7 @@ git clone https://gitlab.in2p3.fr/campagne/LagSHT.git
* Makefile
* Darwin_g++_omp_make.inc : file included into Makefile in case of running on Mac Os X system
* Linux_g++_make.inc : idem but for Linux system
* setup.[c]sh to define some environment variable(s) used in lagsht_testsuite.cc
* **src** directory
* laguerreBuilder.h (.cc) : class of generalized Laguerre function
* laguerreTransform.h (.cc) : class to manage the Laguerre function Transform (Synthesis & Analysis): Single or Multiple
......@@ -31,6 +34,7 @@ git clone https://gitlab.in2p3.fr/campagne/LagSHT.git
* laguerre2bessel.h (.cc) : implement the computation of Fourier-Bessel transform from Fourier-Laguerre coeff.
* lagsht_bessel.h (.cc) : simple Bessel class (use Boost)
* quadinteg.h : Quadrature integration and strategy
* lagsht_func.h : some utility classes for 1D function
* **root** directory
* Makefile that should be tuned to the local platform (*.inc files)
* **doc** directory
......@@ -62,6 +66,11 @@ If the variable `CBLAS` is defined then the code of laguerreTransform.cc use the
This is simply header files which are needed to access to the `"boost/math/special_functions/bessel.hpp"` and related headers of the library.
# Compilation/Installation/Setup
0. setup:
* source setup.sch in Clike-shell
* . setup.sh in SHlike-shell
This is to define **LAGSHTROOTDIR**
1. edit Makefile and adapat to local platform :
* adapt to the type of Machine MacOSX (Darwin) vs Linux
......@@ -99,19 +108,21 @@ This is simply header files which are needed to access to the `"boost/math/speci
Similarly, the library is in <root>/lib but one can switch to
another location by modifying the LIB variable in the
Makefile.
6. execute setup.[c]sh depending on your shell to setup **LAGSHTROOTDIR** environment variable.
# Plateform tested
Mac OS X 10.9.5 + gcc 4.8.4
Mac OS X 10.11.2 + gcc 4.9.3
Linux SLC 6.6 + gcc 4.9.1 20140922
# Running
> ./Objs/lagsht_testsuite -t <test number> [1]
[-n <Nmax value> [5]] [-l <Lmax value> [10]]
[-a <alpha parameter> [2]]
[-n <Nmax value> [128]] [-l <Lmax value> [128]]
[-g <geometry> Gauss|Fejer1|Healpix [Gauss]]
[-ntheta <number of theta rings> [determined by the geometry]
in case of Healpix gives the nside parameter]
[-nphi <number of pixel per rings> [determined by the geometry]]
[-spin <spin value: 0, 1, 2> [0]]
> <test number>:
0: basic test for numerical application
1: to get the nodes & weights of the generalized Gauss-laguerre function quadrature.
......
......@@ -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,45 +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}$. 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
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
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
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$:
$$
......@@ -109,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
......@@ -15,9 +15,9 @@
namespace LagSHT {
void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& flmn,
vector<r_8>& fijk,
vector< complex<r_8> >& flmk) {
void LaguerreSphericalTransform::Synthesis(const std::vector< std::complex<r_8> >& flmn,
std::vector<r_8>& fijk,
std::vector< std::complex<r_8> >& flmk) {
int Ntot3DPix = Npix_*N_;
int NtotCoeff = Nalm_*N_;
if(flmn.size() != (size_t)NtotCoeff) throw LagSHTError("SphLagTransform::Synthesis size error");
......@@ -25,7 +25,6 @@ void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& flmn,
fijk.resize(Ntot3DPix);
//Perform the Inverse Laguerre Transform first
// vector< complex<r_8> > flmk(NtotCoeff);
flmk.resize(NtotCoeff);
......@@ -52,7 +51,7 @@ void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& flmn,
off7k = n*Npix_;
#if DEBUG >= 2
cout << "Synthesis... at shell num: " << n << endl;
std::cout << "Synthesis... at shell num: " << n << std::endl;
#endif
sphere_->alm2map(reinterpret_cast<const r_8*>(&(flmk.data()[off7n])),&fijk.data()[off7k],false); //the false is to forbidd accumulation
......@@ -61,15 +60,15 @@ void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& flmn,
tstack_pop("SHT Synthesis");
#if DEBUG >= 1
cout << "SphLagTransform::Synthesis.... Done" << endl;
std::cout << "SphLagTransform::Synthesis.... Done" << std::endl;
#endif
}//Synthesis
void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk,
vector< complex<r_8> >& flmn,
vector< complex<r_8> >& almk
void LaguerreSphericalTransform::Analysis(const std::vector<r_8>& fijk,
std::vector< std::complex<r_8> >& flmn,
std::vector< std::complex<r_8> >& almk
) {
int Ntot3DPix = Npix_*N_;
int NtotCoeff = Nalm_*N_;
......@@ -79,7 +78,7 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk,
//done in Laguerre routine flmn.resize(NtotCoeff);
#if DEBUG >= 1
cout << "SphLagTransform::Analysis.... Go..." << endl;
std::cout << "SphLagTransform::Analysis.... Go..." << std::endl;
#endif
//Perform the Spherical Hamonic analysis first
......@@ -93,7 +92,7 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk,
off7k = n*Npix_;
#if DEBUG >= 1
cout << "Start Harmonic Spheric... at n = "<< n << endl;
std::cout << "Start Harmonic Spheric... at n = "<< n << std::endl;
#endif
sphere_->map2alm(&(fijk.data()[off7k]),reinterpret_cast<r_8*>(&(almk.data()[off7n])),false); //the false is to forbidd accumulation
......@@ -101,8 +100,8 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk,
}//loop on shell
#if DEBUG >= 1
cout << "End Harmonic Spheric..." << endl;
cout << "Start Laguerre..." << endl;
std::cout << "End Harmonic Spheric..." << std::endl;
std::cout << "Start Laguerre..." << std::endl;
#endif
tstack_pop("SHT Analysis");
......@@ -112,7 +111,7 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk,
lagTrans_->MultiAnalysis(almk, flmn, Nalm_);
#if DEBUG >= 1
cout << "End Laguerre..." << endl;
std::cout << "End Laguerre..." << std::endl;
#endif
tstack_pop("Laguerre MultiAnalysis");
......@@ -120,20 +119,22 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk,
#if DEBUG >= 2
for (int i=0;i<NtotCoeff; i++) {
cout << "flmn("<<i<<"): " <<flmn[i] << endl;
std::cout << "flmn("<<i<<"): " <<flmn[i] << std::endl;
}
#endif
#if DEBUG >= 1
cout << "SphLagTransform::Analysis.... Done" << endl;
std::cout << "SphLagTransform::Analysis.... Done" << std::endl;
#endif
}//Analysis
void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& Elmn, const vector<complex<r_8> >& Blmn,
int spin,
vector<r_8>& fijk_re, vector<r_8>& fijk_im,
vector< complex<r_8> >& Elmk, vector<complex<r_8> >& Blmk
void LaguerreSphericalTransform::Synthesis(const std::vector< std::complex<r_8> >& Elmn,
const std::vector< std::complex<r_8> >& Blmn,
int spin,
std::vector<r_8>& fijk_re, std::vector<r_8>& fijk_im,
std::vector< std::complex<r_8> >& Elmk,
std::vector< std::complex<r_8> >& Blmk
) {
int Ntot3DPix = Npix_*N_;
int NtotCoeff = Nalm_*N_;
......@@ -148,8 +149,8 @@ void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& Elmn, c
Blmk.resize(NtotCoeff);
#if DEBUG >= 1
cout << "SphLagTransform::Synthesis type <r_8> I/O .... Go..." << endl;
cout << "Start Inverse Laguerre...." << endl;
std::cout << "SphLagTransform::Synthesis type <r_8> I/O .... Go..." << std::endl;
std::cout << "Start Inverse Laguerre...." << std::endl;
#endif
tstack_push("Laguerre MultiSynthesis");
......@@ -157,17 +158,9 @@ void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& Elmn, c
lagTrans_->MultiSynthesis(Blmn, Blmk, Nalm_);
tstack_pop("Laguerre MultiSynthesis");
// cout << "Elmk, Blmk" <<endl;
// for(int i=0; i<100; i++){
// cout << "Elmk["<<i<<"]: " << Elmk[i] << "; "
// << "Blmk["<<i<<"]: " << Blmk[i] << endl;
// }
#if DEBUG >= 1
cout << "End Inverse Laguerre...." << endl;
cout << "Start Inverse H.Spherical Transform " << endl;
std::cout << "End Inverse Laguerre...." << std::endl;
std::cout << "Start Inverse H.Spherical Transform " << std::endl;
#endif
tstack_push("SHT Synthesis");
......@@ -178,7 +171,7 @@ void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& Elmn, c
off7k = n*Npix_;
#if DEBUG >= 2
cout << "Synthesis... at shell num: " << n << endl;
std::cout << "Synthesis... at shell num: " << n << std::endl;
#endif
sphere_->alm2map_spin(reinterpret_cast<const r_8*>(&(Elmk.data()[off7n])),
reinterpret_cast<const r_8*>(&(Blmk.data()[off7n])),
......@@ -192,16 +185,19 @@ void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& Elmn, c
tstack_pop("SHT Synthesis");
#if DEBUG >= 1
cout << "SphLagTransform::Synthesis.... Done" << endl;
std::cout << "SphLagTransform::Synthesis.... Done" << std::endl;
#endif
}//Synthesis
void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk_re, const vector<r_8>& fijk_im,
int spin,
vector< complex<r_8> >& Elmn, vector<complex<r_8> >& Blmn,
vector< complex<r_8> >& Elmk, vector<complex<r_8> >& Blmk
) {
void LaguerreSphericalTransform::Analysis(const std::vector<r_8>& fijk_re,
const std::vector<r_8>& fijk_im,
int spin,
std::vector< std::complex<r_8> >& Elmn,
std::vector< std::complex<r_8> >& Blmn,
std::vector< std::complex<r_8> >& Elmk,
std::vector< std::complex<r_8> >& Blmk
) {
int Ntot3DPix = Npix_*N_;
int NtotCoeff = Nalm_*N_;
if(fijk_re.size() != (size_t)Ntot3DPix || fijk_im.size() != (size_t)Ntot3DPix) throw LagSHTError("SphLagTransform::Analysis size error");
......@@ -215,7 +211,7 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk_re, const vect
#if DEBUG >= 1
cout << "SphLagTransform::Analysis.... Go..." << endl;
std::cout << "SphLagTransform::Analysis.... Go..." << std::endl;
#endif
//Perform the Spherical Hamonic analysis first
......@@ -229,7 +225,7 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk_re, const vect
off7k = n*Npix_;
#if DEBUG >= 1
cout << "Start Harmonic Spheric... at n = "<< n << endl;
std::cout << "Start Harmonic Spheric... at n = "<< n << std::endl;
#endif
sphere_->map2alm_spin(&(fijk_re.data()[off7k]),
......@@ -243,8 +239,8 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk_re, const vect
}//loop on shell
#if DEBUG >= 1
cout << "End Harmonic Spheric..." << endl;
cout << "Start Laguerre..." << endl;
std::cout << "End Harmonic Spheric..." << std::endl;
std::cout << "Start Laguerre..." << std::endl;
#endif
tstack_pop("SHT Analysis");
......@@ -255,18 +251,18 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk_re, const vect
lagTrans_->MultiAnalysis(Blmk, Blmn, Nalm_);
#if DEBUG >= 1
cout << "End Laguerre..." << endl;
std::cout << "End Laguerre..." << std::endl;
#endif
tstack_pop("Laguerre MultiAnalysis");
#if DEBUG >= 2
for (int i=0;i<NtotCoeff; i++) {
cout << "flmn("<<i<<"): " << Elmn[i] << ", " << Blmn[i] << endl;
std::cout << "flmn("<<i<<"): " << Elmn[i] << ", " << Blmn[i] << std::endl;
}
#endif
#if DEBUG >= 1
cout << "SphLagTransform::Analysis.... Done" << endl;
std::cout << "SphLagTransform::Analysis.... Done" << std::endl;
#endif
}//Analysis
......
......@@ -54,7 +54,7 @@ class LaguerreSphericalTransform {
// LaguerreSphericalTransform(string spheregeom, int Lmax, int Mmax, int Nmax,
// r_8 R, int Nrings = -1, int Nphi = -1, int alpha=2, bool R2tau=true):
LaguerreSphericalTransform(string spheregeom, int Lmax, int Mmax, int Nmax,
LaguerreSphericalTransform(std::string spheregeom, int Lmax, int Mmax, int Nmax,
r_8 R, int Nrings = -1, int Nphi = -1, r_8 alpha=2, bool R2tau=true):
Lmax_(Lmax), Mmax_(Mmax), R_(R) {
......@@ -66,7 +66,7 @@ class LaguerreSphericalTransform {
} else if (spheregeom == "Healpix") {
sphere_ = new HealpixGeometry(Lmax, Mmax, Nrings); //here Nsides = Nrings
} else {
cout << "FAILURE: LaguerreSphericalTransform::Ctor wrong geometry type" << endl;
std::cout << "FAILURE: LaguerreSphericalTransform::Ctor wrong geometry type" << std::endl;
sphere_ = 0;
throw LagSHTError("LaguerreSphericalTransform::Ctor wrong geometry type");
}
......@@ -97,8 +97,10 @@ class LaguerreSphericalTransform {
//! Destructor
virtual ~LaguerreSphericalTransform() {
if(sphere_) delete sphere_; sphere_ = 0;
if(lagTrans_) delete lagTrans_; lagTrans_ = 0;
if(sphere_) delete sphere_;
sphere_ = 0;
if(lagTrans_) delete lagTrans_;
lagTrans_ = 0;
}
//! Geometry helper
......@@ -119,12 +121,12 @@ class LaguerreSphericalTransform {
//! Get the f(l,m,n) coefficient in the vector collection
//make a Flmn Class
inline complex<r_8> GetFlmn(const vector< complex<r_8> >& flmn, int l , int m, int n) const {
inline std::complex<r_8> GetFlmn(const std::vector< std::complex<r_8> >& flmn, int l , int m, int n) const {
int id= n*Nalm_ + l +m*Lmax_ -m*(m+1)/2;
return flmn[id];
}
//! Set the f(l,m,n) coefficient in the vector collection
inline void SetFlmn(vector< complex<r_8> >& flmn, int l , int m, int n, complex<r_8> value) const {
inline void SetFlmn(std::vector< std::complex<r_8> >& flmn, int l , int m, int n, std::complex<r_8> value) const {
int id= n*Nalm_ + l +m*Lmax_ -m*(m+1)/2;
flmn[id] = value;
}
......@@ -137,18 +139,18 @@ class LaguerreSphericalTransform {
\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,
vector< complex<r_8> >& almk);
void Synthesis(const std::vector< std::complex<r_8> >& flmn,
std::vector<r_8>& fijk,
std::vector< std::complex<r_8> >& almk);
//!Synthesis (spin = 0)
/*! \brief Coeffs -> Pixels with Input/Output using T floating representation
\input flmn: 3D complex spherical-laguerre coefficients
\output fijk: 3D real function values
*/
void Synthesis(const vector< complex<r_8> >& flmn, vector<r_8>& fijk) {