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

(JEC) 6/4/16 add new way of computing the Jnlp for Laguerre->Bessel transform

parent 621b71d4
......@@ -89,7 +89,11 @@ $$\begin{eqnarray}
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
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)
......@@ -97,6 +101,11 @@ f(r,\Omega) = \frac{\sqrt{2 \pi}}{R_{max}^3} \sum_{lm} \left( \sum_{p=0}^{P-1} f
\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}$.
#### 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.
......
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