%ccly_Patterns_Synthesis_Fan
%%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"

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

%变量
global elementpattern lambda D theta M %这里设置全局变量是方便目标函数使用
f=10*1e9;%10GHz
lambda=c/f;%波长
M=10;%阵元个数
d=lambda/2;%阵元间隔，单位：m    
D=0:d:(M-1)*d;%各单元到第一个单元的距离

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

%%本程序为生成一个扇形波束，采用粒子群优化算法进行波束形状的综合设计

%在不考虑有限大阵列的边缘单元效应和精确互耦效应的情况下，对于小规模的阵列，我们可以利用阵列中间的一个单元的有源方向图进行综合设计，这样节省前期的计算时间
%在优化过程中，假如每次都通过Matlab来操控CST利用combine results来合成方向图，计算时间会比较长。
%因此当面对的是一个均匀线阵时，直接采用方向图乘积原理进行方向图的分析可以很快完成优化设计

%显示阵列中间的一个贴片单元的有源方向图并导出，我们取第5个单元的方向图数据
elementnun=num2str(round(M/2));%取阵列中间的一个单元方向图，这里是第5个
Select_name = ['Farfields\farfield (f=10) [', elementnun, ']'];%选中这个方向图
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, 'SetScaleLinear', 'True');%显示的方向图是线性值，方便计算合成方向图
invoke(farfieldplot, 'StoreSettings');
invoke(farfieldplot, 'SetPlotMode', 'realized gain');
invoke(farfieldplot, 'Plot');
File_name = ['elementpattern_', elementnun, '.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);
elementpattern=A.data(:,3);%第三列是我们想要的，打开此文件看头就明白了
theta=A.data(:,1);%CST保存全向方向图时角度只到是从-180到179，刚好就360个点

%这里采用Matlab自带的粒子群算法，对于不是特别复杂的问题，省心又省力
nvars = M;%粒子群优化的变量个数，由于对称性，十个单元只需综合其中5个单元的幅度和相位
lb = zeros(nvars,1);%变量下限都为0
ub = ones(nvars,1);%变量上限都为1  
options = optimoptions('particleswarm','SwarmSize',200,'HybridFcn',@fmincon,'MaxStallIterations',100);  %可以参考help设置
rng(2018) %随机数的种子，纪念2018年
%调用优化算法进行计算了，目标函数PSO_fun在下文，此函数优化到满足终止条件为止并输出最优解。
opt = particleswarm(@PSO_fun,nvars,lb,ub,options) ; %PSO_fun函数另起了一个m文件
amp = [opt(1:M/2) fliplr(opt(1:M/2))];%用最优解得到幅度
pha = [180*double(opt(M/2+1:M)>0.5) fliplr(180*double(opt(M/2+1:M)>0.5))];%用最优解得到相位，这里的相位单位是°，方便写入CST
%利用最优解合成方向图
lablename= 'Patterns_Fan';
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, 'SetScaleLinear', 'False');%显示的方向图是线性值，方便计算合成方向图
invoke(farfieldplot, 'StoreSettings');
invoke(farfieldplot, 'SetPlotMode', 'realized gain');
invoke(farfieldplot, 'Plot');
%优化完成后CST会显示扇形方向图，假如不满足要求可以自己设置微调一下激励的幅度

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



