Fourier-Bessel Transform
Contents
use quad
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
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
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');