Discretization

In order to perform the convolution using FFT, we discretize \(s_n = n \Delta s\) with \(|n| \le N\), and require that both \(F\) and \(G\) go to zero outside of the sampled interval (large \(s_N = N \Delta s\)) and are sampled sufficiently finely (small \(\Delta s\)).

We then pick \(s_N\) so that \(G(s_N) \simeq \epsilon G(0)\), i.e., \[s_N = -s_0\,\log\epsilon + \delta\; ,\] where \(\epsilon\) is a numerical precision parameter that specifies the fractional level at which we truncate \(G(s)\) and we choose \(\delta > 0\) so that \(s_N\) lands at the next zero of the oscillation (in the limit of small \(\epsilon\)) \[e^{\delta} = \lceil Y\rceil/Y \quad ,\quad e^{\delta'} = \frac{\lceil 1/8 +Y'/4\rceil - 1/8}{Y'} \; ,\] where \(\lceil Y\rceil\) is the smallest integer \(\ge Y\) with \[Y \equiv \frac{k_0 r_0}{2\pi} \epsilon^{-s_0} \; .\] We meet the second condition by requiring \(n_s\) samples¹ of the large-\(s\) oscillation at \(s_N\). We accomplish this using \(\Delta s \le \Delta s_{\text{max}}\) with \[\Delta s_{\text{max}} = h\, \epsilon^{s_0} \quad , \quad h \equiv \frac{2\pi}{n_s k_0 r_0} \; . \label{eqn:smax}\] In order to preserve the node at \(s_N\) and ensure that \(N\) is integral, we pick \[N_G = \lceil s_N/\Delta s_{\text{max}} \rceil \quad , \quad \Delta s = s_N/N_G \; . \label{eqn:nds}\] Fig. \ref{fig:Ggrid} shows examples of \(G(s)\) tabulated with these prescriptions. Note that the symmetrized \(G(s)\) peaks at small positive \(s\) but that the oscillations for \(s > 0\) lead to a negative first moment.

We now turn to the discretization of \(F(s)\), which depends on an unknown \(f\). Since \(s = \kappa\) when evaluating \(F(s)\), our choice of \(\Delta s\) above directly determines the sampling of \(\tilde{f_{\ell m}}(k)\). To ensure that \(f\) is sufficiently sampled, we take \[\Delta s_{\text{max}} = \min(\frac{\log 10}{N_{10}}, c \epsilon^{s_0})\] where \(N_{10}\) is the minimum allowed sampling per decade of \(\tilde{f_{\ell m}}(k)\).

In practice, we want to know \(f_{\ell m}(r)\) over some range \(r_{\text{min}} < r < r_{\text{max}}\). Therefore, we set \[r_0 = \sqrt{r_{\text{min}} r_{\text{max}}}\] and calculate the corresponding \(k_0\) using the value of \(k_0 r_0\) obtained in eqn. (\ref{eqn:symm}) or eqn. (\ref{eqn:symmH}). We want a result that is free of aliasing artifacts over \(2 N_F\) samples where \[2 N_F = \log(r_{\text{max}}/r_{\text{min}})/\Delta s \; ,\] so we pad the tabulated \(G(s_n)\) with \(2 N_F\) zeros and perform transforms of length \(2(N_F + N_G)\). The resulting range of \(k\) used in the calculation has limits \[k_{\text{min}} = k_0 \exp(-(N_F + N_G)\Delta s) \quad , \quad k_{\text{max}} = k_0 \exp(+(N_F + N_G)\Delta s) \; .\]