%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 功能说明：MATLAB与STK互联
% 作者：
% 版本：v1.0
% 时间：2023/7/10
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 1.初始化需要的接口函数
stkInit;                                % 建立连接
remMachine = stkDefaultHost;
conid = stkOpen(remMachine);            % 得到连接句柄（用于发送指令）
% 判断场景是否存在
scen_open = stkValidScen;
if scen_open == 1
    rtn = questdlg('Close the current scenario?');
    if ~strcmp(rtn,'Yes')
        stkClose(conid);
        return;
    else
        stkUnload('/*');
    end
end

%% 2.建立场景、航天器和观测站
% 建立场景
disp('创建场景中......');
disp('-------------------------');
stkNewObj('/','Scenario','TestScenario');

% 设定场景时间
disp('设定场景时间......');
disp('-------------------------');
% 设定场景开始时间
% startDate = '27'; startMonth = 'Nov'; startYear = '2019';
% startHour = '00'; startMin = '00';     startSec = '00.000';
% strStartTime = strcat(startDate,32,startMonth,32,startYear,32,startHour,':',startMin,':',startSec);% 逗号之间的“32”代表“空格”之意
strStartTime = '27 Nov 2019 00:00:00.000';

% 设定场景结束时间
% stopDate = '27'; stopMonth = 'Nov'; stopYear = '2019';
% stopHour = '02'; stopMin = '00'; stopSec = '00.000';
% strEndTime = strcat(stopDate,32,stopMonth,32,stopYear,32,stopHour,':',stopMin,':',stopSec);% 逗号之间的“32”代表“空格”之意
strEndTime = '28 Nov 2019 02:00:00.000';
stkSetTimePeriod(strStartTime,strEndTime,'GREGUTC');
stkSetEpoch(strStartTime,'GREGUTC');
stkSyncEpoch;

% 设定场景动画开始时间
strQuteBegTime = ['"' strStartTime '"']; % 时间写入命令是得加双引号，比较刻板
rtn = stkConnect(conid,'Animate','Scenario/TestScenario',['SetValues ' strQuteBegTime ' 10 0.1']);
% 下面这个用来参考格式
% rtn = stkConnect(conid,'Animate','Scenario/TestScenario','SetValues "27 May 2019 06:14:00.000" 60 0.1');
% 设定动画时间回到起始点
rtn = stkConnect(conid,'Animate','Scenario/TestScenario','Reset');

% 建立航天器
disp('建立航天器......');
disp('-------------------------');
stkNewObj('*/','Satellite','AoXiangI');

%% 3.设定星座中各卫星的属性
%--------------------------------------------------------------%
%               注意：这里的单位都是弧度(rad)
%                     千万别用错单位！
%--------------------------------------------------------------%
disp('仿真运行中......');
disp('-------------------------');
rtn = stkConnect(conid,'Animate','Scenario/TestScenario','Start End');
d2r = pi/180;
% 卫星积分起止时间、轨道历元时间、步长
t_start=0; t_stop=2*60*60; orbitEpoch=t_start; dt=10;
tSimu = t_start+1:t_stop/dt;
% 卫星初始轨道根数
a = 6928.137*1000;  e = 0.0; i = 53.0*d2r;
w = 0*pi/180; Raan = 160*pi/180; M = 0*d2r;

% 设置卫星轨道属性
stkSetPropClassical('*/Satellite/AoXiangI','J2Perturbation','J2000',t_start,t_stop,dt,orbitEpoch,a,e,i,w,Raan,M);
% 设置卫星姿态属性(下面写的这个是错的)
Qua = [sqrt(2)/2;0;0;sqrt(2)/2];  QuaMat = repmat(Qua,1,length(tSimu));
stkSetAttitudeCBI('*/Satellite/AoXiangI','AoXiangI',tSimu,QuaMat);

tic;
global satPos satVel
satPos = zeros(t_stop/dt,3);
satVel = zeros(t_stop/dt,3);
time   = zeros(1,t_stop/dt);
k = 1;
for t = 1:dt:t_stop
    [satPos1,satVel1] = stkPosVelCBI('*/Satellite/AoXiangI',t);
    satPos(k,:) = satPos1;
    satVel(k,:) = satVel1;
    k = k+1;
end
toc;

% 保存数据至csv中
col_str = {'sat_pos_x','sat_pos_y','sat_pos_z','sat_vel_x','sat_vel_y','sat_vel_z'};
res_table = table(satPos(:,1),satPos(:,2),satPos(:,3),satVel(:,1),satVel(:,2),satVel(:,1),'VariableNames',col_str);
writetable(res_table, 'satellite_data.csv');
stkClose(conid);