clc
clear
lambda=1.064;%单位um
d=204;%穿透深度
d_p=600;%实际加工深度
n1=1;
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(-100,100,800);
% z=D_z*2.88+(2.88-1)*d;
z=linspace(-d,600,1000);
v=2*pi*n1/lambda*abs(x)*(NA/n1);
u=2*pi*n2/lambda*z*(NA/n1)^2;
[u,v]=meshgrid(u,v);
 f=4000; %f焦距um
%%  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));
Ts=@(phi1,phi2)2*n1*cos(phi1)/(n1*cos(phi1)+n2*cos(phi2));
Tp=@(phi1,phi2)2*n1*cos(phi1)/(n2*cos(phi1)+n1*cos(phi2));
% Ts=@(phi1,phi2)0.5;
% Tp=@(phi1,phi2)0.5;
%% Sphere Aberration
Phi0=@(phi1,d)phi1+asin(cos(phi1)/f*(d_p*n1*sin(phi1)/sqrt(n2^2-n1^2*sin(phi1)^2)-d*tan(phi1)));
Phi_d=@(phi1,phi2,d)(-d*(n1*cos(phi1)-n2*cos(phi2)))-(n1*(f*cos(Phi0(phi1,d))-d)/cos(phi1)+n2*d_p/cos(phi2))+n1*(f-d)+n2*d_p;
% Phi_d=@(phi1,phi2,d)d*n2*cos(phi2)*(1-1/0.3472); %Salter算法校正
% Phi_d=@(phi1,phi2,d)-1.8904*d*n2*cos(phi2); %中科院郑炜算法校正
% Phi_d=@(phi1,phi2,d)0;

% I Transformed Integrals
I0_in = @(phi1,phi2) integral(@(phi1) I0_integrand(phi1, phi2), -pi/2, pi/2);
I1_in = @(phi1,phi2) integral(@(phi1) I1_integrand(phi1, phi2), -pi/2, pi/2);
I2_in = @(phi1,phi2) integral(@(phi1) I2_integrand(phi1, phi2), -pi/2, pi/2);

% 辅助函数定义
I0_integrand = @(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_integrand = @(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_integrand = @(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));

I0 = I0_in(phi1_value, phi2_value);
I1 = I1_in(phi1_value, phi2_value);
I2 = I2_in(phi1_value, phi2_value);

 l0=1;%l_0 is an amplitude factor
 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+d,x,intensity);
  xlabel('z/um');
  ylabel('x/um');
% axis([500,750,-3,3])
figure(2)
plot(z+d,intensity(51,:))
% [c,h]=contour(z,x,intensity/max(intensity(:)),[0.025 0.05 0:0.1:1],'-');
% [c,h]=contour(z,x,intensity/max(intensity(:)),[1.1e-1, 5.6e-2, 2.9e-2, 1.5e-2, 8.1e-3, 4.3e-3, 2.2e-3, 1.2e-3, 6.2e-4, 3.3e-4, 1.7e-4, 9.1e-5, 4.8e-5, 2.5e-5, 1.3e-5, 6.9e-6, 3.7e-6, 1.9e-6, 1e-6, 5.3e-7]/1.3e-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;