clc
clear
lambda = 1.064; %um
n1 = 1;
n2 = 2.585;
d_nom = 590; %um
NA = 0.65;
alpha = asin(NA);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%定义参数

Phi_sa = @(rho) -2 * pi * d_nom * (sqrt(n2 ^ 2 - (NA * rho) .^ 2) - sqrt(n1 ^ 2 - (NA * rho) .^ 2)) / lambda;
Dn2 = @(rho) 2 * pi * d_nom * sqrt(n2 ^ 2 - (NA * rho) .^ 2) / lambda;

phi_ave = integral(@(rho) rho .* Phi_sa(rho), 0, 1) / integral(@(rho) rho, 0, 1);
dn2_ave = integral(@(rho) rho .* Dn2(rho), 0, 1) / integral(@(rho) rho, 0, 1);

phi_sa_p = @(rho) Phi_sa(rho) - phi_ave;
dn2_p = @(rho) Dn2(rho) - dn2_ave;

s = 1 / (1 + integral(@(rho) phi_sa_p(rho) .* dn2_p(rho) .* rho, 0, 1) / integral(@(rho) dn2_p(rho) .* dn2_p(rho) .* rho, 0, 1));
R = 1 / s;
Phi_SA_hat = @(rho) 2 * pi * d_nom * (s * sqrt(n1 ^ 2 - (NA * rho) .^ 2) - sqrt(n2 ^ 2 - (NA * rho) .^ 2)) ./ (lambda * s);
Phi_SA_hat = @(rho) Phi_SA_hat(rho) - Phi_SA_hat(0) + mod(Phi_SA_hat(0), 2 * pi) - 2 * pi;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%不同rho条件下相位补偿值（列出文献公式）

pixelx = 1280;
pixely = 1024;
F = 1;
pixel = F * pixely;
%%%%%%%%%%%%%%%%%%%%%%%设置相位图像素数目，设置调制区域线性占比F(填充因子) slm分辨率 1280*1024，

[x, y] = meshgrid(1:pixelx, 1:pixely);
x = x - (pixelx + 1) / 2; %坐标移到中心
y = y - (pixely + 1) / 2;
D = (pixel - 1) / 2; % +-都可以 长度刚好等于上面矩阵边界的大小
Aperture = (x .^ 2 + y .^ 2) .* 4 ./ pixel .^ 2; % 画一个环形的有效区域
Aperture(Aperture <= 1) = 1;
Aperture(Aperture > 1) = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 设置约束孔径Aperture
rho = (sqrt(x .^ 2 + y .^ 2) / D) .* Aperture; % 关键； 上面对rho的定义

Phase = Phi_SA_hat(rho);
offset = -min(min(Phi_SA_hat(rho)));
Phase = Phi_SA_hat(rho) + offset; %相位数值转变为全正
Phase = mod(Phase, 2 * pi); %相位折叠
PhaseHolo = round((Phase * 255 / (2 * pi))) .* Aperture; %相位归一化到0-255

ApertureN = (x .^ 2 + y .^ 2) .* 4 ./ pixel .^ 2;
ApertureN(ApertureN <= 1) = 0;
ApertureN(ApertureN > 1) = 1;
Noise = round(255 * rand(pixely, pixelx)) .* ApertureN;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%生成0-255间随机噪声，并设置相反约束孔径ApertureN ； 增加噪声的实验

PhaseHolo = PhaseHolo + Noise; %添加随机噪声后的相位图

figure(1)
imshow(PhaseHolo, [])
figure(2)
r1 = linspace(-1, 1, pixelx);
% plot(r1,PhaseHolo(512,:))
plot(r1, PhaseHolo(512, :) * 2 * pi / 255, 'r')
axis([-1, 1, 0, 6 * pi])
%%%%%%%%%%%%%%%%%%%%%%%%%%水平直径上的相位分布
filename = ['C:\Users\simin\Desktop\球差校正_SiC加工相位图1280X1024_8位\', 'd=', num2str(d_nom), 'um', '_F=', num2str(F), '.bmp'];
imwrite(uint8(PhaseHolo), filename); % % % % % % % % % % % % % % % % % % % % % % % % % %导出8位相位图存储