stktest.m

STK的二次开发

function STKTest

delete(get(0,'children'));

% Matlab连接STK，首先要获得STK的地址，示例中采用了默认地址stkDefaultHost。
% 得到地址后，就可以利用语句stkOpen打开默认地址，建立Matlab与STK的连接了。
remMachine = stkDefaultHost;
conid=stkOpen(remMachine);

dtr = pi/180;
rtd = 180/pi;

% 在STK中创建新场景前，需要检查STK中的当前场景，如果目前STK中已经存在一个场景，
% 就可以利用语句stkUnload关闭场景，或者利用语句stkClose关闭当前连接。
scen_open = stkValidScen;
if scen_open == 1
   rtn = questdlg('Close the current scenario?');
   if ~strcmp(rtn,'Yes')
      stkClose(conid)
      return
   else
      stkUnload('/*')
   end
end

% 利用语句stkNewObj创建新场景。这里需要注意一个问题，stkNewObj在STK5.0中创建新场景时，
% 语法为stkNewObj('Scenario','','Scenario22-1')，这应该是STK5.0的程序错误所致。
% stkNewObj的正确语法应该是stkNewObj('ParentPath','Class','Name'),
% 其中ParentPath为场景或对象的路径，Class为新建对象的类别，Name为新建对象的名称。
disp('Create a new scenario');
stkNewObj('Scenario','','Scenario22-1');

% 场景建立完成后，需要为场景设置时间参数。
% 其中场景的历元时刻设置这是通过语句stkSetEpoch和stkSyncEpoch实现的。
% 而设置场景运行方式时是通过语句stkConnect直接调用STK连接命令的方式实现的。
% STK调用连接命令语句stkConnect非常重要，
% 它的正确语法是rtn=stkConnect(conid,'<Command>','<ObjectPath>','<Parameters>')
% 其中，conid代表通过语句stkOpen打开的当前连接，Command代表要执行的命令，
% ObjectPath代表对象名称，Parameters代表相应命令的执行参数。
disp('Set scenario time period');
stkSetEpoch('1 May 2000 00:00:00.0','GREGUTC');
stkSyncEpoch;
rtn = stkConnect(conid,'Animate','Scenario/Scenario22-1',...
    'SetValues "1 May 2000 00:00:00.0" 60 0.1');
rtn = stkConnect(conid,'Animate','Scenario/Scenario22-1','Reset');

% 创建飞机“unknowplane”和“lanjiefeiji”。
% 创建对象时仍采用语句stkNewObj，同时利用stkConnect语句创建的飞行设置三维模型。
% 这里需要注意的语句是stkSetWaypoints，它是STK的航迹设置语句。
stkNewObj('*/','Aircraft','unknownplane');
stkSetWaypoints('Scenario/Scenario22-1/Aircraft/unknownplane','1 May 2000 00:00:00.0',...
    [26.15580056*pi/180,25.77953895*pi/180,25.59306383*pi/180,25.30582979*pi/180,...
     25.19944681*pi/180,24.99731915*pi/180,24.84838298*pi/180,24.63561702*pi/180,...
     24.46540426*pi/180,24.39093617*pi/180,24.21008511*pi/180,24.03987234*pi/180,...
     24.06114893*pi/180,23.88029787*pi/180,23.71008510*pi/180;121.54734890*pi/180,...
     120.31490293*pi/180,120.09594603*pi/180,120.03523778*pi/180,119.76205059*pi/180,...
     119.76205059*pi/180,119.49898145*pi/180,119.42815514*pi/180,119.48886341*pi/180,...
     119.22579427*pi/180,119.34721079*pi/180,119.25614840*pi/180,119.02343339*pi/180,...
     118.92225295*pi/180,119.09425969*pi/180;...
     3048,3048,3048,3048,3048,3048,3048,3048,3048,3048,3048,3048,3048,3048,3048,],...
     [100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,]);
