%parameters
r                   = 2.5e-3;  %radius of pupil
z0                  = 50e-3;   %distance from pupil to image
rho_0_vec           = linspace(0,50e-6,500);  %vector of observation-plane radii
wavelength          = 500e-9;
k                   = 2*pi/wavelength;

U_0                 = [];

bessel_arg_vec      = 2*pi*rho_0_vec*r/wavelength/z0;
W                   = @(rho) 1*rho.^2;  %aberration as a function of rho, the normalized pupil radius

disp('Working on calculation')
for ii  = 1:length(rho_0_vec)

    rho_0           = rho_0_vec(ii);
    if rho_0 == 0
        rho_0 = 1e-10;
    end
    bessel_arg      = bessel_arg_vec(ii);
    mybessel        = @(rhoprime) rhoprime.*exp(j*2*pi*W(rhoprime)).*besselj(0,bessel_arg*rhoprime);
    field_value     = quad(mybessel,0,1);
    U_0             = [U_0 field_value];

end

I_0                 = abs(U_0).^2;

figure(100);plot(rho_0_vec*1e6,I_0)

disp('Calculation complete')