1. 首页
  2. 文章中心
  3. matlab计数函数

matlab计数函数_Matlab调用STK用例两则

69 阅读 0 评论
我是靠谱客的博主 有魅力抽屉，最近开发中收集的这篇文章主要介绍matlab计数函数_Matlab调用STK用例两则，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

ea4073c5b9ab5e73a149edc44e849a4d.png

在上一篇“基于STK/Matlab的Starlink星座仿真分析”（https://zhuanlan.zhihu.com/p/68385977）中，我们介绍了使用Matlab调用STK进行仿真分析的原理，这里进一步给出两个示例，主要使用STK Objects Library，若想撇除STK的界面，可以换成STK X Library。

示例1：将指定双行根数转换为J2000根数

% 本函数实现双行根数转换为J2000开普勒根数的功能
% 戴正旭 2019-12-14

function transTLE
% 读入的3TLE文件目标和输出结果目录
str3TLEFilePath = 'F:TLE3le20191027_part.txt';
strOutFilePath = 'F:TLE3le20191027_J2000.txt';

% 指定转换的卫星NORAD编号
matNoradID = [41765];

% 指定创建STK场景名称
strScenarioName = 'TranslateTLE';
% 指定创建STK场景时间
strBegTime = '9 Dec 2019 23:00:00.000';
strEndTime = '10 Dec 2019 12:00:00.000';

% 处理编号
nStarCount = size(matNoradID,2);
% 读入所有3TLE根数
fidin = fopen(str3TLEFilePath);
nReadNum = 0;% 读文件计数
vecTLE = [];
matReadNum = [];
while ~feof(fidin) && nStarCount>nReadNum
    % 读取根数
    strNameLine = fgetl(fidin);
    % 避免空格
    if isempty(strNameLine)
        continue;
    end
    % 处理名字
    strNameLine = regexp(strNameLine,' ','split');
    strName = strNameLine(2);
    
    strTLELine1 = fgetl(fidin);
    strTLELine2 = fgetl(fidin);
    
    % 获取NoradID
    vecTLELine2 = regexp(strTLELine2,' ','split');
    vecTLELine2(strcmp(vecTLELine2,''))=[];
    if find(matNoradID==str2double(vecTLELine2(2)))
        % 增加计数
        nReadNum = nReadNum+1;
        
        tmpTLE.Name = [char(strName) '_' char(vecTLELine2(2))];
        tmpTLE.Code = str2double(vecTLELine2(2));
        tmpTLE.Line1 = strTLELine1;
        tmpTLE.Line2 = strTLELine2;
        
        vecTLE = [vecTLE,tmpTLE];
        matReadNum = [matReadNum,str2double(vecTLELine2(2))];
    end
end
   
fprintf('指定卫星数：%d，读入卫星数：%dn',nStarCount,nReadNum);
disp('未找到卫星列表：');
disp(setdiff(matNoradID,matReadNum));
fclose(fidin);

%% 创建STK场景
disp('打开STK程序，创建一个场景');
try
    % 获取运行中的STK实例句柄
    uiapp = actxGetRunningServer('STK11.application');
    root = uiapp.Personality2;
    checkempty = root.Children.Count;
    if checkempty == 0
        % 如果未发现场景，新建一个
        uiapp.visible = 1;
        root.NewScenario(strScenarioName);
        scenario = root.CurrentScenario;
    else
        % 如果发现打开着的场景，询问是否关闭重建
        rtn = questdlg({'Close the current scenario?',' ','(WARNING: If you have not saved your progress will be lost)'});
        if ~strcmp(rtn,'Yes')
            return
        else
            root.CurrentScenario.Unload
            uiapp.visible = 1;
            root.NewScenario(strScenarioName);
            scenario = root.CurrentScenario;
        end
    end

catch
    % STK没有启动，新建实例获取句柄
    uiapp = actxserver('STK11.application');
    root = uiapp.Personality2;
    uiapp.visible = 1;
    root.NewScenario(strScenarioName);
    scenario = root.CurrentScenario;
