Hi S,
here is one way. First of all, the transform pair is
q*F(q) = Int r*F(r) sin(q*r) dr
r*F(r) = (1/(2*pi))*Int q*F(q) sin(q*r) dq
so the basic functions are qFq = q*F(q) and rFr = r*F(r). If you want F(r) itself you can always obtain it from rFr ./ r and similarly with q. The code below takes r*F(r) for the domain (0,rmax) and for mathematical convenience creates an odd (antisymmetric) function on the domain (-rmax,rmax). That means the sine terms in the fft will survive because they are odd, and the cosine terms will not because they are even. The fft produces an antisymmetric function of q in the domain (-qmax,qmax). For a final result you can just ignore the functions for r<0 and q<0.
The code starts with rFr, does a transform to find qFq. Then, to test the inverse transform there is a transform back to find rFr2 and compare it with rFr. It's the same.
% create a function for positive r
n = 2e4;
delr = .001;
r = (-n:n)*delr;
rpos = r(r>=0);
Frpos = exp(-5*rpos); % F(r), r >=0
rFrpos = rpos.*Frpos;
% create antisymmetric function with r=0 at the center
% see plot 2
rFr = [-fliplr(rFrpos(2:end)) rFrpos];
% transform to q, multiply by delr to approximate the Riemmann integral
N = length(r);
qFq = -imag(fftshift(fft(ifftshift(rFr))))*delr;
delq = 2*pi/(N*delr); % golden rule for fft
q = (-n:n)*delq;
figure(1)
plot(q,qFq)
grid on
% transform back to r, multiply by delq for Riemann integral
% factor of N because the ifft divides by N and that must be canceled out
rFr2 = (1/(2*pi))*N*imag(fftshift(ifft(ifftshift(qFq))))*delq;
figure(2)
plot(r,rFr,r,rFr2)
grid on
