MATLAB与zemax怎么联调?
Zemax软件的本质是一个进行光线追迹的计算器，软件中分析（Analysis）、优化（Optimization）及公差（Tolerance）等功能，均是利用追迹后的光线数据所进行的后续的数据处理。但这些内置的数据处理功能，并不能完全满足对光学系统各种各样的分析需求。鉴于这一不足，我们希望获取到追迹后的全部光线数据，并根据需要自行编写程序，对光线数据进行更灵活的处理。相比Zemax内置的宏编程语言（ZPL），Matlab编程在数据处理上的功能更为全面，而Zemax中也提供了关于Matlab的API可供调用，因此将Zemax与Matlab进行交互，就能对光线数据执行任意的数学处理，最终完成对光学系统的分析和评估。
本次我们以一个公差分析的实例操作来说明整个Zemax与Matlab的交互流程。如下图所示，此处公差分析的需求是，对像面上的光线坐标按照一定间隔进行分组和计数，然后将统计结果绘制直方图及曲线，最后求出计数曲线的半高宽FWHM（大家可以猜猜这种分析方法是用在哪种仪器设计中）。由于Zemax没有相应的操作数来对数据的进行分组统计，就需要将光线数据导入到Matlab中数据集中，进而通过Matlab的函数执行数据处理（如histcounts函数可执行分组统计）；另一方面，由于Matlab中处理后的结果不能写入到Zemax的评价函数中去，因而在调用Zemax公差分析的API时，无法将Matlab计算出的FWHM数值作为评价函数的数值，导致公差分析输出不了FWHM的结果（可参见该帖子http://www.optzmx.com/thread-1957-2-1.html，Matlab与Zemax的交互暂时实现不了这个评价函数的互传，这是个比较大的bug，但C#和C++可以实现，这里为了数据处理上的方便和编程上的低门槛，还是采用Matlab来做），所以就需要先执行公差分析的API，分别保存出现了元件偏心/倾斜等偏差的各个镜头文件（.zmx），然后逐个加载各镜头文件并执行Matlab中的数据处理，将全部结果汇总后即为公差分析的最终结论。
综上所述，本次公差分析的完整流程为：
（1）建立Matlab与Zemax间的连接；
（2）在Matlab中加载光学系统的源文件，调用Zemax公差分析相关的API，分别保存灵敏度分析和Monte Carlo分析下，已产生表面/元件偏差的各镜头文件；
（3）依次加载上述各镜头文件，调用Zemax的API进行光线追迹，将追迹后的光线数据保存到Matlab的数据工作区中；
（4）利用MATLAB的函数命令进行数据处理，给出最终的公差分析结果。
以下是各部分的详细代码及解释：
一、 首先要建立Matlab与Zemax间的连接。Zemax提供了如下两种方式：（1）交互扩展，此时仅建立了Zemax与Matlab的连接，在连接成功后，Zemax中的任意操作（如设置光瞳-视场-波长、执行光线追迹、优化、公差等）既可直接在Zemax界面中设置，也可在Matlab中通过代码来实现（2）独立应用程序，此时不需要打开Zemax软件，只需在Matlab中执行程序代码，控制Zemax在后台执行相应操作。在Zemax的软件界面中，可找到交互扩展和独立应用程序的代码块。
两种方式没有实质区别，最终的分析结果都将通过Matlab来计算和呈现。本次选择独立应用程序的代码形式。打开独立应用程序的代码块，此时代码主体是一个函数function，代码主体格式不用改，找到 % Add your custom code here...的位置，从此处开始编写自己的代码块。
function [ r ] = Customized_tolerance( args )
if ~exist('args', 'var')
args = [];
end
% Initialize the OpticStudio connection
TheApplication = InitConnection();
if isempty(TheApplication)
% failed to initialize a connection
r = [];
else
try
r = BeginApplication(TheApplication, args);
CleanupConnection(TheApplication);
catch err
CleanupConnection(TheApplication);
rethrow(err);
end
end
end
function [r] = BeginApplication(TheApplication, args)
import ZOSAPI.*;
TheSystem = TheApplication.PrimarySystem;
% Add your custom code here...
% 此处填写特定的代码
function app = InitConnection()
import System.Reflection.*;
% Find the installed version of OpticStudio.
zemaxData = winqueryreg('HKEY_CURRENT_USER', 'Software\Zemax', 'ZemaxRoot');
NetHelper = strcat(zemaxData, '\ZOS-API\Libraries\ZOSAPI_NetHelper.dll');
% Note -- uncomment the following line to use a custom NetHelper path
% NetHelper = 'C:\Users\FXGX\Documents\Zemax\ZOS-API\Libraries\ZOSAPI_NetHelper.dll';
% This is the path to OpticStudio
NET.addAssembly(NetHelper);
success = ZOSAPI_NetHelper.ZOSAPI_Initializer.Initialize();
% Note -- uncomment the following line to use a custom initialization path
% success = ZOSAPI_NetHelper.ZOSAPI_Initializer.Initialize('C:\Program Files\OpticStudio\');
if success == 1
LogMessage(strcat('Found OpticStudio at: ', char(ZOSAPI_NetHelper.ZOSAPI_Initializer.GetZemaxDirectory())));
else
app = [];
return;
end
% Now load the ZOS-API assemblies
NET.addAssembly(AssemblyName('ZOSAPI_Interfaces'));
NET.addAssembly(AssemblyName('ZOSAPI'));
% Create the initial connection class
TheConnection = ZOSAPI.ZOSAPI_Connection();
% Attempt to create a Standalone connection
% NOTE - if this fails with a message like 'Unable to load one or more of
% the requested types', it is usually caused by try to connect to a 32-bit
% version of OpticStudio from a 64-bit version of MATLAB (or vice-versa).
% This is an issue with how MATLAB interfaces with .NET, and the only
% current workaround is to use 32- or 64-bit versions of both applications.
app = TheConnection.CreateNewApplication();
if isempty(app)
HandleError('An unknown connection error occurred!');
end
if ~app.IsValidLicenseForAPI
HandleError('License check failed!');
app = [];
end
end
function LogMessage(msg)
disp(msg);
end
function HandleError(error)
ME = MException('zosapi:HandleError', error);
throw(ME);
end
function CleanupConnection(TheApplication)
% Note - this will close down the connection.
% If you want to keep the application open, you should skip this step
% and store the instance somewhere instead.
TheApplication.CloseApplication();
end
二、 加载源文件，调用公差分析API，保存各镜头文件。注意，这里将灵敏度分析和蒙特卡洛分析的镜头文件分别进行保存，其中灵敏度分析可以评估单个公差变量对评价函数造成的改变，而蒙特卡洛分析则综合评估各公差变量造成的整体影响。在保存灵敏度分析的镜头文件时，需在关注的公差变量的操作数下方插入一个SAVE操作数，这样在执行公差分析时会额外保存两个镜头文件（文件命名为TSAV_最大/最小_0003），对应着该公差变量取最大和最小值时的镜头数据；而在蒙特卡洛分析中，可设置保存的文件数，这里设为保存20个文件，分别命名为MC_T0001~MC_T0020。执行整段代码后，生成的全部镜头文件如下方列表所示。
% Add your custom code here...
% 打开光学系统的源文件，执行公差分析，分别保存灵敏度分析和Monte Carlo分析的镜头文件，最后保存和关闭源文件。
Filepath = pwd;
FileNameSeq = System.String.Concat(Filepath, '\Spectrometer.zmx');
TheSystem.LoadFile(char(FileNameSeq),false);
% 公差的灵敏度分析可以评估单个变量对偏差的影响，通过Tolerance中的Save操作数
% 可以将灵敏度分析中生成的镜头文件进行保存。
TheTDE = TheSystem.TDE;
Operand = TheTDE.AddOperand(); % 新建一个操作数列表的行
Operand.ChangeType(ZOSAPI.Editors.TDE.ToleranceOperandType.SAVE); % 添加某个操作数，此处为SAVE
Operand.GetCellAt(2).IntegerValue = 3; % 将操作数中某列设置为相应的数值，具体可参见Zemax中ZOS-API帮助文档的实例3
% Run the Tolerancing analysis
% 建立公差分析并保存Monte Carlo分析的各样本文件
tol = TheSystem.Tools.OpenTolerancing();
% Select Sensitivity mode
tol.SetupMode = ZOSAPI.Tools.Tolerancing.SetupModes.Sensitivity;
% Select Criterion and related settings
tol.Criterion = ZOSAPI.Tools.Tolerancing.Criterions.RMSSpotRadius;
tol.CriterionSampling = 3;
tol.CriterionComp = ZOSAPI.Tools.Tolerancing.CriterionComps.None;
tol.CriterionCycle = 2;
tol.CriterionField = ZOSAPI.Tools.Tolerancing.CriterionFields.UserDefined;
% Select the number or MC runs and files to save
tol.NumberOfRuns = 20;
tol.NumberToSave = 20;
% Run the Tolerancing analysis
tol.RunAndWaitForCompletion();
% TheSystem.Tools.CurrentTool=[];
tol.Close();
% 保存并关闭当前的zmx文件
TheSystem.Save();
TheSystem.Close(false);
三、依次加载各镜头文件，调用Zemax的API进行光线追迹，将光线数据保存到Matlab数据区中。其中，光线追迹需要设置光线的光瞳、视场、波长等参数；其他注意项可参见代码中的注释；
% 加载已经保存的某个Monte Carlo文件中的某个文件
Filepath = pwd;
FileNameSeq = System.String.Concat(Filepath, '\MC_T0001.zmx');
TheSystem.LoadFile(char(FileNameSeq),false);
% 进行光线追迹并导出光线数据进行分辨率计算
% Set up Batch Ray Trace
% 设置批量光线追迹的参数
raytrace = TheSystem.Tools.OpenBatchRayTrace();
nsur = TheSystem.LDE.NumberOfSurfaces;
max_rays = 30;
normUnPolData = raytrace.CreateNormUnpol((max_rays + 1) * (max_rays + 1),ZOSAPI.Tools.RayTrace.RaysType.Real,nsur);
% define batch ray trace constants
% 设置波长和视场的限制条件
hx = 0;
max_wave = TheSystem.SystemData.Wavelengths.NumberOfWavelengths;
hy_ary = [0, 0.707, 1];
% initialize x/y image plane arrays
% 初始化视场、波长和光线数的高维矩阵，用于存放追迹后的光线坐标数据
x_ary = zeros(TheSystem.SystemData.Fields.NumberOfFields, max_wave,(max_rays+1) * (max_rays+1));
y_ary = zeros(TheSystem.SystemData.Fields.NumberOfFields, max_wave,(max_rays+1) * (max_rays+1));
% determine maximum field in Y only
max_field = 0;
for i=1:TheSystem.SystemData.Fields.NumberOfFields
if TheSystem.SystemData.Fields.GetField(i).Y > max_field; max_field = TheSystem.SystemData.Fields.GetField(i).Y ; end;
end
% 按视场、波长和光瞳的位置逐条进行光线追迹
for field = 1:length(hy_ary)
for wave = 1:max_wave
% Adding Rays to Batch, varying normalised object height hy
normUnPolData.ClearData();
waveNumber=wave;
for i = 1:((max_rays + 1) * (max_rays + 1))
px = 0.22*(rand() * 2 - 1);
py = 0.22*(rand() * 2 - 1);
while (px^2 + py^2 > 1) % 若光瞳的坐标pxpy超出归一化，重新对py随机取值
py = rand() * 2 - 1;
end
normUnPolData.AddRay(waveNumber, hx, hy_ary(field), px, py, ZOSAPI.Tools.RayTrace.OPDMode.None);
end
% Run Batch Ray Trace
raytrace.RunAndWaitForCompletion();
% Read batch raytrace and display results
normUnPolData.StartReadingResults();
% ReadNextResult(success,raynumber,errcode,vigcode,x,y,z,L,M,N,L2,M2,N2,opd,intensity)
% 其中x,y,z表示目标表面上的坐标；L,M,N表示光线经过目标面折射到其后介质中的方向余弦；
% L2,M2,N2表示目标表面法线的方向余弦；OPD代表光程差；
% Intensity 表示光线的相对透射强度，包括已定义的任何光瞳或表面切趾。
[success, rayNumber, errCode, vigCode, x, y, ~, ~, ~, ~, ~, ~, ~, ~, ~] = normUnPolData.ReadNextResult();
while success
if ((errCode == 0 ) && (vigCode == 0))
x_ary(field, wave, rayNumber) = x;
y_ary(field, wave, rayNumber) = y;
end
[success, rayNumber, errCode, vigCode, x, y, ~, ~, ~, ~, ~, ~, ~, ~, ~] = normUnPolData.ReadNextResult();
end
end
end
四、利用MATLAB的函数命令进行数据处理，给出最终的公差分析结果。注意，此时可以将光线数据另外备份到excel中，方便以后随时进行调取。这里对光线梳理的处理方法也不做具体展开，各位可根据需要编写合适的代码来处理数据。
% 将光线数据导出到excel中
Excelname = 'Raytrace.xlsx';
% squeeze函数将数组降维
raydata = [squeeze(x_ary(1,1,:)),squeeze(y_ary(1,1,:))];
writematrix(raydata,Excelname);
% 采用histcounts函数可将光线坐标按照CCD每个像素点的范围进行分组统计
% 具体使用方法可参照官网资料
% https://ww2.mathworks.cn/help/matlab/ref/histcounts.html#d124e602688
% X为各光线的坐标；edges为CCD各像素的区间；N为分组后的计数值。
N = histcounts(X,edges);