test.m

Fast Hankel Transform Algorithm

%
% This script tests the Hankel transform code using the examples from
% Anderson (1979).
%
%
% First, generate the vector of B values.
%
bvec=0.05:0.05:5.0;
bvec=[0.0001 0.001 0.01 0.02 0.03 0.04 bvec];
truevec=zeros(size(bvec));
ourvec=zeros(size(bvec));
diffs=zeros(size(bvec));

%
%  Test 1, int(g*exp(-g*g)*J0(g*b),g=0..infinity)=exp(-b^2/4)/2.
%
figure(1)
for i=1:length(bvec),
  B=bvec(i);
  ourvec(i)=hankel0('c1',B);
  truevec(i)=exp(-B^2/4)/2;
end;
subplot(2,1,1);
plot(bvec,truevec,'k');
hold on
plot(bvec,ourvec,'k:');
xlabel('B');
ylabel('T(B)');
subplot(2,1,2);
semilogy(bvec,abs(ourvec-truevec)./abs(truevec),'k');
xlabel('B');
ylabel('Relative Error');
print -deps c1.eps

%
%  Test 2, int(g*g*exp(-g*g)*J1(g*b),g=0..infinity)=b*exp(-b*b/4)/4
%
figure(2);
for i=1:length(bvec),
  B=bvec(i);
  ourvec(i)=hankel1('c2',B);
  truevec(i)=B*exp(-B*B/4)/4;
end;
subplot(2,1,1);
plot(bvec,truevec,'k');
hold on
plot(bvec,ourvec,'k:');
xlabel('B');
ylabel('T(B)');
subplot(2,1,2);
semilogy(bvec,abs(ourvec-truevec)./abs(truevec),'k');
xlabel('B');
ylabel('Relative Error');
print -deps c2.eps


%
%  Test 3, int(exp(-g)*J1(g*B),g=0..infinity)=(sqrt(1+B^2)-1)/(B*sqrt(1+B^2))
%
figure(3);
for i=1:length(bvec),
  B=bvec(i);
  ourvec(i)=hankel1('c3',B);
  truevec(i)=(sqrt(1+B^2)-1)/(B*sqrt(1+B^2));
end;
subplot(2,1,1);
plot(bvec,truevec,'k');
hold on
plot(bvec,ourvec,'k:');
xlabel('B');
ylabel('T(B)');
subplot(2,1,2);
semilogy(bvec,abs(ourvec-truevec)./abs(truevec),'k');
xlabel('B');
ylabel('Relative Error');
print -deps c3.eps



%
%  Test 4, int(g*exp(-2*g)*J1(g*B),g=0..infinity)=
%
figure(4);
for i=1:length(bvec),
  B=bvec(i);
  ourvec(i)=hankel1('c4',B);
  truevec(i)=B/(4+B^2)^(3/2);
end;
subplot(2,1,1);
plot(bvec,truevec,'k');
hold on
plot(bvec,ourvec,'k:');
xlabel('B');
ylabel('T(B)');
subplot(2,1,2);
semilogy(bvec,abs(ourvec-truevec)./abs(truevec),'k');
xlabel('B');
ylabel('Relative Error');
print -deps c4.eps

%
%  Test 5, int(exp(-2*g)*J0(g*B),g=0..infinity)=1/(sqrt(4+B^2))
%
figure(5);
for i=1:length(bvec),
  B=bvec(i);
  ourvec(i)=hankel0('c5',B);
  truevec(i)=1/(sqrt(4+B^2));
end;
subplot(2,1,1);
plot(bvec,truevec,'k');
hold on
plot(bvec,ourvec,'k:');
xlabel('B');
ylabel('T(B)');
subplot(2,1,2);
semilogy(bvec,abs(ourvec-truevec)./abs(truevec),'k');
xlabel('B');
ylabel('Relative Error');
print -deps c5.eps