end

% 设定场景时间
disp('设定场景时间');
scenario.SetTimePeriod(strBegTime,strEndTime);
% 设定场景动画开始时间
scenario.Epoch = strBegTime;
% 设定动画时间回到起始点，这里用的是Connect命令
root.ExecuteCommand('Animate * Reset');

% 记录根数转换结果
fidout = fopen(strOutFilePath,'w');
fprintf(fidout,'卫星名称;轨道历元(BJT);半长轴（km）;偏心率(°);轨道倾角(°);升焦点赤经(°);近地点幅角(°);平近点角(°)n');
% 循环创建卫星
for i = 1:size(vecTLE,2)
   tmpStarTLE = vecTLE(i);
   strProSat = tmpStarTLE.Name;
   fprintf('Create a new satellite %sn',char(strProSat));
   
   tmpSat = root.CurrentScenario.Children.New('eSatellite',char(strProSat));
   strCommendTLE = ['SetState */Satellite/' char(strProSat) ' TLE "' char(tmpStarTLE.Line1) '" "' char(tmpStarTLE.Line2) '"'];
   root.ExecuteCommand(strCommendTLE);
   propagator = tmpSat.Propagator;
   % 卫星数据结束时间设为1小时以后
   propagator.StopTime = datestr(datenum(propagator.StartTime)+1/24);
   propagator.Step = 60;
   propagator.Propagate;
   
   % 获取根数转换结果
   keplerian = tmpSat.DataProviders.Item('Classical Elements').Group.Item('J2000').Exec(propagator.StartTime,propagator.StopTime,60);
   mAxis = cell2mat(keplerian.Interval.Item(cast(0,'int32')).DataSets.GetDataSetByName('Semi-major Axis').GetValues);
   mEcce = cell2mat(keplerian.Interval.Item(cast(0,'int32')).DataSets.GetDataSetByName('Eccentricity').GetValues);
   mIncl = cell2mat(keplerian.Interval.Item(cast(0,'int32')).DataSets.GetDataSetByName('Inclination').GetValues);
   mRAAN = cell2mat(keplerian.Interval.Item(cast(0,'int32')).DataSets.GetDataSetByName('RAAN').GetValues);
   mArgP = cell2mat(keplerian.Interval.Item(cast(0,'int32')).DataSets.GetDataSetByName('Arg of Perigee').GetValues);
   mMeaA = cell2mat(keplerian.Interval.Item(cast(0,'int32')).DataSets.GetDataSetByName('Mean Anomaly').GetValues);
   
   % 格式化时间
   sTime = datestr(datenum(propagator.StartTime)+8/24,'yyyy-mm-dd HH:MM:SS.FFF');
   dAxis = mAxis(1);
   dEcce = mEcce(1);
   dIncl = mIncl(1);
   dRAAN = mRAAN(1);
   dArgP = mArgP(1);
   dMeaA = mMeaA(1);
   
   % 只输出第一行即为根数转换结果
   fprintf(fidout,'%s;%s;%.4f;%.6f;%.4f;%.4f;%.4f;%.4fn',strProSat,sTime,dAxis,dEcce,dIncl,dRAAN,dArgP,dMeaA);   
end

fclose(fidout);
uiapp.Quit;
clear uiapp root
end

示例2：使用双行根数，计算指定目标相对于测站的跟踪性能

% 本函数实现从文件读入指定卫星根数，并计算每天8:00到22:00可视弧段
% 戴正旭 2019.12.14

function TrackSectionCalculate
% 读入的3TLE文件目标和输出结果目录
str3TLEFilePath = 'F:TLE3le20191027.txt';
strOutFilePath = 'F:TLE3le20191027_track.txt';

% 指定转换的卫星NORAD编号
matNoradID = [11765];

% 指定创建STK场景名称
strScenarioName = 'TrackSection';
% 指定创建STK场景时间
strBegTime = '11 Dec 2019 12:00:00.000';
strEndTime = '12 Dec 2019 12:00:00.000';
% 本地时与UTC的差异
nFiltHour = 8;% 过滤用的工作时区UTC+8
nAddHour = 8;% 生成报告时区UTC+8

