clear all;
close all;

R = 2.5/2; % m
lambda = 0.5e-6; % m

theta = 0.25*linspace(-1,1,1e3); % arcsecond
theta_d = theta/3600; % degree
theta_r = pi/180*theta_d; % radian

k = 2*pi/lambda;
I = abs(besselj(1,k*R*theta_r)./(k*R*theta_r)).^2;
I = I/max(I);
plot(theta,I);
grid on;
hold on;
xlim([min(theta) max(theta)]);
xlabel('Angle (arcseconds)');
ylabel('Intensity (a.u.)');

% Bessel-weighted aperture
r = linspace(0,10,1e3); % Integration grid
dr = r(2)-r(1);
R = 1.25; % m
f = besselj(0,2.4048/R*r).*(r<R);

rho = theta_r/lambda;
F = zeros(size(rho));
for i = 1:length(rho)
    F(i) = 2*pi*sum(r.*f.*besselj(0,2*pi*r*rho(i)))*dr;
end
I2 = abs(F/max(F)).^2;
plot(theta,I2,'r--');

legend('Hard-edge circular aperture',...
       'J0-weighted circular aperture',...
       'Location','NorthWest');