%ccly_Patterns_Synthesis_lowSLL
%%MATLAB控制CST的思想基础：CST的建模仿真的命令都是VB命令，而要想使用MATLAB来控制CST就需要将应用之间的通信接口作为一个桥梁架接在MATLAB与CST之间，COM组件就是做这样的事情。
%%由于CST程序中的帮助文件关于VB语言的使用非常详细，所以我们需要做的就是将这座桥梁搭好。注：理论上所有手动操作CST的VB命令能在帮助文件里面找到对应代码
%%将VB命令转化为在MATLAB中能够运行的代码是MATLAB自动控制CST进行建模仿真的关键
%%典型的VB命令转化为CST代码如下所示[1]:
% With Brick                                    brick = invoke(mws, 'Brick')
%     .Reset                                    invoke(brick, 'Reset');
%     .Name ("brick1")                          invoke(brick, 'Name', "brick1");
%     .Component ("component1")                 invoke(brick, 'Component', "component1");
%     .Material ("PEC")               ==》      invoke(brick, 'Material', "PEC");
%     .Xrange (0, 2)                            invoke(brick, 'Xrange', 0, 2);
%     .Yrange (0, 3)                            invoke(brick, 'Yrange', 0, 3);
%     .Zrange (0, "a+3")                        invoke(brick, 'Zrange', '0', num2str(a+3));
%     .Create                                   invoke(brick, 'Create');
% End With                                      release(brick);
%基本思路就是找出mws中需要操作的对象，然后对其所有参数进行赋值或者执行命令，然后释放对象

%%参考文献：
%%[1]CST. "Application Note - Calling CST Studio from Matlab"
clc;
clear;
close all;

cst = actxserver('CSTStudio.application');%首先载入CST应用控件
mws = invoke(cst, 'NewMWS');%新建一个MWS项目
app = invoke(mws, 'GetApplicationName');%获取当前应用名称
ver = invoke(mws, 'GetApplicationVersion');%获取当前应用版本号
fullfilename=which('Patch_Array_1x10.cst');%打开一个已经有仿真结果的CST文件,此CST文件要和本m文件在同一文件目录下
invoke(mws, 'OpenFile', fullfilename);%打开一个存在的CST文件

%%模型文件初始化结束
%参数声明
e0=8.85418782e-12;
mu0=1.25663706e-6;
c=1/sqrt(e0*mu0);%光速常量，m/s
f=10*1e9;%10GHz
lambda=c/f;
M=10;%阵元个数
d=lambda/2;%阵元间隔，m   

%一些会用到的对象
combineresults = invoke(mws, 'CombineResults');%方向图合成操作
farfieldplot = invoke(mws, 'FarfieldPlot');
export = invoke(mws, 'ASCIIExport');
plot1D = invoke(mws, 'Plot1D');

%%本程序综合设计一个低副瓣的方向图
%最便捷的设计方法就是采用现成的窗函数来生成幅度激励，如采用Taylor分布
nbar = 5; sll = -20;
taper=taylorwin(M,nbar,sll).';%泰勒窗函数
amp=taper/max(taper);%归一化幅度值，不用归一化也可以的
pha=zeros(1,10);%同相激励

%有源方向图合成，将生成的幅度和相位代入计算，然后比较方向图是否满足副瓣要求。
lablename= 'Patterns_SLL_20';
invoke(combineresults, 'Reset');
invoke(combineresults, 'SetMonitorType', 'frequency');
invoke(combineresults, 'SetOffsetType', 'phase');
invoke(combineresults, 'EnableAutomaticLabeling', 'False');
invoke(combineresults, 'SetLabel', lablename);
invoke(combineresults, 'SetNone');
invoke(combineresults, 'AddFilter' ,':farfield:');
invoke(combineresults, 'OnlyAtFrequency' ,'10');%只对10GHz的频点进行分析，也可以设置成对所有的远场观测点进行综合，参照help文档里面VBA部分关于combineresults的说明即可
for i=1:M
    invoke(combineresults, 'SetExcitationValues' ,'port', i, 1, amp(i), pha(i));%给每个端口设置激励
end
invoke(combineresults, 'Run')%计算合成的方向图
%显示并导出贴片阵列的沿组阵方向的方向图
Select_name = ['Farfields\farfield (f=10) [', lablename, ']'];
invoke(mws, 'SelectTreeItem', Select_name);
invoke(farfieldplot, 'Reset');
invoke(farfieldplot, 'Plottype', 'cartesian');
invoke(farfieldplot, 'Vary', 'angle1');
invoke(farfieldplot, 'Phi', '0');
invoke(farfieldplot, 'Step', '1');
invoke(farfieldplot, 'SetTheta360', 'True');
invoke(farfieldplot, 'SymmetricRange', 'True');
invoke(farfieldplot, 'StoreSettings');
invoke(farfieldplot, 'SetPlotMode', 'realized gain');
invoke(farfieldplot, 'Plot');
File_name = ['Gain_', lablename, '.txt'];
%在当前文件夹内新建一个txt文档用来保存方向图
fid=fopen(File_name,'w+');
fclose(fid);
%获取方向图的文件路径，方便下面保存结果时用
File_path=which(File_name);
invoke(export, 'Reset');
invoke(export, 'FileName', File_path);
invoke(export, 'Execute');
A=importdata(File_path);
Pattern_SLL=A.data(:,3);
pks=findpeaks(Pattern_SLL,'NPeaks',2,'SortStr','descend');%从大到小寻找到两个极值点
synthesized_sll=pks(2)-pks(1)%求出副瓣电平
%至此，如果得到的副瓣电平满足设计要求，则终止；要是不满足，则可以调整taylorwin生成的窗使得SLL更小一点

release(combineresults);
release(mws);
release(cst);
