function [nx, ny, dx, dy, isPol, lambda, EX, EY, PilotRays] = read_ZBF(filename)
% 接收一个文件名，返回采样点数（nx，ny），网格距离（dx，dy），偏振状态（ispol），
% 波长（lambda），电场数据（Ex，Ey），引导光线的数据（PilotRays），
% 以及照度和相位的图像（intensityImage, phaseImage）
% 打开文件
fid = fopen(filename, 'r');
% 读取第一部分（文件格式，采样点数以及偏振状态）
vec = fread(fid, 9, 'int32'); % 假设整数是32位的，根据实际情况可能需要调整
fileFormat = vec(1);
nx = vec(2);
ny = vec(3);
isPol = vec(4);
% 读取第二部分（根据文件格式）
if fileFormat == 1
vec = fread(fid, 20, 'double');
lambda = vec(9);
PilotRays.posX = vec(3);
PilotRays.posY = vec(6);
PilotRays.rayleighX = vec(4);
PilotRays.rayleighY = vec(7);
PilotRays.waistX = vec(5); % 注意：原代码中waist可能是waistX的笔误，但这里保留原样
PilotRays.waistY = vec(8); % 同上
else
vec = fread(fid, 11, 'double');
lambda = vec(5);
% 对于旧格式，假设PilotRays的某些字段是相同的（这里只是示例，可能需要根据实际情况调整）
PilotRays.posX = vec(3);
PilotRays.posY = vec(3); % 通常posY应该有不同的值，但这里只是示例
PilotRays.rayleighX = vec(4);
PilotRays.rayleighY = vec(4); % 通常rayleighY应该有不同的值，但这里只是示例
PilotRays.waistX = vec(6); % 同上，可能是waistX的笔误
PilotRays.waistY = vec(6); % 同上
end
dx = vec(1);
dy = vec(2);
% 读取电场数据EX
vecEX = fread(fid, nx * ny * 2, 'double');
realEX = reshape(vecEX(1:2:end), [ny, nx]);
imagEX = reshape(vecEX(2:2:end), [ny, nx]);
EX = realEX + 1i * imagEX;
EX = EX'; % 转置以匹配通常的图像格式（行表示y，列表示x）
% 如果光束是偏振的，读取电场数据EY
if isPol == 1
vecEY = fread(fid, nx * ny * 2, 'double');
realEY = reshape(vecEY(1:2:end), [ny, nx]);
imagEY = reshape(vecEY(2:2:end), [ny, nx]);
EY = realEY + 1i * imagEY;
EY = EY'; % 转置
else
EY = zeros(size(EX));
end
% 数据大小为 Nx x Ny
Nx = 1024; % 从文件中解析
Ny = 1024; % 从文件中解析
dataSize = Nx * Ny;
% 形成复数形式的光场
complexField = fread(fid, nx * ny * 2, 'double');
% 执行傅里叶变换
fftResult = fftshift(fft2(complexField));
% 计算傅里叶变换后，在频域中得到的维度信息
[X, Y] = meshgrid(1:Nx, 1:Ny);
% 中心点转换为波数向量
centerX = ceil(Nx / 2);
centerY = ceil(Ny / 2);
kx = X - centerX;
ky = Y - centerY;
kr = sqrt(kx.^2 + ky.^2); % 横向波数幅值
% 假定最大波数为标量单位
k0 = max(kr( );
% 计算与z轴的夹角（出射角）
theta = asin(min(kr / k0, 1)); % 确保比例在有效区间内
% 划分角度范围
thetaBins = linspace(0, pi/2, 90); % 将0到90度分为89个区间
angleEnergy = zeros(1, length(thetaBins)-1);
for i = 1:length(thetaBins)-1
mask = (theta >= thetaBins(i)) & (theta < thetaBins(i+1));
angleEnergy(i) = sum(energyDistribution(mask));
end
% 绘制能量随角度的分布（极坐标图）
figure;
polarplot(thetaBins(1:end-1), angleEnergy);
title('Energy Distribution over Angle (0 to 90 degrees)');
% 归一化能量分布
totalEnergy = sum(angleEnergy);
angleEnergyNormalized = angleEnergy / totalEnergy;
% 绘制能量百分比分布（极坐标图）
figure;
polarplot(thetaBins(1:end-1), angleEnergyNormalized);
title('Normalized Energy Percentage over Angle (0 to 90 degrees)');
% 输入一个角度（以度为单位）
inputAngleDegrees = 15; % 示例：计算从 0 到 15 度的能量
% 将输入角度转换为弧度
inputAngleRadians = deg2rad(inputAngleDegrees);
% 创建角度掩码，筛选出从 0 到输入角度的范围
angleMask = (theta >= 0) & (theta <= inputAngleRadians);
% 计算从 0 到输入角度范围内的能量总和
energyInRange = sum(energyDistribution(angleMask));
% 计算总能量
totalEnergy = sum(energyDistribution();
% 计算范围内的能量占总能量的百分比
percentageEnergy = (energyInRange / totalEnergy) * 100;
% 输出结果
fprintf('The energy from 0 to %d degrees is %.2f.\n', inputAngleDegrees, percentageEnergy);
|
|
