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
% Read complex field data for EX
vecEX = fread(fid, nx * ny * 2, 'double');
% Validate dimensions
if numel(vecEX) ~= nx * ny * 2
fclose(fid);
error('The data size does not match expected dimensions.');
end
% Reshape and create complex representation
realEX = reshape(vecEX(1:2:end), ny, nx);
imagEX = reshape(vecEX(2:2:end), ny, nx);
EX = realEX + 1i * imagEX;
fclose(fid); % Close file
% Perform Fourier Transform
fftResult = fftshift(fft2(EX));
energyDistribution = abs(fftResult).^2;
% Calculate frequency domain kx, ky
[kx, ky] = meshgrid(linspace(-pi, pi, nx), linspace(-pi, pi, ny));
% Wave number magnitude and angle
kr = sqrt(kx.^2 + ky.^2);
theta = atan2(ky, kx);
% Angle range bins
thetaBins = linspace(0, pi/2, 90); % 89 bins
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
% Plot energy as a function of angle (polar coordinates)
figure;
polarplot(thetaBins(1:end-1), angleEnergy);
title('Energy Distribution over Angle (0 to 90 degrees)');
% Normalize energy distribution
totalEnergy = sum(angleEnergy);
angleEnergyNormalized = angleEnergy / totalEnergy;
% Plot normalized energy percentage (polar coordinates)
figure;
polarplot(thetaBins(1:end-1), angleEnergyNormalized);
title('Normalized Energy Percentage over Angle (0 to 90 degrees)');
% Example analysis for a specified angle
inputAngleDegrees = 15;
inputAngleRadians = deg2rad(inputAngleDegrees);
angleMask = (theta >= 0) & (theta < inputAngleRadians);
energyInRange = sum(energyDistribution(angleMask));
fprintf('Energy from 0 to %d degrees: %f\n', inputAngleDegrees, energyInRange);
end就是这个代码,报错显示The data size does not match expected dimensions. |
|