Fourier-Bessel Transform

Contents

$$G(\rho) = 2\pi \int_0^\infty g(r)J_0(2\pi\rho r) r dr$$

use quad

$$g(r)=1/r$$

$$H(r,\rho) = (1/r) J_0(2\pi\rho r) r$$

rho = 0.1:0.1:2;

% note that \rho in the denominator cancels \rho in the numerator
hankel1 = @(r,rho)  besselj(0,2*pi*rho*r);

g = zeros(size(rho));
for k=1:length(rho)
    % use quad (adaptive Simpson's rule)
    % results are noisy and very sensitive to second limit
    g(k) = 2*pi*quad(hankel1,0,20,[],[],rho(k));
end

u = linspace(0.1,2,101);
v = 1./u;
plot(u,v,rho,g,'ko');
xlabel('radial freq \rho');
ylabel('value');

use quadl

g = zeros(size(rho));
for k=1:length(rho)
    % use quadl (recursive adaptive Lobatto quadrature)
    g(k) = 2*pi*quadl(hankel1,0,50,[],[],rho(k));
end

u = linspace(0.1,2,101);
v = 1./u;
plot(u,v,rho,g,'ko');
xlabel('radial freq \rho');
ylabel('value');

gaussian example

$$g(r)=\mbox{exp}(-\pi r^2)$$

rho = linspace(0,2,21);
hankel2 = @(r,rho) exp(-pi*r.^2).*besselj(0,2*pi*rho*r).*r;

g = zeros(size(rho));
for k=1:length(rho)
    % try second limit equal 50, does not work
    g(k) = 2*pi*quad(hankel2,0,5,[],[],rho(k));
end

u = linspace(0,2,101);
v = exp(-pi*u.^2);
plot(u,v,rho,g,'ko');
xlabel('radial freq \rho');
ylabel('value');

exponential example

$$g(r)=\frac{e^{-r}}{r}$$

rho = linspace(0,2,21);
hankel3 = @(r,rho) exp(-r).*besselj(0,2*pi*rho*r);

g = zeros(size(rho));
for k=1:length(rho)
    g(k) = 2*pi*quad(hankel3,0,5,[],[],rho(k));
end

u = linspace(0,2,101);
v = 2*pi./sqrt(4*(pi*u).^2+1);
plot(u,v,rho,g,'ko');
xlabel('radial freq \rho');
ylabel('value');

demonstrate nested function

use previous example

g = zeros(size(rho));
for k=1:length(rho)
g(k) = quadp(hankel3,0,5,rho(k));
end

u = linspace(0,2,101);
v = 2*pi./sqrt(4*(pi*u).^2+1);
plot(u,v,rho,g,'ko');
xlabel('radial freq \rho');
ylabel('value');


Published with MATLAB® 7.13