% 处理编号
nStarCount = size(matNoradID,2);
% 读入所有3TLE根数
fidin = fopen(str3TLEFilePath);
nReadNum = 0;% 读文件计数
vecTLE = [];
matReadNum = [];
while ~feof(fidin) && nStarCount>nReadNum
    % 读取根数
    strNameLine = fgetl(fidin);
    % 避免空格
    if isempty(strNameLine)
        continue;
    end
    % 处理名字
    strNameLine = regexp(strNameLine,' ','split');
    strName = strNameLine(2);
    
    strTLELine1 = fgetl(fidin);
    strTLELine2 = fgetl(fidin);
    
    % 获取NoradID
    vecTLELine2 = regexp(strTLELine2,' ','split');
    vecTLELine2(strcmp(vecTLELine2,''))=[];
    if find(matNoradID==str2double(vecTLELine2(2)))
        % 增加计数
        nReadNum = nReadNum+1;
        
        tmpTLE.Name = [char(strName) '_' char(vecTLELine2(2))];
        tmpTLE.Code = str2double(vecTLELine2(2));
        tmpTLE.Line1 = strTLELine1;
        tmpTLE.Line2 = strTLELine2;
        
        vecTLE = [vecTLE,tmpTLE];
        matReadNum = [matReadNum,str2double(vecTLELine2(2))];
    end
end

fprintf('指定卫星数：%d，读入卫星数：%dn',nStarCount,nReadNum);
disp('未找到卫星列表：');
disp(setdiff(matNoradID,matReadNum));
fclose(fidin);

%% 创建STK场景
disp('打开STK程序，创建一个场景');
try
    % 获取运行中的STK实例句柄
    uiapp = actxGetRunningServer('STK11.application');
    root = uiapp.Personality2;
    checkempty = root.Children.Count;
    if checkempty == 0
        % 如果未发现场景，新建一个
        uiapp.visible = 1;
        root.NewScenario(strScenarioName);
        scenario = root.CurrentScenario;
    else
        % 如果发现打开着的场景，询问是否关闭重建
        rtn = questdlg({'Close the current scenario?',' ','(WARNING: If you have not saved your progress will be lost)'});
        if ~strcmp(rtn,'Yes')
            return
        else
            root.CurrentScenario.Unload
            uiapp.visible = 1;
            root.NewScenario(strScenarioName);
            scenario = root.CurrentScenario;
        end
    end
    
catch
    % STK没有启动，新建实例获取句柄
    uiapp = actxserver('STK11.application');
    root = uiapp.Personality2;
    uiapp.visible = 1;
    root.NewScenario(strScenarioName);
    scenario = root.CurrentScenario;
end

% 设定场景时间
disp('设定场景时间');
scenario.SetTimePeriod(strBegTime,strEndTime);
% 设定场景动画开始时间
scenario.Epoch = strBegTime;
% 设定动画时间回到起始点，这里用的是Connect命令
root.ExecuteCommand('Animate * Reset');

% 创建移动测站
ship = root.CurrentScenario.Children.New('eShip','MyShip');
ship.SetRouteType('ePropagatorGreatArc');
route = ship.Route;
route.Method = 'eDetermineVelFromTime';
route.SetAltitudeRefType('eWayPtAltRefWGS84');

waypoints = {
    {-5.0000, 120.0000, 0.0, '11 Dec 2019 12:00:00.000'},...
    {-4.0000, 121.0000, 0.0, '12 Dec 2019 12:00:00.000'},...
    {-3.0000, 122.0000, 0.0, '13 Dec 2019 12:00:00.000'},...
    };

% Remove any previous waypoints
route.Waypoints.RemoveAll();

