matlab产生OAM相位图

clc;
clear all;

%读取单元相位数据
filename = ‘phase.xlsx’;
serial_num = xlsread(filename, ‘A2:A9’); % 读取序列号，
phase_lib = xlsread(filename, ‘B2:B9’); % 读取相位
angleall = xlsread(filename, ‘C2:C9’); % 读取angleall
angle = xlsread(filename, ‘D2:D9’); % 读取angle
r = xlsread(filename, ‘E2:E9’); % 读取r
w = xlsread(filename, ‘F2:F9’); % 读取w
phase_lib = mod(phase_lib + 360, 360); %加360，保证相位为正

freq = 0.9; %频率
lam = 300/freq; %波长
k = 2*pi/lam; %波数

p = 110; %超表面单元周期
num = 11; %一排单元数量num+1

x = -num/2p:p:nump/2;
y = x;
[X,Y] = meshgrid(x,y); %生成坐标网格

fx=0;
fy=0;
fz=-1500; %焦点位置

l = 1; %OAM模式

phase_need = (sqrt((X-fx).^2+(Y-fy).2+fz^2)-fz)k/pi180+latan2(Y,X)/pi180+l360; %每个单元所需补偿的相位
phase_need = phase_need - (phase_need > 360). floor(phase_need/360) *360; %归一到0~360°范围内

f1 = figure;
imagesc(phase_need); %相位图绘图

%以下是生成一个单元数据信息，记录相位图每个位置需要那个单元结构
SERIAL_NUM = zeros(size(phase_need)); %生成和phase_need大小相同的零矩阵
ANGLEALL = zeros(size(phase_need));
ANGLE = zeros(size(phase_need));
R = zeros(size(phase_need));
W = zeros(size(phase_need));

%生成矩阵

for i = 1:1:(num+1)^2
SERIAL_NUM(i) = serial_num(find((abs(phase_need(i) - phase_lib) == min(abs(phase_need(i) - phase_lib)))));
%在对应位置记录所需相位和相位表最接近的那个单元序号
ANGLE(i) = angle(find((abs(phase_need(i) - phase_lib) == min(abs(phase_need(i) - phase_lib)))));
ANGLEALL(i) = angleall(find((abs(phase_need(i) - phase_lib) == min(abs(phase_need(i) - phase_lib)))));
R = r(find((abs(phase_need(i) - phase_lib) == min(abs(phase_need(i) - phase_lib)))));
W = w(find((abs(phase_need(i) - phase_lib) == min(abs(phase_need(i) - phase_lib)))));
end

f2 = figure;
imagesc(SERIAL_NUM); %绘图

ANGLE=reshape(ANGLE,[],1);
ANGLEALL=reshape(ANGLEALL,[],1);
R=reshape(R,[],1);
W=reshape(W,[],1);
%数据整形

xlswrite([‘./单元数据’ num2str(num+1) ‘x’ num2str(num+1) ‘.xlsx’],ANGLE,‘A1:A144’);
xlswrite([‘./单元数据’ num2str(num+1) ‘x’ num2str(num+1) ‘.xlsx’],ANGLEALL,‘B1:B144’);
xlswrite([‘./单元数据’ num2str(num+1) ‘x’ num2str(num+1) ‘.xlsx’],R,‘C1:C144’);
xlswrite([‘./单元数据’ num2str(num+1) ‘x’ num2str(num+1) ‘.xlsx’],W,‘D1:D144’);
%写入xlsx文档