close all; clear; clc;
%%
x = linspace(-1,1,1e3);
[X,Y] = meshgrid(x);

R = sqrt(X.^2+Y.^2);
Airy_spot = @(x) (2*besselj(1,x)./x).^2;

r = linspace(0,1,1e3); phi = linspace(0,2*pi,1e3);
[R,PHI] = meshgrid(r,phi);
A = Airy_spot(R*10);

% Plot
figure(1); clf;
h = surfl(R.*cos(PHI), R.*sin(PHI), A,'light');
h(1).EdgeColor = 'none';
cm = lines; alpha = 0.3;
h(2).Color = alpha*cm(3,:) + (1-alpha)*[1 1 1];
h(2).Color = h(2).Color.*0.5;
h(1).DiffuseStrength = 1;
h(2).Position = [.8 -.8 0.35];
h(2).Style = 'local';
shading interp;
axis equal;
view(0,35);
axis off;