stkConnect(conid,'VO','Scenario/Scenario22-1/Aircraft/unknownplane','Model "D:\f16.mdl"');
stkNewObj('*/','Aircraft','lanjiefeiji');
stkSetWaypoints('Scenario/Scenario22-1/Aircraft/lanjiefeiji','1 May 2000 00:30:00.0',...
    [24.98155918*pi/180,24.56978218*pi/180,24.24752192*pi/180,24.10429514*pi/180,...
     24.19381187*pi/180,24.42655540*pi/180,24.71444728*pi/180;118.10581246*pi/180,...
     118.2760905*pi/180,118.51447977*pi/180,118.87206367*pi/180,...
     119.34884219*pi/180,119.58723145*pi/180,120.04183343*pi/180;3048,...
     3048,3048,3048,3048,3048,3048],[100,100,100,100,100,100,100]);
stkConnect(conid,'VO','Scenario/Scenario22-1/Aircraft/lanjiefeiji','Model "D:\an30.mdl"');

% 创建雷达站“leidazhan1”和“leidazhan2”，并重新设置两个雷达站的位置。
% 同时为两个雷达站分别设置侦察雷达“haifang1”和“haifang2”。
% 这里需要注意的语句是stkSetFacPosLLA，它专门用于设置地面站的位置。
disp('Create a leidazhan1');
stkNewObj('*/','Facility','leidazhan1');
stkSetFacPosLLA('Scenario/Scenario22-1/Facility/leidazhan1',...
    [26.841*pi/180;119.802*pi/180;0]);
disp('Create a hangfangleida1 radar');
stkNewObj('Scenario/Scenario22-1/Facility/leidazhan1','Radar','haifang1');

disp('Create a leidazhan2');
stkNewObj('*/','Facility','leidazhan2');
stkSetFacPosLLA('Scenario/Scenario22-1/Facility/leidazhan2',...
    [24.742*pi/180;118.472*pi/180;0]);
disp('Create a haifangleida2 radar');
stkNewObj('Scenario/Scenario22-1/Facility/leidazhan2','Radar','haifang2');

% 创建一个导弹对象
disp('Create a new missle');
stkNewObj('Scenario/Scenario22-1','Missile','missile');

% 建立两部雷达“haifang1”和“haifang2”对飞机“unkownplane”的侦察关系。
% 这种关系的建立是通过语句stkAccess实现的。
stkAccess('Scenario/Scenario22-1/Facility/leidazhan1/Radar/haifang1','*/Aircraft/unknownplane');
stkAccess('Scenario/Scenario22-1/Facility/leidazhan2/Radar/haifang2','*/Aircraft/unknownplane');

% 创建两架飞机“unknownplane”和“lanjiefeiji”的动态数据结构.
% 这主要是通过语句stkDynState实现的。该语句能为一个机动目标产生一个动态数据结构，
% 从中可以获取动态目标的位置、速度、俯仰角等信息值。
cm =  stkDynState('*/Aircraft/unknownplane',10);
intx = stkDynState('*/Aircraft/lanjiefeiji',10);

cmMin = min(cm.times);
cmMax = max(cm.times);
cmTimes = linspace(cmMin, cmMax, 1000);

% 计算在指定时间段内，飞机“unknownplane”的飞行状态，
% 从中可以获取飞机每一时刻的飞行数据。
% 这主要是通过语句atbGeometry.m实现的，它能产生一个很大的数组用来存储相关数据。
cmCurData = atbGeometry(cmTimes, cm);

% 利用语句atbFlatten降低飞机“unknownplane”的“pos”数据的维数。
cmPos= atbFlatten(cmCurData, 'pos')';
class(cmPos)
cmPos1=atbFlatten(cmCurData,'pos')';

% 利用语句atbCbfToLla将CBF坐标转换为LLA坐标。
cmLLA = atbCbfToLla(cm.cb, cmPos);

intxMin = min(intx.times);
intxMax = max(intx.times);

intxTimes = linspace(intxMin, intxMax, 1000);
intxCurData = atbGeometry(intxTimes, intx);
intxPos= atbFlatten(intxCurData, 'pos')';
intxLLA = atbCbfToLla(intx.cb, intxPos);

