# How to calculate efficiently and accurately the Fourier transform of a radial function in Fortran

As my question states, I want to calculate the Fourier transform $F(q)$ of a radial function $f(r)$ (defined on $[0,\infty)$ and which decays like an exponential $\exp(-Ar+b)$ at large $r$) as accurately as possible in Fortran. The function values come from a data file (which I can easily interpolate through cubic interpolation for example and extrapolate since the behavior at large $r$ is known).

I'm using the "physics" definition of the Fourier transform in 3D, which gives (because $f$ is radial):

$$\int d\mathbf{r}\,f(\mathbf{r})\exp(i\mathbf{q}\mathbf{r}) = 4\pi\int\limits_0^\infty r^2\frac{\sin(qr)}{qr}f(r)\,dr$$

I first tried to calculate this integral for some chosen values of $q$ by using Gauss-Legendre quadrature, by generating some 60 or 100 abscissas and weights via the NAG routine D01BCF (D01BCF link). In the case of Gauss Legendre quadrature, the problem is to choose the interval $[0,B]$ on which to integrate. While the function $f$ loses 4 to 5 orders of magnitude from $r=10$ to $r=20$ (example), the choice of $B$ as a strong influence on the result of the calculation... When I compared the result I get to a "nearly exact" calculation (made with MATLAB but with a very long computation time), I saw that in fact this was only valid for small values of $q$ (of the order of 5, when I have to deal with values as large as 150). A Gauss-Laguerre quadrature does not give any better result, probably because of the oscillatory part of the integrand.

I then tried to compute this Fourier transform for some given values of $q$ with the routine D01ASF (D01ASF link). It is a "one-dimensional quadrature, adaptive, semi-infinite interval, weight function $\cos(ωx)$ or $\sin(ωx)$", which is exactly what I need. The results are quite convincing for $q$ up to 80 or 100 if I input absolute error tolerances of 10E-5. Problems are: I would need to go at larger $q$, and the Fourier transform $F(q)$ oscillates with a magnitude of ~10E-6 at such $q$'s. Lowering the tolerance to 10E-5 already takes some time and even makes the whole thing to output some error message from the subroutine so I don't know if 10E-6 would be feasible.

I'm thus currently wondering if trying to calculate this Fourier transform with FFT wouldn't be a good idea? The problems I face are that I don't know how to calculate radial wave functions with FFT (and also that I don't even know how to use FFT properly either since the definition of the transform is not even the same (exponent sign and argument) and that I never used it before).

Would you have ideas?

EDIT: this comes from https://stackoverflow.com/q/39127340/2320757, I've been advised to put this here on stackoverflow.

EDIT: I tried by FFT (using the routine C06FAF from NAG library). It works quite well up to some large values of $q$. The problem I face is that there is always some constant normalising factor to account for. I don't get why. This normalising factor evolves with the number $N$ of points used in the mesh. It has the for of a power law: Normalising Factor $F = N^(-0.5) \exp(9.9)$ approximately (see figure where the black line is the "exact" Fourier Transform and the green, magenta, blue, red and yellow lines are the FFT calculated for different values of $N$).

• One approach might be to expand rf(r) numerically as a sum of some basis functions (with known Fourier-sine transform) and perform radial Fourier transform analytically. But before that, isn't it sufficient to use a simple 1D FFT with sufficiently fine mesh...? Aug 26 '16 at 15:50
• @roygvib I just tried FFT. What I get is quite good. I still get some problem with normalisation but I should find it quite soon I guess. Aug 29 '16 at 16:32
• @roygvib It's weird... the normalisation factor F depends on the number of points of the mesh N. It seems to follow the trend F = N^(-0.5) x exp(9.9). I added a figure to the post to show it. Aug 30 '16 at 12:04
• Is this data for numerically given f(r)? What happens for some analytical (test) function for f(r), say, f(r) = exp(-r) etc (for which an analytical result of the radial Fourier transform is known and can be compared with the FFT result). I think we have to multiply the FFT result by some multiplicative factor (which depends on the normalization of FFT routine). We can check the result by testing with some analytical function (like above). Aug 30 '16 at 12:49
• @roygvib No, it comes from some data points. I interpolate between points with a cubic spline to be able to have a sufficiently small mesh... I'll try that and come back to you. Aug 30 '16 at 12:53

I have tried the following code in Julia for $f(r) = exp(-r)$ with 1D FFT, which seems to be working somehow. So, could you compare your code with it and see if there is some difference of normalization? (Please note that fft() below makes no normalization, i.e., $1/sqrt(N)$ is not present, while the FFT routine in NAG library may be multiplying by some normalization factor.)

function test( N, rmax )

dr = rmax / N
r = linspace( 0.0, rmax - dr, N )       # r[1:N] = 0, dr, 2dr, ..., (N-1)dr

rexp = r .* exp( - r )                  # rexp[i] = r[i] * exp(-r[i]), i=1:N
integral = - imag( fft( rexp ) ) * dr   # \int_0^rmax dr sin(qr) r exp(-r)

file = open( "test.dat", "w" )

for m = 1 : N

if m < div( N, 2 )
q = (2 * pi / rmax) * (m - 1)
else
q = -(2 * pi / rmax) * (N - m + 1)
end

H = integral[ m ]                   # obtained from FFT
H_exact = 2 * q / ( q^2 + 1 )^2     # analytical solution

if 0 <= q <= 200.0
@printf( file, "%20.10e %20.10e %20.10e\n",
q, H, H_exact )
end
end

close( file )
end

test( 10^6, 50.0 ) • The normalisation factor seems to be dr*sqrt(N). I thought I had tested that... Thank you for your time :) Aug 31 '16 at 9:53

I work with these types of transforms all the time. My tool of choice is not the full FFT, but rather a discrete Fourier sine transform (DST). Not all Fourier analysis packages have one, but it looks like NAG's C06RAF will do the trick. The example code supplied by NAG should be a good starting point, but here's the outline I follow:

1. Make sure my $q$ and $r$ arrays satisfy the Nyquist condition $\Delta q \ge \frac{\pi}{r_{max}}$. (*see footnote)
2. Call the DST routine on $r f(r)$.
3. Scale the result. You will need a least a factor $\frac{\Delta r}{2q}$ and can determine the rest empirically. The $1/2$ is because your integral runs from $0$ to $\infty$, not $-\infty$ to $\infty$.

*There is factor of $2$ difference from the ordinary FFT condition. This is because your actual $r$ domain is $[0, r_{max}]$, but the DST treats the transform as though it were $[-r_{max}, r_{max}]$, so you get double the resolution in $q$ for free.