16参数卫星轨道位置计算（基于Matlab）

一、代码实现

clear,clc;

fid=fopen('d:/Users/HP/Desktop/cccc/xyz.txt','r');

for i = 1:25
    tline = fgetl(fid);
end

%按行读取星历文件的每一个元素
data=textscan(fid,'%f %f %f %f %f %f %f %f %f %f\n','HeaderLines',1);
B1=cell2mat(data);
data=textscan(fid,'%f %f %f %f\n','HeaderLines',1);
B2=cell2mat(data);
data=textscan(fid,'%f %f %f %f\n','HeaderLines',1);
B3=cell2mat(data);
data=textscan(fid,'%f %f %f %f\n','HeaderLines',1);
B4=cell2mat(data);
data=textscan(fid,'%f %f %f %f\n','HeaderLines',1);
B5=cell2mat(data);
data=textscan(fid,'%f %f %f %f\n','HeaderLines',1);
B6=cell2mat(data);
data=textscan(fid,'%f %f %f %f\n','HeaderLines',1);
B7=cell2mat(data);
data=textscan(fid,'%f\n','HeaderLines',1);
B8=cell2mat(data);

%定义常量GM,Ve
GM=3.986004e14;
We=7.2921151467e-5;

%求出轨道长半轴a
sqrta=B3(4);
a=sqrta^2;
e=B3(2);

%求出平均角速度n0
n0=sqrt(GM/(a^3));



XYZ=[];

%计算从2004年1月30日16：00开始，每隔1分钟的卫星坐标
%计算从2023年6月1日0：00开始，每隔1分钟的卫星坐标
for i=0:1:20
%     计算当前时刻的儒略日
%     JulianDay=fix(365.25*(2023-1))+fix(30.6001*(6+13))+1+(0+i/60+0/3600)/24+1720981.5;
%     %JulianDay=fix(365.25*(2023+4716))+fix(30.6001*(6+1))+1+0-1524.5;
%     %将儒略日转换成BDS时
% 	  t2=(JulianDay-2444244.5)/7-fix((JulianDay-2444244.5)/7);
%     t=t2*604800-14;
% 
%     %计算卫星钟差改正后的当前时刻的GPS时
% 	t=t-(B1(8)+B1(9)*(t-B4(1))+B1(10)*(t-B4(1))^2);
% 	
%     %计算参考时刻的儒略日
% 	JulianDayoe=fix(365.25*(2023-1))+fix(30.6001*(6+13))+1+(0+0/60+0/3600)/24+1720981.5;
% 
%      %将儒略日转换成GPS时
%    	t2oe=(JulianDayoe-2444244.5)/7-fix((JulianDayoe-2444244.5)/7);
%     toe=t2oe*604800;
% 
%     %计算时间差tk
%     tk=t-toe;
    tk=-14.00;
    %tk=0;
    %计算改正平角速度
    n=n0+B2(3);

    %计算平近点角Ek
    Mk=B2(4)+n*tk;

    %迭代计算，求出偏近点角Ek
    Ek=Mk;
    Ek0=inf;

    while abs(Ek0-Ek)>1e-100
        Ek0=Ek;
        Ek=Mk+e*sin(Ek0);
    end

    %计算真近点角
	Sinfk=(sqrt(1-e^2)*sin(Ek))/(1-e*cos(Ek));
    Cosfk=(cos(Ek)-e)/(1-e*cos(Ek));
	fk=atan2(Sinfk,Cosfk);
    %计算升交角距
    phik=fk+B5(3);
    %从卫星星历文件中读取各个方向的振幅Crs,Cuc,Cus,Cic,Cis,Crc
    Crs=B2(2);
    Cuc=B3(1);
    Cus=B3(3);
    Cic=B4(2);
    Cis=B4(4);
    Crc=B5(2);

    %计算卫星轨道摄动项改正数
    deltauk=Cus*sin(2*phik)+Cuc*cos(2*phik);
	deltark=Crs*sin(2*phik)+Crc*cos(2*phik);
	deltaik=Cis*sin(2*phik)+Cic*cos(2*phik);

    %计算改正后的升交角距uk
    uk=phik+deltauk;

    %计算改正后的向径rk
	rk=a*(1-e*cos(Ek))+deltark;

    %计算改正后的倾角ik
    ik=B5(1)+deltaik+B6(1)*tk;
    
    %计算卫星轨道平面坐标系的坐标
    xk=rk*cos(uk);
    yk=rk*sin(uk);
 
    %计算观测瞬间升交点的经度Lk
	Lk=B4(3)+B5(4)*tk-We*B4(1);

    %计算卫星在轨道平面内的坐标x,y,z
	xyz=[rk*cos(uk);rk*sin(uk);0];

    A=[xk*cos(Lk)-yk*cos(-ik)*sin(Lk);xk*sin(Lk)+yk*cos(-ik)*cos(Lk);yk*sin(ik)];
    
    Rx=[1,0,0;0,cos(-5*pi/180),sin(-5*pi/180);0,-sin(-5*pi/180),cos(-5*pi/180)];
    Rz=[cos(We*tk),sin(We*tk),0;-sin(We*tk),cos(We*tk),0;0,0,1];
    
    
   
    XYZ=vpa([Rz*Rx*A],20);
end

XYZ

二，广播星历文件格式展示

三、广播星历文件解读

 四、代码详解

• 打开文件

• 跳过广播星历中不需要的数据

• 读文件 

 

•  从第26行开始按行读取数据，并把他们存储在B#里

这样转换的原因：

1.存储方式不同：单元格数组一般存储混合数据类型的表格或其他非均匀数据集，而数值数组以规则的网格形式存储数据，所有元素都按照相同的规则排列，转换有利于下标来索引。

2.运算上区别：单元格数组不能直接进行数学运算，数值数组支持各种数学运算，例如加法、乘法、矩阵运算等。

• 定义常量

• 求出轨道长半轴

• 求出平均角速度，时间差，改正平均角速度，平近点角 

•  计算升交距角Phik和真近点角fk

•  将广播星历中的数据读取到变量中

• 计算摄动改正量

• 计算摄动改正后的 升交距角uk、卫星矢径rk和轨道倾角ik 

• 计算卫星在轨道平面坐标系的坐标xk和yk，升交点经度Lk

• 计算GEO卫星在自定义惯性系中的坐标

• 计算GEO卫星在CGCS2000中的坐标 

使用vpa函数来设置有效数字的位数 

• 结果