% 利用drawnow语句分别画出飞机“unknownplane”和“lanjiefeiji”的飞行轨迹。
h1 = plot(intxLLA(2,:)*180/pi,intxLLA(1,:)*180/pi,'b',...
     cmLLA(2,:)*180/pi,cmLLA(1,:)*180/pi,'r');
  xlabel('Longitude')
  ylabel('Latitude')
  title('Interceptor/Target Trajectories')
set(h1,'LineWidth',1);
drawnow;
  
aimTimeU = linspace(max([cmMin intxMin]),min([cmMax intxMax]),1000);
[cmCurData,intxCurData] = atbGeometry(aimTimeU, cm, intx);

% 利用Matlab继续进行相关的数据处理。
% 并计算出两架飞机在空中的交接处。
disp(' ');
disp('Compute the local horizontal-referenced elevation angle ...');
disp(' ');

vvlhRelPos=zeros(3,1000);
relAER = zeros(3,1000);

for k=1:1000
   t = atbCbfToVvlhMtx(cmCurData(k).pos, cmCurData(k).vel);
   vvlhRelPos(:,k) = t*cmCurData(k).relPos;
   relAER(1,k) = atan2(-vvlhRelPos(1,k), vvlhRelPos(2,k));
   relAER(3,k) = norm(vvlhRelPos(:,k));
   relAER(2,k) = acos(vvlhRelPos(3,k)/relAER(3,k)) - pi/2;
end   

ind = find(relAER(2,:) > 0);

aimTime = linspace(min(aimTimeU(ind)),max(aimTimeU(ind)),1000);
[cmCurData,intxCurData] = atbGeometry(aimTime, cm, intx);

vvlhRelPos=zeros(3,1000);
relAER = zeros(3,1000);
for k=1:1000
   t = atbCbfToVvlhMtx(cmCurData(k).pos, cmCurData(k).vel);
   vvlhRelPos(:,k) = t*cmCurData(k).relPos;
   relAER(1,k) = atan2(-vvlhRelPos(1,k), vvlhRelPos(2,k));
   relAER(3,k) = norm(vvlhRelPos(:,k));
   relAER(2,k) = acos(vvlhRelPos(3,k)/relAER(3,k)) - pi/2;
end   


intxPosIV = atbFlatten(intxCurData, 'pos')';
intxVelIV = atbFlatten(intxCurData, 'vel')';
intxRelVelIV = atbFlatten(intxCurData, 'relVel')';
cmPosIV = atbFlatten(cmCurData, 'pos')';
cmVelIV = atbFlatten(cmCurData, 'vel')';

intxLLAIV = atbCbfToLla(intx.cb, intxPosIV);
cmLLAIV = atbCbfToLla(cm.cb, cmPosIV);

% 在Matlab图形窗口中，用粗线条表示两架飞机在空中的交接过程。
figure(1);
hold on;
h2 = plot(intxLLAIV(2,:)*180/pi,intxLLAIV(1,:)*180/pi,'r',...
     cmLLAIV(2,:)*180/pi,cmLLAIV(1,:)*180/pi,'b');
  
set(h2,'LineWidth',2);
drawnow;

% 计算由拦截飞机“lanjiefeiji”所发射导弹的轨道和速度。
disp(' ');
disp('Compute a long range air intercept missile trajectory...');
disp(' ');

lInd = 10;

dT = aimTime(2) - aimTime(1);  
curPos = intxPosIV(:,lInd);
curVel = intxVelIV(:,lInd);
speed = norm(curVel);

longAccel = 20; 
maxSpeed = 500;

curPosSave = zeros(3,1000);
curVelSave = zeros(3,1000);
rangeSave = zeros(1,1000);
timeSave = zeros(1,1000);

rDot = -norm(intxRelVelIV(:,lInd));
rDotW = rDot;
prevRange = relAER(3,lInd);
k = lInd;
j = 0;