% Insert the waypoints
for i = 1:length(waypoints)
    waypoint = route.Waypoints.Add();
    waypoint.Latitude  = waypoints{1,i}{1,1};
    waypoint.Longitude  = waypoints{1,i}{1,2};
    waypoint.Altitude = waypoints{1,i}{1,3};
    waypoint.Time  = waypoints{1,i}{1,4};    
end

route.Propagate;
% 设置传感器
sensor = ship.Children.New('eSensor','MySensor');
% 设置测站限制条件
sensor.CommonTasks.SetPatternSimpleConic(85,0.1);

% 记录跟踪性能
fidout = fopen(strOutFilePath,'w');
fprintf(fidout,'目标;开始时间;结束时间;跟踪时长;起始方位角;最高仰角;最小距离n');
% 循环创建卫星
for i = 1:size(vecTLE,2)
    tmpStarTLE = vecTLE(i);
    strProSat = tmpStarTLE.Name;
    fprintf('Create a new satellite %sn',char(strProSat));
    
    tmpSat = root.CurrentScenario.Children.New('eSatellite',char(strProSat));
    strCommendTLE = ['SetState */Satellite/' char(strProSat) ' TLE "' char(tmpStarTLE.Line1) '" "' char(tmpStarTLE.Line2) '"'];
    root.ExecuteCommand(strCommendTLE);
    propagator = tmpSat.Propagator;
    % 卫星数据结束时间设为1小时以后
    propagator.StartTime = scenario.StartTime;
    propagator.StopTime = scenario.StopTime;
    propagator.Step = 60;
    propagator.Propagate;
    
    % 计算并获取access数据
    access = sensor.GetAccessToObject(tmpSat);
    access.ComputeAccess();
    
    % 获取跟踪AER数据
    accessDP = access.DataProviders.Item('Access Data').Exec(scenario.StartTime, scenario.StopTime);
    % 如果没有可视弧段则继续下一个目标
    if accessDP.DataSets.count<1
        continue;
    end
    
    % 获取数据
    accessStartTimes = (accessDP.DataSets.GetDataSetByName('Start Time').GetValues);
    accessStopTimes = (accessDP.DataSets.GetDataSetByName('Stop Time').GetValues);
    accessDuration = cell2mat(accessDP.DataSets.GetDataSetByName('Duration').GetValues);
    
    accessAER = access.DataProviders.Item('AER Data').Group.Item('Default').Exec(scenario.StartTime, scenario.StopTime, 10);
    RzMin = min(cell2mat(accessAER.Interval.Item(cast(0,'int32')).DataSets.GetDataSetByName('Range').GetValues));
    ElMax = max(cell2mat(accessAER.Interval.Item(cast(0,'int32')).DataSets.GetDataSetByName('Elevation').GetValues));
    matAl = cell2mat(accessAER.Interval.Item(cast(0,'int32')).DataSets.GetDataSetByName('Azimuth').GetValues);
    AlStart = matAl(1);
    % 记录文件，时间段每天8:00到22:00
    % 目标 开始时间  结束时间  跟踪时长  最高仰角  最小距离
    dStartTime = str2double(datestr(datenum(accessStartTimes(1,:))+nFiltHour/24,'HH'));
    dStopTime = str2double(datestr(datenum(accessStopTimes(1,:))+nFiltHour/24,'HH'));
    
    if(ElMax>10 && ((dStartTime>=8 && dStarTime<=22) || (dStopTime>=8 && dStopTime<=22)))
        sStartTime = datestr(datenum(accessStartTimes(1,:))+nAddHour/24,'yyyy-mm-dd HH:MM:SS.FFF');
        sStopTime = datestr(datenum(accessStartTimes(1,:))+nAddHour/24,'yyyy-mm-dd HH:MM:SS.FFF');
        % 获取航向数据
        matShipHeading = ship.DataProviders.Item('Heading').Group.Item('Fixed').Exec(char(accessStartTimes(1,:)),char(accessStopTimes(1,:)),10);
        matShipAzi = cell2mat(matShipHeading.Interval.Item(cast(0,'int32')).DataSets.GetDataSetByName('Azimuth').GetValues);
        ShipAlStart = matShipAzi(1);
        AlStart = AlStart + ShipAlStart;
        AlStart = mod(AlStart,360);
        
        fprintf(fidout,'%s;%s;%s;%.0f;%.3f;%.3f;%.0fn',strProSat,sStartTime,sStopTime,accessDuration(1),AlStart,ElMax,RzMin);
    end
    
    for i = 1:1:accessAER.Interval.Count-1
        RzMin = min(cell2mat(accessAER.Interval.Item(cast(i,'int32')).DataSets.GetDataSetByName('Range').GetValues));
        ElMax = max(cell2mat(accessAER.Interval.Item(cast(i,'int32')).DataSets.GetDataSetByName('Elevation').GetValues));
        matAl = cell2mat(accessAER.Interval.Item(cast(i,'int32')).DataSets.GetDataSetByName('Azimuth').GetValues);
        AlStart = matAl(1);
        % 目标 开始时间  结束时间  跟踪时长  最高仰角  最小距离
        dStartTime = str2double(datestr(datenum(accessStartTimes(i+1,:))+nFiltHour/24,'HH'));
        dStopTime = str2double(datestr(datenum(accessStopTimes(i+1,:))+nFiltHour/24,'HH'));
        % 获取航向数据
        matShipHeading = ship.DataProviders.Item('Heading').Group.Item('Fixed').Exec(char(accessStartTimes(i+1,:)),char(accessStopTimes(i+1,:)),10);
        matShipAzi = cell2mat(matShipHeading.Interval.Item(cast(0,'int32')).DataSets.GetDataSetByName('Azimuth').GetValues);
        ShipAlStart = matShipAzi(1);
        AlStart = AlStart + ShipAlStart;
        AlStart = mod(AlStart,360);
        
        if(ElMax>10 && ((dStartTime>=8 && dStartTime<=22) || (dStopTime>=8 && dStopTime<=22)))
            sStartTime = datestr(datenum(accessStartTimes(i+1,:))+nAddHour/24,'yyyy-mm-dd HH:MM:SS.FFF');
            sStopTime = datestr(datenum(accessStopTimes(i+1,:))+nAddHour/24,'yyyy-mm-dd HH:MM:SS.FFF');
            fprintf(fidout,'%s;%s;%s;%.0f;%.3f;%.3f;%.0fn',strProSat,sStartTime,sStopTime,accessDuration(i+1),AlStart,ElMax,RzMin);
        end
    end
