内容简介：推导见正文，代码见代码，另外一个很有用的Shack-Hartmann波前探测器资源链接见其他资源
逼逼叨：做Shack-Hartmann波前探测器，需要波前重构，懒得自己写，想找个现成的，发现有限差分最小二乘居然要卖一百多巨款！！！matlab有自带的函数（ lsqr - 求解线性方程组 - 最小二乘法），矩阵自己手动写一下，发现要不了多久，放出来。
推导
目的：需要根据测量得到的波前斜率重构波前相位。 方程：写在word里懒得敲了，直接截图了
结果如下：
代码
%测量单个二次相位透镜，看不同的波前斜率的输出的位置
close all;
%% 参数设置
M_len=55;N_len=55;%单个透镜包含的像素数量
num_len_x=15;num_len_y=15;%x,y方向包含的透镜数量
M=num_len_y*M_len;
N=num_len_x*N_len;
x=(-(N-1)/2:(N-1)/2).*dx;
y=(-(M-1)/2:(M-1)/2).*dy;
[Y,X]=ndgrid(x,y);%z=0面的坐标
%% 生成相位分布，计算斜率
phase_test_set=-k0/2*(X.*X+Y.*Y)./(r_max*num_len_y*sqrt(2));
x_center=zeros(num_len_y,num_len_x);%中心位置
y_center=zeros(num_len_y,num_len_x);
%中心位置坐标
for i=1:num_len_x
for j=1:num_len_y
x_center(i,j)=(-(num_len_x+1)/2+j)*N_len*dx;
y_center(i,j)=(-(num_len_y+1)/2+i)*M_len*dy;
end
end
S_x_set=-k0*x_center./(r_max*num_len_y*sqrt(2));%中心斜率
S_y_set=-k0*y_center./(r_max*num_len_y*sqrt(2));
r2_center=x_center.*x_center+y_center.*y_center;
phi_center_set=-k0/2*(r2_center)./(r_max*num_len_y*sqrt(2));
%% 根据dW/dx,DW/dy波前重构
%认为恢复的相位的坐标在透镜的中心位置采用southwell
Px=N_len*dx;Py=M_len*dy;
%构造矩阵
%Sx和系数矩阵的上半部分
SA_x=zeros(num_len_y*(num_len_x-1),1);
B_h=zeros(num_len_y*(num_len_x-1),2*num_len_y*num_len_x);
B_l=zeros((num_len_y-1)*num_len_x,2*num_len_y*num_len_x);
for i=1:num_len_y
for j=1:num_len_x-1
SA_x((i-1)*(num_len_x-1)+j)=S_x_set(i,j)+S_x_set(i,j+1);
%B的上半部分
hang=(i-1)*(num_len_x-1)+j;%方程编号
lie_1=(i-1)*(num_len_x)+j;%对应的phi位置
lie_2=lie_1+1;
B_h(hang,lie_1)=-1;
B_h(hang,lie_2)=1;
end
end
%Sy和系数矩阵的下半部分
SA_y=zeros(num_len_x*(num_len_y-1),1);
for i=1:num_len_y-1
for j=1:num_len_x
SA_y((i-1)*num_len_x+j)=S_y_set(i,j)+S_y_set(i+1,j);
%B的下半部分
hang=(i-1)*(num_len_x)+j;%方程编号
lie_1=(i-1)*(num_len_x)+j;%对应的phi位置
lie_2=i*(num_len_x)+j;
B_l(hang,lie_1)=-1;
B_l(hang,lie_2)=1;
end
end
%组合SA矩阵
SA=zeros(2*num_len_y*(num_len_x-1),1);
SA(1:num_len_y*(num_len_x-1),1)=Px/2*SA_x;
SA(1+num_len_y*(num_len_x-1):2*num_len_y*(num_len_x-1),1)=Py/2*SA_y;
B=zeros((num_len_y-1)*num_len_x+num_len_y*(num_len_x-1),2*num_len_y*num_len_x);
B(1:num_len_y*(num_len_x-1),:)=B_h;
B(1+num_len_y*(num_len_x-1):2*num_len_y*(num_len_x-1),:)=B_l;
phi_matrix = lsqr(B,SA,1e-10,500);
phi_solve=zeros(num_len_y,num_len_x);
for i=1:num_len_y
for j=1:num_len_x
phi_solve(i,j)=phi_matrix((i-1)*(num_len_x)+j,1);
end
end
phi_solve=phi_solve-mean(mean(phi_solve))+mean(mean(phi_center_set));%减去相差的常数
figure;
meshc(x,y,phase_test_set);%原相位
hold on;
plot3(x_center,y_center,phi_solve,"*");%恢复的相位
其他资源
matlab社区上的沙克-哈特曼模拟器 基本每一步都有，要考虑的都考虑到了，而且都写成函数了，直接调用就行，非常好用。 引用格式：SergioB (2023). Shack-Hartmann Simulator (https://github.com/SergioBonaqueGonzalez/Shack-Hartmann-Simulator), GitHub. 检索来源 2023/9/12. 资源链接：https://ww2.mathworks.cn/matlabcentral/fileexchange/63851-shack-hartmann-simulator
淘宝直通车怎么开最省钱？直通车需要注意哪些事项？
拳击手多大适龄 一般从多大年龄开始练拳击较合适？