%% 参数设置
% 点源位移网格
row = 1;
col = 1;
delta = 2e-6; % 点源间距
delta_x = ((0 : (row-1)) - floor((row-1)/2)) * delta; % x方向点源位移
delta_y = ((0 : (col-1)) - floor((col-1)/2)) * delta;% y方向点源位移
[delta_xx, delta_yy] = meshgrid(delta_x, delta_y);
delta_xx = delta_xx(:); % 展平为一维数组
delta_yy = delta_yy(:);
lambda = 355e-9;
NA = 0.9;
k = 2*pi / lambda; %真空波束
n1 = 1;
n2 = 1;
f2 = 0.3e-3; % 待测物镜焦距
d2 = 101.*f2; % 照明透镜到待测透镜距离
d1 = d2/5; % 点源到照明透镜距离 5倍放大，均匀照明物镜光瞳面 。 也就是确定倍率后，根据光瞳面大小确定光源间距
f1 = 5/6.*d1; % 照明透镜焦距
d3 = 101/100.*f2; % 待测物镜到像面光栅距离
d4 = 3e-3; %像面光栅到ccd
d_objgrat = 2e-6; % 物光栅面采样间距 %增大可以增加 物镜在衍射面尺寸占比
r = 1e-3; % 物光栅 光阑半径
px = 200e-6; % 物面光栅周期
pi_ = 2e-6; % 像面光栅周期（MATLAB中pi为内置函数，用pi_表示）
fill_factor = 0.5; % 占空比
% 和矢量仿真有关的参数
L = 4096; % 全部采样数（需为2的幂）
z=0; % 观察平面位于焦平面（z=0）
m0=1/2; % 坐标偏移（避免原点奇异性）
xmax=0.9; % 输出图像半宽（0.9个波长）
lens_radius = d3 *NA; % 透镜半径 0.0091
%lens_radius = d3 ./ sqrt(1-(1./NA.^2)); % 透镜直径 0.0232
%
N = 128/2; %透镜半径所占的像素个数
d_obj=lens_radius/N; % 透镜面的空间分辨率（μm/像素） ≈ 142 μm/px
% 物光栅面坐标网格生成
x1 = linspace(-L*d_objgrat/2, L*d_objgrat/2 - d_objgrat, L);
y1 = x1;
[X1, Y1] = meshgrid(x1, y1);
aper = square_aper(X1, Y1, r, r); % 物面方形光阑
obj_grat_x = rect_grating(X1, px, fill_factor); % 物面光栅
phi_1 = exp(1i * (-k/(2*f1) * (X1.^2 + Y1.^2))); % 准直透镜相位
%
figure;imagesc(obj_grat_x.*aper),colorbar;title('obj_grat_x')
%
%待测透镜 出瞳面空间坐标系生成 before lens (x_inf, y_inf)
% m0=1/2偏移避免原点位于网格中心，防止奇异性。 *d_obj转化为真实数据坐标
[x_inf,y_inf]=meshgrid(linspace(-L/2+m0,L/2-m0,L)*d_obj,-linspace(-L/2+m0,L/2-m0,L)*d_obj);
[theta_inf,rho_inf] = cart2pol(x_inf,y_inf);
% zernike_coefficients = generate_custom_zernike_coefficients();
% a = calculate_zernike_phase(zernike_coefficients, lens_radius, x_inf, y_inf); % 直接在光阑区域归一化生成像差
% %
% obj_aper = circ_aper(x_inf, y_inf, lens_radius); %物镜面光瞳大小网格
% phi = exp(1i * (-k/(2*f2)*(x_inf.^2 + y_inf.^2) + k*a)); % 含像差的透镜相位; 前面e4 ；ka才0.04左右
% figure;imagesc(a),colorbar;title('a')
% 角谱计算 利用透镜面坐标参数算出 角谱
k1= 1*k; % 输出面介质波数
dk=k*NA/N; % 角谱分辨率
kx=x_inf*k1/f2; % x方向波数 ，根据论文公式4。x_inf/f 应该是透镜上某一点到焦点的角度。tan
ky=y_inf*k1/f2;
dxf = d_objgrat/100; % *********傅里叶后输出平面(像面光栅处)空间分辨率（μm/px）*********
kM1=k1*ones(L,L); % 波数 矩阵
kz1=sqrt(kM1.^2 - kx.^2 - ky.^2); % z方向波数分量
%
xmax2=round(xmax*lambda/dxf); % 像面光栅处聚焦光斑 的 物理尺寸（波长）→ 像素坐标
%相位校正
[X,Y]=meshgrid(linspace(-L/2,L/2,L),linspace(-L/2,L/2,L)); %%%校正辅助坐标
PhaseShift=m0*2*pi*X/L+m0*2*pi*Y/L; %correction phase
%软边孔径函数 生成平滑孔径，避免硬边界衍射。tanh函数的斜率因子 1.5 控制边缘过渡宽度。就是论文公式11
Theta=0.5*(1+tanh((1.5/1)*(N-sqrt((x_inf/d_obj).^2+(y_inf/d_obj).^2))));
cx=ones(L,L); cy=zeros(L,L); % x lin horizontal
% 像面光栅坐标网格 dxf
x2 = linspace(-L*dxf/2, L*dxf/2 - dxf, L);
y2 = x2;
[X2, Y2] = meshgrid(x2, y2);
img_grat_x = rect_grating(X2, pi_, fill_factor); % 位移后的像面光栅
aper2 = square_aper(X2, Y2, r/100, r/100); % 物面方形光阑
% 矢量角谱
[FX_img, FY_img] = meshgrid(linspace(-1/(2*dxf), 1/(2*dxf) - 1/(L*dxf), L), linspace(-1/(2*dxf), 1/(2*dxf) - 1/(L*dxf), L));
% 像面光栅处的波数分量
kx_img = 2 * pi * FX_img;
ky_img = 2 * pi * FY_img;
kz_img = sqrt(kM1.^2 - kx_img.^2 - ky_img.^2);
obliquity = kz_img/k; % 矢量角谱传播的倾斜因子
% CCD 面的空间坐标网格
[X_ccd, Y_ccd] = meshgrid(linspace(-L*10*dxf/2, L*10*dxf/2 - 10*dxf, L), linspace(-L*10*dxf/2, L*10*dxf/2 - 10*dxf, L));
%
delta_g = [0.5].*pi_; %pi_像面光栅周期
% 主循环：遍历光栅位移
n_dg = numel(delta_g);
for n = 1:n_dg
current_time = datetime('now', 'Format', 'HH:mm:ss.SSS');
% 显示迭代次数和时间 +
fprintf('迭代: %d / %d | 当前时间: %s\n', n, n_dg, current_time);
I_x = zeros(size(aper));
I_y = zeros(size(aper));
I_z = zeros(size(aper));
E_field_x = zeros(size(aper));
E_field_y = zeros(size(aper));
E_field_z = zeros(size(aper));
I = zeros(size(aper));
img_grat_x = rect_grating(X2, pi_, fill_factor);
if n == 1
I_tong = zeros(size(aper), 'single'); % 物瞳面光强
I_grat = zeros(size(aper), 'single'); % 光栅面光强
end
for s = 1:numel(delta_xx) %将二维空间转化为一维序列，遍历不同的光源点
% 添加进度显示：每2000次迭代或到达终点时显示
if mod(s, 1000) == 0 || s == 10201
current_time_inner = datetime('now', 'Format', 'HH:mm:ss.SSS');
fprintf(' 点源迭代: %d / %d | 当前时间: %s | 剩余: %d 次\n', ...
s, 10201, current_time_inner, 10201 - s);
end
dx_batch = delta_xx(s);
dy_batch = delta_yy(s);
% 计算点源到输入面的距离
R = sqrt((X1 + dx_batch).^2 + (Y1 + dy_batch).^2 + d1^2); % 需要改成相减嘛
E0 = (1./R) .* exp(1i*k*R); % 球面波
E0_x = E0 .* aper .* obj_grat_x .* phi_1; % 调制光阑、光栅和透镜相位
[E1_x, dx1] = fresnel_batch(E0_x, d_objgrat,lambda , d2); % E1_x 就是物镜前表面入射场
if n == 1 % 仅n=1时计算
I_tong = I_tong + abs(E1_x.*Theta).^2;
end
E_incx = E1_x;
%input x 论文里面公式5：假设入射场为x偏振光，远场输出的3个偏振分量。就是3*1矩阵的每一行代表一个分量
%x component
E_infxx=sqrt(n1/n2)*sqrt(kz1./kM1).*Theta.*E_incx.*(ky.^2+kx.^2.*kz1./kM1)./(kx.^2+ky.^2);
%y component
E_infxy=sqrt(n1/n2)*sqrt(kz1./kM1).*Theta.*E_incx.*(-kx.*ky+kx.*ky.*kz1./kM1)./(kx.^2+ky.^2);
% %z component
E_infxz=sqrt(n1/n2)*sqrt(kz1./kM1).*Theta.*E_incx.*(-(kx.^2+ky.^2).*kx./kM1)./(kx.^2+ky.^2);
%y input
E_incy = E_incx;
%x component
E_infyx=sqrt(n1/n2)*sqrt(kz1./kM1).*Theta.*E_incy.*(-ky.*kx+kx.*ky.*kz1./kM1)./(kx.^2+ky.^2);
E_infyy=sqrt(n1/n2)*sqrt(kz1./kM1).*Theta.*E_incy.*(kx.^2+ky.^2.*kz1./kM1)./(kx.^2+ky.^2);
E_infyz=sqrt(n1/n2)*sqrt(kz1./kM1).*Theta.*E_incy.*(-(kx.^2+ky.^2).*ky./kM1)./(kx.^2+ky.^2);
%factors to assemble E_inf 公式10前面的系数
CF2=-1i*f2*exp(-1i*f2*2*pi*n2/lambda)*exp(1i*kz1.*z)./(2*pi*kz1);
%%E_inf 傅里叶的线性性质，FFT(af)=aFFT(f)；这里CF2就是a
Fieldx=CF2.*(cx.*E_infxx +cy.*E_infyx);
Fieldy=CF2.*(cx.*E_infxy +cy.*E_infyy);
Fieldz=CF2.*(cx.*E_infxz +cy.*E_infyz);
Ex0=ifftshift(ifft2(fftshift((Fieldx)))).*exp(1i*PhaseShift);
Ey0=ifftshift(ifft2((fftshift(Fieldy)))).*exp(1i*PhaseShift);
Ez0=ifftshift(ifft2((fftshift(Fieldz)))).*exp(1i*PhaseShift);
% 2. 矢量场通过像面光栅调制 ******再加上像面光阑
Ex1 = Ex0 .* img_grat_x.* aper2; % x偏振分量与x光栅调制
Ey1 = Ey0 .* img_grat_x.* aper2; % y偏振分量与x光栅调制（假设光栅对x/y偏振影响相同）
Ez1 = Ez0 .* img_grat_x.* aper2; % z偏振分量与x光栅调制
% figure();
% imagesc(abs(Ex1));
% 空间频率域中的传播核
H_prop = exp(1i * kz_img * d4);
% 使用角谱法传播每个偏振分量
% 1. 对像面上的场量进行 FFT
FFT_Ex_after_grating = fftshift(fft2(Ex1));
FFT_Ey_after_grating = fftshift(fft2(Ey1));
FFT_Ez_after_grating = fftshift(fft2(Ez1));
% 2. 乘以传播核
FFT_Ex_at_ccd = FFT_Ex_after_grating .* H_prop;
FFT_Ey_at_ccd = FFT_Ey_after_grating .* H_prop;
FFT_Ez_at_ccd = FFT_Ez_after_grating .* H_prop;
% 3. IFFT 得到 CCD 平面上的场量
Ex_at_ccd = ifft2(ifftshift(FFT_Ex_at_ccd));
Ey_at_ccd = ifft2(ifftshift(FFT_Ey_at_ccd));
Ez_at_ccd = ifft2(ifftshift(FFT_Ez_at_ccd));
I_x = I_x + abs(Ex_at_ccd).^2; % ccd累加光强
I_y = I_y + abs(Ey_at_ccd).^2;
I_z = I_z + abs(Ez_at_ccd).^2;
end
end
I = I_x + I_y + I_z ; 分析总结我的代码
最新发布
本文展示了一个使用C语言编程的实例，通过迭代计算的方法来近似求解指数函数e^x的值。代码中利用了数学库中的fabs函数来判断迭代误差，确保计算精度达到1e-8级别。
扫一扫
2839
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