end

fclose(fidout);
uiapp.Quit;
clear uiapp root
end

最后

以上就是有魅力抽屉为你收集整理的matlab计数函数_Matlab调用STK用例两则的全部内容，希望文章能够帮你解决matlab计数函数_Matlab调用STK用例两则所遇到的程序开发问题。

如果觉得靠谱客网站的内容还不错，欢迎将靠谱客网站推荐给程序员好友。

本图文内容来源于网友提供，作为学习参考使用，或来自网络收集整理，版权属于原作者所有。
点赞(46)

相关文章

matlab设置脚本,MATLAB脚本和功能
matlab设置脚本,MATLAB脚本和功能
matlab中的常见的统计函数的使用与范例介绍
matlab中的常见的统计函数的使用与范例介绍
MATLAB实现各种概率密度函数（概率密度/分布/逆概率分布函数）
MATLAB实现各种概率密度函数（概率密度/分布/逆概率分布函数）
matlab求各个点的函数值,matlab常用函数集锦
matlab求各个点的函数值,matlab常用函数集锦
matlab计数函数_Matlab调用STK用例两则
matlab计数函数_Matlab调用STK用例两则
matlab中的count函数,Excel 中COUNT函数的使用详解，详情介绍
matlab中的count函数,Excel 中COUNT函数的使用详解，详情介绍
Matlab 7数理统计
Matlab 7数理统计
Matlab中sum函数的使用
Matlab中sum函数的使用

评论列表共有 0 条评论

立即
投稿
返回
顶部