lambda = 1.064;
d = 590;
n1 = 1.00;
n2 = 2.585;
NA = 0.65;
alpha = asin(NA / n1);
k0 = 2 * pi / lambda;
x = linspace(-3, 3, 100); %linspace(-6,6,200);
D_z = linspace(-50, 50, 800);
z = D_z * n2 / n1 + (n2 / n1 - 1) * d;
%v=2*pi*n1/lambda*sqrt(2*x.^2)*(NA/n1);
v = 2 * pi * n1 / lambda * abs(x) * (NA / n1);
u = 2 * pi * n2 / lambda * z * (NA / n1) ^ 2;
[u, v] = meshgrid(u, v);
%%  Fresnel transmission coefficients
Rs = @(phi1, phi2)((n1 * cos(phi1) - n2 * cos(phi2)) / (n1 * cos(phi1) + n2 * cos(phi2)));
Rp = @(phi1, phi2)((n1 * cos(phi2) - n2 * cos(phi1)) / (n1 * cos(phi2) + n2 * cos(phi1)));
Ts = @(phi1, phi2)(1 - Rs(phi1, phi2)); %可能错误
Tp = @(phi1, phi2)(1 - Rp(phi1, phi2));
%% Sphere Aberration
Phi_d = @(phi1, phi2, d)(-d * (n1 * cos(phi1) - n2 * cos(phi2)));

%% I transfored integrals
I0_in = @(phi1, phi2)cos(phi1) ^ 0.5 * sin(phi1) * (Ts(phi1, phi2) + Tp(phi1, phi2) * cos(phi2)) ...
    * exp(1i * u .* cos(phi2) / sin(alpha) ^ 2) * exp(1i * k0 * Phi_d(phi1, phi2, d)) .* besselj(0, v * sin(phi1) / sin(alpha));
I1_in = @(phi1, phi2)cos(phi1) ^ 0.5 * sin(phi1) * (Tp(phi1, phi2) * sin(phi2)) ...
    * exp(1i * u .* cos(phi2) / sin(alpha) ^ 2) * exp(1i * k0 * Phi_d(phi1, phi2, d)) .* besselj(1, v * sin(phi1) / sin(alpha));
I2_in = @(phi1, phi2)cos(phi1) ^ 0.5 * sin(phi1) * (Ts(phi1, phi2) - Tp(phi1, phi2) * cos(phi2)) ...
    * exp(1i * u .* cos(phi2) / sin(alpha) ^ 2) * exp(1i * k0 * Phi_d(phi1, phi2, d)) .* besselj(2, v * sin(phi1) / sin(alpha));

%%  integree
intnum = 500;
Dphi = alpha / intnum;
I0 = 0;
I1 = 0;
I2 = 0;

for ii = 0:intnum
    disp(ii);
    phi1 = ii * Dphi;
    phi2 = asin(sin(phi1) * n1 / n2);
    I0 = I0 + I0_in(phi1, phi2) * Dphi;
    I1 = I1 + I1_in(phi1, phi2) * Dphi;
    I2 = I2 + I2_in(phi1, phi2) * Dphi;
end

%%
l0 = 1; %l_0 is an amplitude factor
f = 250; %f����
K = pi * n1 * f * l0 / lambda;
%��xz��һ�������� theta_p=0
e_2x = -1i * K * (I2 + I0);
e_2y = 0;
e_2z = -2 * K * I1;

intensity = (e_2x .* conj(e_2x) + e_2y .* conj(e_2y) + e_2z .* conj(e_2z));

figure;
imagesc(z, x, intensity);
xlabel('z/um');
ylabel('x/um');

% [c,h]=contour(z,x,intensity/max(intensity(:)),[0.025 0.05 0:0.1:1],'-');
% clabel(c,h);

%%
rr = mat2gray(intensity(round(length(x) / 2):end, round(length(D_z) / 2)));
r = linspace(0, x(end), length(rr));
theta = linspace(0, 2 * pi, 100);
[R, Theta] = meshgrid(r, theta);
zz = meshgrid(rr, theta);
xx = R .* cos(Theta);
y = R .* sin(Theta);
figure
pcolor(xx, y, zz);
axis square;
shading interp;