% program to plot the difference between two first order Bessel
% functions -- this gives the IFT of the circularly symmetric bandpass 
% filter
% run from c:\jl\I8_imageprocessing\imageprocessing3\imageproc3mat\imageproc3demo

function runbessel

c1=(pi/5)*[3.0,3.0,3.0,3.0];
c2 =(pi/5)*[2.99,2.5,2.0,0.0];
for kk=1:4
for i=1:101
 i1 = -2*pi + (i-1)*(4*pi/100); 
 for j=1:101
   j1 = -2*pi + (j-1)*(4*pi/100); 
   k = sqrt(i1^2+j1^2);
   h(i,j) = (1/(2*pi*k))*(c1(kk)*besselj(1,c1(kk)*k)-c2(kk)*besselj(1,c2(kk)*k));
 end
end
AXIS([0 150 0 150 -0.1  0.1]) 
colormap(gray(256))
subplot(2,2,kk);
mesh(h)
AXIS([0 150 0 150 -0.1  0.1]) 
if(kk==1)
 title('Upper limit=3.0*(pi/5),  lower limit=2.99*(pi/5)');
end
if(kk==2)
 title('Upper limit=3.0*(pi/5),  lower limit=2.5*(pi/5)');
end
if(kk==3)
 title('Upper limit=3.0*(pi/5),  lower limit=2.0*(pi/5)');
end
if(kk==4)
 title('Upper limit=3.0*(pi/5),  lower limit=0.2');
end
rotate3d on;

end