while (k < 1000)
    k = k+1;
    j = j+1;
    
    curSpeed = norm(curVel);
    curVelNorm = curVel/curSpeed;
    speed = min([speed+longAccel*dT maxSpeed]);
    
    curVel = curVelNorm * speed;	
	curPos  = curPos + curVel*dT;
 	cmPos = cmPosIV(:,k);
	
    relPos = cmPos - curPos;
	range = norm(relPos);

	cmCurVel = cmVelIV(:,k);
	cmSpeed = norm(cmCurVel);

	psi = acos(-relPos'*cmCurVel/cmSpeed/range);
	lead = asin(cmSpeed / curSpeed * sin(psi));

	xUnit = relPos / norm(relPos);
	yUnit = cross(cross(xUnit, cmCurVel/cmSpeed), xUnit);

	curVel = (xUnit*cos(lead) + yUnit*sin(lead)) * speed;

    curPosSave(:,j) = curPos;
	curVelSave(:,j) = curVel;
	rangeSave(j) = range;
    timeSave(j) = aimTime(k);
   
	rDotW = rDotW*0.9 + rDot;
	rDot = (range - prevRange)/dT;
	prevRange = range;
	
	if (range < dT*speed)
		break
	end
end

curPosSave(:,j+1:1000) = [];
curVelSave(:,j+1:1000) = [];
rangeSave(:,j+1:1000) = [];
timeSave(:,j+1:1000) = [];

aimLLA = atbCbfToLla(cm.cb, curPosSave);

% 在Matlab图形窗口中，用粉红色线条来表示导弹的轨迹。
h3 = plot(aimLLA(2,:)*180/pi,aimLLA(1,:)*180/pi,'m');

set(h3,'LineWidth',1);
drawnow;

% 将导弹的轨道数据以文件的形式保存起来。
% 这主要是通过语句stkSetEphemerisCBF实现的，该语句专门用来运载器类目标的位置数据。
disp(' ');
disp('Create an external state file and tell STK to load it.')
disp(' ');

tmpDir = getenv('TEMP');
if isempty(tmpDir)
    tmpDir = pwd;
end
if (strcmp(computer,'PCWIN'))
	ephem_path = ['D:\ephem.e'];
else
	ephem_path = ['D:\ephem.e'];
end

stkSetEphemerisCBF('*/Missile/missile',cm.cb,timeSave,curPosSave,curVelSave,ephem_path);

fid=fopen('curPosSave.txt','w');
fprintf(fid,'%6.2f %12.8f\n',curPosSave);
fclose(fid);

fid=fopen('cmPos1.txt','w');
fprintf(fid,'%6.2f %12.8f\n',cmPos1);
fclose(fid);

% 利用stkAccReport语句产生两部雷达对飞机“unknownplane”的截获报告。
[data1,names1]=stkAccReport('Scenario/Scenario22-1/Facility/leidazhan2/Radar/haifang2',...
    '*/Aircraft/unknownplane','AER');
[data2,names2]=stkAccReport('Scenario/Scenario22-1/Facility/leidazhan1/Radar/haifang1',...
    '*/Aircraft/unknownplane','AER');

% 利用stkFindData语句从报告中导出需要的数据文件。
time1 = stkFindData(data1{1},'Time');
time2 = stkFindData(data2{1},'Time');

Range1 = stkFindData(data1{1},'Range');
Range2 = stkFindData(data2{1},'Range');

% 利用Matlab绘制出两部雷达对飞机“unknownplane”的跟踪过程示意图。
% 其中用红色代表雷达“haifang1”的跟踪效果，用蓝色代表雷达“haifang2”的跟踪效果。
figure(2);
subplot(1,2,1);
plot(time1,Range1,'r');

title('Plot haifang1 radar to unknownplane Range')
xlabel('Time (s)')
ylabel('Range (m)')
subplot(1,2,2);
plot(time2,Range2,'b');
title('Plot haifang2 radar to unknownplane Range')
xlabel('Time (s)')
ylabel('Range (m)')

% 利用stkConnect语句调整STK三维显示窗口的设置。
stkConnect(conid, 'VO * 3dView Eye FromTo', 'Aircraft/unknownplane', ' ');
stkConnect(conid, 'Zoom',...
    'Scenario/Scenario22-1 Object Scenario/Scenario22-1/Facility/leidazhan2', '1.5');
       
% 结束Matlab与STK的连接关系，并关闭此次连接。
stkClose(conid)
stkClose