新浪博客

[转载]卫星的 两行星历 及转化为轨道六根数的MATLAB程序

2012-09-02 23:31阅读：

http://blog.sina.cn/dpool/blog/u/2714001757

原文地址：卫星的 两行星历 及转化为轨道六根数的MATLAB程序
原文作者：clover

两行式轨道数据是美国的北美防空联合司令部(North American Aerospace Defense Command ，NORAD)创立的用于描述卫星位置和速度的表达式。具体格式为：
Data for each satellite consists of three lines in the following format:

AAAAAAAAAAAAAAAAAAAAAAAA
1 NNNNNU NNNNNAAA NNNNN.NNNNNNNN +.NNNNNNNN +NNNNN-N +NNNNN-N N NNNNN
2 NNNNN NNN.NNNN NNN.NNNN NNNNNNN NNN.NNNN NNN.NNNN NN.NNNNNNNNNNNNNN

第0行是一个24字符的名字，后两行分别代表的含义为：
第 1行
列
描述
01
 行号
03-07
 卫星编号
08
 保密分级 (U=非保密的)
10-11
 国际标志符(发射年份后两位数字)
12-14
 国际标志符(那一年的发射编号)
15-17
 国际标志符(那次发射的件编号)
19-20
 TLE 历时(年份后两位数)
21-32
 TLE 历时(用一个十进制小数表示的一年中的第几日和日中的小数部分)
34-43
 平均运动的一阶时间导数
45-52
 平均运动的二阶时间导数(小数点的位置已确定)
54-61
 BSTAR阻力系数（小数点位置已确定）
63
 星历表类型
65-68
 星历编号
69
 校验和(以10为模)
(对于非数字部分：字母, 空格, 句点, 正号 = 0; 负号 = 1)

第2行
列
 描述
01
 行号
03-07
 卫星编号
09-16
 轨道的交角[度数]
18-25
 升交点赤经 [度数]
27-33
 离心率 （小数）
35-42
 近地点角距[度数]
44-51
 平近点角[度数]
53-63
 平均运动 [每日绕行圈数]
64-68
 在轨圈数
69
 校验和(以10为模)
(对于非数字部分：字母, 空格, 句点, 正号 = 0; 负号 = 1)


例如，2007年6月23号升空的中星6B卫星（由法国阿尔卡特阿莱尼亚宇航公司基于SB4000卫星平台研制的广播电视卫星,装载38个转发器，中星6B卫星将主要应用于广播电视传输，可传送300套电视节目。），其星历数据为：
CHINASAT 6B
1 31800U 07031A 07186.84787415 -.00000110 00000-0 00000+0 0 15
2 31800 24.2268 102.4131 7891523 179.2824 183.7579 1.56363120 16

很多的卫星数据都可以在http://celestrak.com/查询
星历的具体的详细的说明，见：http://celestrak.com/columns/v04n03/

将两行星历拷贝到一个文件中。转换成轨道六根数的MATLAB程序如下：

% function [oe,epoch,yr,M,E,satname] = TLE2oe(fname);
% fname is a filename string for a file containing
% a two-line element set (TLE)
% oe is a 1/6 matrix containing the orbital elements
% [a e i Om om nu]
% yr is the two-digit year
% M is the mean anomaly at epoch
% E is the eccentric anomaly at epoch
% satname is the satellite name
%
% Calls Newton iteration function file EofMe.m
function [oe,epoch,yr,M,E,satname] = TLE2oe(file1.txt);
% Open the file up and scan in the elements
fid = fopen(fname, 'r');
A = fscanf(fid,'%13c%*s',1);
B = fscanf(fid,'%d%6d%*c%5d%*3c%2d%f%f%5d%*c%*d%5d%*c%*d%d%5d',[1,10]);
C = fscanf(fid,'%d%6d%f%f%f%f%f%f',[1,8]);
fclose(fid);
satname=A;
% The value of mu is for the earth
mu = 3.986004415e5;
% Calculate 2-digit year (Oh no!, look out for Y2K bug!)
yr = B(1,4);
% Calculate epoch in julian days
epoch = B(1,5);
%ndot = B(1,6);
% n2dot = B(1,7);
% Assign variables to the orbital elements
i = C(1,3)*pi/180; % inclination
Om = C(1,4)*pi/180; % Right Ascension of the Ascending Node
e = C(1,5)/1e7; % Eccentricity
om = C(1,6)*pi/180; % Argument of periapsis
M = C(1,7)*pi/180; % Mean anomaly
n = C(1,8)*2*pi/(24*3600); % Mean motion
% Calculate the semi-major axis
a = (mu/n^2)^(1/3);
% Calculate the eccentric anomaly using mean anomaly
E = EofMe(M,e,1e-10);
% Calculate true anomaly from eccentric anomaly
cosnu = (e-cos(E)) / (e*cos(E)-1);
sinnu = ((a*sqrt(1-e*e)) / (a*(1-e*cos(E))))*sin(E);
nu = atan2(sinnu,cosnu);
if (nu<0), nu=nu+2*pi; end
% Return the orbital elements in a 1x6 matrix
oe = [a e i Om om nu];


另外一个需要的程序为：
% this function solves Kepler's equation,
% computing E as a function of M and e
%
function E = EofMe(M,e,tol)
if ( nargin<3 ), tol=1e-11; end
En = M;
En1 = En - (En-e*sin(En)-M)/(1-e*cos(En));
while ( abs(En1-En) > tol )
En = En1;
En1 = En - (En-e*sin(En)-M)/(1-e*cos(En));
end;
E = En1;

我的更多文章

下载客户端阅读体验更佳

APP专享