基于广播星历的北斗定位解算原理（基于C语言和MATLAB实现）

已于 2023-12-15 15:31:36 修改
阅读量4.1k 收藏 87
点赞数 38
分类专栏： GNSS 文章标签： matlab 算法 开源 1024程序员节
于 2023-10-19 14:26:44 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.youkuaiyun.com/weixin_63030308/article/details/133915787
版权

一、利用广播星历解算卫星位置

网上已经有很多解算北斗卫星位置的资料和算法，此文章的目的只是记录自己的学习过程，并且积累一些写文章博客的经验。本文先用C语言解算卫星位置，再用MATLAB绘出卫星三维坐标图。本篇博客所使用的资料和文件都是网络上公开发表且可以找到的资料文件。16参数所解算的文件是RINEX3.04标准；我做的18参数模型采用的是国内的格式，故解算的文件是RINEX3.03标准。

以下是实现方法

我所参考的是《2017中国卫星导航系统管理办公室》和《2019中国卫星导航系统管理办公室》文件，链接放在下面，大家感兴趣的可自行下载，本文所引用的公式及截图均来自此文件。

2017版：http://www.beidou.gov.cn/yw/xwzx/201712/W020171226734109234149.pdf

2019版：http://www.beidou.gov.cn/xt/gfxz/201902/P020190227592987952674.pdf

一、16参数算法

一、前言

16参数解算步骤每次只能解算其中一颗卫星的位置；18参数结算步骤可以直接导入完整的文件，并对其中的北斗卫星进行解算并存入文件中，方便后续导入其他软件做数据研究。

无论是16参数或者是18参数版本，在main函数中参数的顺序是按照（1+6+9）即（时间参数+6个开普勒轨道根数+9个摄动参数）声明的，18参数只是加了两个变化率参数。

卫星型号具体是MEO/GEO/IGSO，可以参考北斗卫星导航系统办公室发表的文件。如下图所示

27af87669afe4875b1ad88db067bdfec.png

16/18参数在广播星历中各自所对应的位置在网上均可找到，此处不做赘述。

 

二、实现步骤（具体每一步讲解会在代码中体现并带有注释）

7299ac3b92dc4394afaafa0c970277b1.png

880f2f40301b49f4b9634e2eb784cc2c.png

c90af0156d97471ba91f965de7d2dc88.png

9606031f37bc4c4e9ad74f185eacffbc.png

三、C语言源码

//**************************************** Message ***********************************//
//技术交流：folkore0118@163.com
//关注优快云博主：“彼稷”
//Author: 彼稷 
//All rights reserved	                               
//----------------------------------------------------------------------------------------
// Target Devices: 		DELL-G15
// Tool Versions:       DEVC++5.11
// File name:           Calculation_BDS_Position
// Last modified Date:  2023年10月12日20:00:00
// Last Version:        
// Descriptions:        广播星历&北斗卫星位置 
//----------------------------------------------------------------------------------------
//****************************************************************************************//
#include <stdio.h>
#include <string.h>
#include <math.h>

#define	Mu			3.986004418e14
#define	Omega_e		7.292115e-5
#define	Pi			3.1415926535898

typedef	struct{
	char	PRN[10];
	int		year;
	int		month;
	int		day;
	int		hour;
	int		minute;
	double	second;
}TIME;
	
	/*--------------------------------------------------------------------------------
	get the BDS time
	*/

int GetsBDSTime(char PRN[10],int year,int month,int day,int hour,int minute,double second,int *WeekNum,double *SecondOfWeek){
	int	DayOfYear	=	0;
	int	DayOfMonth	=	0;
//	/*-------------------------------------
//	year
//	
	for(int DayOfBDS = 2006;DayOfBDS < year;DayOfBDS ++){
		if(((DayOfBDS % 4 == 0)&&(DayOfBDS % 100 != 0))||(DayOfBDS % 400 == 0)){
			DayOfYear	+=	366;
		}
		else{
			DayOfYear	+=	365;
		}
	}
//	/*--------------------------------------
//	month
//	
	for(int DayOfBDS = 1;DayOfBDS < month;DayOfBDS ++){
		if((DayOfBDS == 1)||(DayOfBDS == 3)||(DayOfBDS == 5)||(DayOfBDS == 7)||(DayOfBDS == 8)||(DayOfBDS == 10)||(DayOfBDS == 12)){
			DayOfMonth	+=	31;
		}	
		else if((DayOfBDS == 4)||(DayOfBDS == 6)||(DayOfBDS == 9)||(DayOfBDS == 11)){
			DayOfMonth	+=	30;
		} 
		else{
			if(((year % 4 == 0)&&(year % 100 != 0))||(year % 400 == 0)){
				DayOfMonth	+=	29;
			}
			else{
				DayOfMonth	+=	28;
			}
		}
	}
	int	VariableDay	=	0;
	int	ExtraDay	=	0;
	VariableDay		=	DayOfYear	+	DayOfMonth	+	day -	1;			//总天数 
	*WeekNum		=	VariableDay/7;										//总周数 
	ExtraDay		=	VariableDay - (7*(VariableDay/7));//VariableDay%7
	*SecondOfWeek	=	(double)((ExtraDay * 86400) + (hour * 3600) + (minute * 60) + second);
	//return DayOfBDS;
}


int	main(){
	
	TIME time;
	/*-------------------------------------------------------------------------------------
	parameters
	*/
	/*-----------------------------------------------
	timing
	*/
	double	a0	=	0;				//钟差 
	double	a1	=	0;				//钟速 
	double	a2	=	0;				//钟飘 
	double	toe	=	0;				//星历表参考历元 
	double	toc	=	0;				//参考时间 
	/*-----------------------------------------------
	Kepler orbit parameters
	*/
	double	i0		=	0;			//轨道倾角 
	double	e		=	0;			//轨道偏心率 
	double	OMEGA	=	0;			//Ω升交点赤经 
	double	omega	=	0;			//w近地点角距 
	double	M0		=	0;			//M0toe时刻的平近点角 
	double	sqrt_a	=	0;			//轨道长半径的平方根 
	/*-----------------------------------------------
	perturbation parameters
	*/
	double	OMEGADOT=	0;			//升交点赤经变化率 
	double	IDOT	=	0;			//轨道倾角变化率 
	double	Cuc		=	0;			//纬度幅角的余弦调和项 
	double	Cic		=	0;			//轨道倾角的余弦调和项 
	double	Crc		=	0;			//轨道半径的余弦调和项 
	double	Cus		=	0;			//纬度幅角的正弦调和项 
	double	Cis		=	0;			//轨道倾角的正弦调和项 
	double	Crs		=	0;			//轨道半径的正弦调和项 
	double	DELTA_n	=	0;			//平均角速度摄动量 
	/*-----------------------------------------------
	extra parameters
	*/
	double	AODE	=	0;
	
	/*----------------------------------------------------------------------------------------
	read file
	*/
	FILE *pfile			=	fopen("文件路径（相对路径或绝对路径均可）","r");
	FILE *pfileWrite	=	fopen("文件路径（相对路径或绝对路径均可）","w");
	if(!pfile || !pfileWrite){
		printf("/-----------------------------------------------------------------------------------------/\n");
		printf("\topen the file error!\n");
	}
	else{
		printf("/-----------------------------------------------------------------------------------------/\n");
		printf("\topen the file success!\n");
	}
	printf("/-----------------------------------------------------------------------------------------/\n");
	/*---------------------------------------------
	first line
	*/
	fscanf(pfile,"%s %d %d %d %d %d %lf %lf %lf %lf\n",&time.PRN,&time.year,&time.month,&time.day,&time.hour,&time.minute,&time.second,&a0,&a1,&a2);
//	printf("%s %d %d %d %d %d %.2lf %.15lf %.15lf %.15lf\n",time.PRN,time.year,time.month,time.day,time.hour,time.minute,time.second,a0,a1,a2);
	printf("\t%s %d %d %d %d %d %.2lf %e\t%e\t%e\t\n",time.PRN,time.year,time.month,time.day,time.hour,time.minute,time.second,a0,a1,a2);
	/*----------------------------------------------
	second line
	*/
	fscanf(pfile,"%lf %lf %lf %lf\n",&AODE,&Crs,&DELTA_n,&M0);
//	printf("    %.15lf %.15lf %.15lf %.15lf\n",AODE,Crs,DELTA_n,M0);
	printf("\t    %e     %e\t%e\t%e\t\n",AODE,Crs,DELTA_n,M0);
	/*----------------------------------------------
	third line
	*/
	fscanf(pfile,"%lf %lf %lf %lf\n",&Cuc,&e,&Cus,&sqrt_a);
//	printf("    %.15lf %.15lf %.15lf %.15lf\n",Cuc,e,Cus,sqrt_a);
	printf("\t    %e     %e\t%e\t%e\t\n",Cuc,e,Cus,sqrt_a);
	/*----------------------------------------------
	fourth line
	*/	
	fscanf(pfile,"%lf %lf %lf %lf\n",&toe,&Cic,&OMEGA,&Cis);
//	printf("    %.15lf %.15lf %.15lf %.15lf\n",toe,Cic,OMEGA,Cis);
	printf("\t    %e     %e\t%e\t%e\t\n",toe,Cic,OMEGA,Cis);
	/*----------------------------------------------
	fifth line
	*/
	fscanf(pfile,"%lf %lf %lf %lf\n",&i0,&Crc,&omega,&OMEGADOT);
//	printf("    %.15lf %.15lf %.15lf %.15lf\n",i0,Crc,omega,OMEGADOT);
	printf("\t    %e     %e\t%e\t%e\t\n",i0,Crc,omega,OMEGADOT);
	/*----------------------------------------------
	sixth
	*/
	fscanf(pfile,"%lf\n",&IDOT);
//	printf("    %.15lf\n",IDOT);
	printf("\t    %e\n",IDOT);
	printf("/-----------------------------------------------------------------------------------------/\n");
	
	/*-------------------------------------------------------------------------------------
	calculate values
	*/
	int year,month,day,hour,minute,WeekNum;
	char PRN[10];
	double second,SecondOfWeek;
	GetsBDSTime(time.PRN,time.year,time.month,time.day,time.hour,time.minute,time.second,&WeekNum,&SecondOfWeek);
	
	printf("\t周内秒	= %.20lf\n",SecondOfWeek);
	/*---------------------------------------------
	一、计算平均角速度n 
	*/
	
	double	n0,n;
	n0	=	sqrt(Mu)/pow(sqrt_a,3);
	n	=	n0	+	DELTA_n;
	printf("\tn0	= %.20lf\n",n0);
	printf("\t平均角速度n	= %.20lf\n",n);
	/*---------------------------------------------
	二、计算归化时间 
	*/
	double	t,tk;
//	double	delta_t;
	double trans_t	=	0;
	printf("\t输入参考时间t:(小于60分钟)\t====>");
	scanf("%lf",&t);
//	delta_t	=	a0	+	a1 * ((t * 60 + SecondOfWeek) - toc)	+	a2 * pow(((t * 60 + SecondOfWeek) - toc),2);
	trans_t	=	t * 60	+	SecondOfWeek;
	tk	=	trans_t	-	toe	-	14;
	if(tk > 302400){
		tk	=	tk	-	604800;
	}
	else if(tk < -302400){
		tk	=	tk	+	604800;
	}
	printf("\t归化时间tk	= %.20lf\n",tk);	
	/*---------------------------------------------
	三、计算平近点角Mk 
	*/
	double	Mk;
	Mk	=	M0	+	n * tk;
	printf("\t平近点角Mk	= %.20lf\n",Mk);
	/*---------------------------------------------
	四、计算偏近点角（迭代） 
	*/
	double	Ek;
	int	i	=	0;
	Ek	=	Mk;
	do{
		Ek	=	Mk	+	e * sin(Ek);
		i++;
	}while(i < 7);
	printf("\t偏近点角Ek	= %.20lf\n",Ek);
	/*---------------------------------------------
	五、计算真近点角 
	*/
	double	Vk;
	Vk	=	atan(((sqrt(1 - e * e )*sin(Ek)))/(cos(Ek) - e));		//弧度 
	printf("\t真近点角Vk	= %.20lf\n",Vk);
	/*---------------------------------------------
	六、计算升交距角 
	*/
	double	Phik;
	Phik	=	Vk	+	omega;
	printf("\t升交距角φk	= %.20lf\n",Phik);
	/*---------------------------------------------
	七、计算摄动改正量 
	*/
	double	delta_u,delta_r,delta_i;
	delta_u	=	Cuc * cos(2 * Phik)	+	Cus * sin(2 * Phik); 
	delta_r	=	Crc * cos(2 * Phik)	+	Crs * sin(2 * Phik);
	delta_i	=	Cic * cos(2 * Phik)	+	Cis * sin(2 * Phik);
	printf("\t摄动改正量δu	= %.20lf\n",delta_u);
	printf("\t摄动改正量δr	= %.20lf\n",delta_r);
	printf("\t摄动改正量δi	= %.20lf\n",delta_i);
	/*---------------------------------------------
	八、计算摄动改正后的升交距角uk、卫星矢径rk、轨道倾角ik 
	*/
	double	uk,rk,ik;
	uk	=	Phik	+	delta_u;
	rk	=	(sqrt_a * sqrt_a) * (1 - e * cos(Ek))	+	delta_r;
	ik	=	i0	+	delta_i	+	IDOT * tk;	
	printf("\t摄动改正后的升交距角uk	= %.20lf\n",uk);
	printf("\t摄动改正后的卫星矢径rk	= %.20lf\n",rk);
	printf("\t摄动改正后的轨道倾角ik	= %.20lf\n",ik);
	/*---------------------------------------------
	九、计算轨道平面坐标系的坐标 
	*/
	double	Xk,Yk;
	Xk	=	rk * cos(uk);
	Yk	=	rk * sin(uk);
	printf("\t轨道平面坐标系X	= %.20lf\n",Xk);
	printf("\t轨道平面坐标系Y	= %.20lf\n",Yk);

	if(strcmp(time.PRN,"C01")==0||strcmp(time.PRN,"C02")==0||strcmp(time.PRN,"C03")==0||strcmp(time.PRN,"C04")==0
	||strcmp(time.PRN,"C05")==0||strcmp(time.PRN,"C59")==0||strcmp(time.PRN,"C60")==0){
		/*---------------------------------------------
		十、计算观测时刻的升交点经度 
		*/
		printf("\tPRN	= %s,为GEO卫星\n",time.PRN);
		double	OMEGA_k;
		OMEGA_k	=	OMEGA	+	OMEGADOT * tk - Omega_e * toe;
		printf("\t升交点经度Ωk	= %.20lf\n",OMEGA_k);
		/*--------------------------------------------
		十一、计算自定义惯性系中的坐标 
		*/
		double	ApostOfX,ApostOfY,ApostOfZ; 
		double	AApostOfX,AApostOfY,AApostOfZ; 
		double	ResultX,ResultY,ResultZ,Resulta0;
		ApostOfX	=	Xk * cos(OMEGA_k)	-	Yk * cos(ik) * sin(OMEGA_k);
		ApostOfY	=	Xk * sin(OMEGA_k)	+	Yk * cos(ik) * cos(OMEGA_k);
		ApostOfZ	=	Yk * sin(ik);
		/*--------------------------------------------
		十一、计算自CGCS2000的坐标 
		*/
		double	RMatrixZ[3][3]	=	{{cos(Omega_e*tk),sin(Omega_e*tk),0},{-sin(Omega_e*tk),cos(Omega_e*tk),0},{0,0,1}};
		double	RMatrixX[3][3]	=	{{1,0,0},{0,cos(-Pi/180*5),sin(-Pi/180*5)},{0,-sin(-Pi/180*5),cos(-Pi/180*5)}};
		double	RMatrixGk[3][1]	=	{{ApostOfX},{ApostOfY},{ApostOfZ}};
		double	ZX[3][3];
		double	ZXGk[3][1];
		/*第一次矩阵*/
		for(int i = 0;i < 3;i++){
			for(int j = 0;j < 3;j++){
				ZX[i][j]	=	0;
				for(int k = 0;k < 3;k++){
					ZX[i][j]	+=	RMatrixZ[i][k] * RMatrixX[k][j];
				}	
			}	
		}
		printf("\t矩阵RMatrixZ和矩阵RMatrixX相乘的结果为：\n");
		for(int i = 0;i < 3;i++){
			for(int j = 0;j < 3;j++){
				printf("\t%.10lf\t",ZX[i][j]);
			}
			printf("\n");
		}
		/*第二次矩阵*/
		for(int i = 0;i < 3;i++){
			for(int j = 0;j < 1;j++){
				ZXGk[i][j]	=	0;
				for(int k = 0;k < 3;k++){
					ZXGk[i][j]	+=	ZX[i][k] * RMatrixGk[k][j];
				}	
			}	
		}
		printf("\t矩阵ZX和矩阵RMatrixGk相乘的结果为：\n");
		for(int i = 0;i < 3;i++){
			for(int j = 0;j < 1;j++){
				printf("\t%.10lf\t",ZXGk[i][j]);
			}
			printf("\n");
		}
		AApostOfX	=	ZXGk[0][0];
		AApostOfY	=	ZXGk[1][0];
		AApostOfZ	=	ZXGk[2][0];
		printf("\tCoordinateOfX	=	%.20lf\n",AApostOfX);
		printf("\tCoordinateOfY	=	%.20lf\n",AApostOfY);
		printf("\tCoordinateOfZ	=	%.20lf\n",AApostOfZ);
		ResultX		=	AApostOfX / 1000.0;
		ResultY		=	AApostOfY / 1000.0;
		ResultZ		=	AApostOfZ / 1000.0;
		Resulta0	=	a0 * 1000000.0;
		printf("\tResultX		=	%.6lf\n",ResultX);
		printf("\tResultY		=	%.6lf\n",ResultY);
		printf("\tResultZ		=	%.6lf\n",ResultZ);
		printf("\tResulta0	=	%.6lf\n",Resulta0);
			/*-----------------------------------------------------------------------------------------------
		将数据写入文件并保存 
		*/
		printf("/-----------------------------------------------------------------------------------------/\n");
		fprintf(pfileWrite,"%s %.6lf %.6lf %.6lf %.6lf\n",time.PRN,ResultX,ResultY,ResultZ,Resulta0);
		printf("/-----------------------------------------------------------------------------------------/\n");
	}	
	
	else if(strcmp(time.PRN,"C11")==0||strcmp(time.PRN,"C12")==0||strcmp(time.PRN,"C14")==0||strcmp(time.PRN,"C19")==0||
	strcmp(time.PRN,"C20")==0||strcmp(time.PRN,"C21")==0||strcmp(time.PRN,"C22")==0||strcmp(time.PRN,"C23")==0||
	strcmp(time.PRN,"C24")==0||strcmp(time.PRN,"C25")==0||strcmp(time.PRN,"C26")==0||strcmp(time.PRN,"C27")==0||
	strcmp(time.PRN,"C28")==0||strcmp(time.PRN,"C29")==0||strcmp(time.PRN,"C30")==0||strcmp(time.PRN,"C32")==0||
	strcmp(time.PRN,"C33")==0||strcmp(time.PRN,"C34")==0||strcmp(time.PRN,"C35")==0||strcmp(time.PRN,"C36")==0||
	strcmp(time.PRN,"C37")==0||strcmp(time.PRN,"C41")==0||strcmp(time.PRN,"C42")==0||strcmp(time.PRN,"C43")==0||
	strcmp(time.PRN,"C44")==0||strcmp(time.PRN,"C45")==0||strcmp(time.PRN,"C46")==0||strcmp(time.PRN,"C06")==0||
	strcmp(time.PRN,"C07")==0||strcmp(time.PRN,"C08")==0||strcmp(time.PRN,"C09")==0||strcmp(time.PRN,"C10")==0||
	strcmp(time.PRN,"C13")==0||strcmp(time.PRN,"C16")==0||strcmp(time.PRN,"C38")==0||strcmp(time.PRN,"C39")==0||
	strcmp(time.PRN,"C40")==0){
		/*---------------------------------------------
		十、计算观测时刻的升交点经度 
		*/
		printf("\tPRN = %s,为MEO/IGSO卫星\n",time.PRN);
		double	OMEGA_k;
		OMEGA_k	=	OMEGA	+	(OMEGADOT	-	Omega_e) * tk	-	Omega_e * toe;
		printf("\t升交点经度Ωk = %.20lf\n",OMEGA_k);
		printf("/-----------------------------------------------------------------------------------------/\n");
		/*---------------------------------------------
		十一、计算CGCS2000的坐标 	apostrophe--->撇 	coordinate--->坐标 
		*/
		double	ApostOfX,ApostOfY,ApostOfZ; 
		double	ResultX,ResultY,ResultZ,Resulta0;
		ApostOfX	=	Xk * cos(OMEGA_k)	-	Yk * cos(ik) * sin(OMEGA_k);
		ApostOfY	=	Xk * sin(OMEGA_k)	+	Yk * cos(ik) * cos(OMEGA_k);
		ApostOfZ	=	Yk * sin(ik);
		printf("\tCoordinateOfX	=	%.20lf\n",ApostOfX);
		printf("\tCoordinateOfY	=	%.20lf\n",ApostOfY);
		printf("\tCoordinateOfZ	=	%.20lf\n",ApostOfZ);
		ResultX		=	ApostOfX / 1000.0;
		ResultY		=	ApostOfY / 1000.0;
		ResultZ		=	ApostOfZ / 1000.0;
		Resulta0	=	a0 * 1000000.0;
		printf("\tResultX		=	%.6lf\n",ResultX);
		printf("\tResultY		=	%.6lf\n",ResultY);
		printf("\tResultZ		=	%.6lf\n",ResultZ);
		printf("\tResulta0		=	%.6lf\n",Resulta0);
	/*-----------------------------------------------------------------------------------------------
	将数据写入文件并保存 
	*/
	printf("/-----------------------------------------------------------------------------------------/\n");
	fprintf(pfileWrite,"%s %.6lf %.6lf %.6lf %.6lf\n",time.PRN,ResultX,ResultY,ResultZ,Resulta0);
	printf("/-----------------------------------------------------------------------------------------/\n");
	}
	fclose(pfile);
	fclose(pfileWrite);
	return 0;
}

值得注意的是，toc转到BDT是从2006年1月1日（周日）作为起始点；因为BDT和GPST相差14秒，所以计算归化时间tk时，需要人为减去14秒，调整到和GPST同步。计算tk时，t相当于是对卫星位置做预测，且北斗卫星每隔一个小时更新一次星历，故t的取值小于60分钟。

二、18参数算法

一、前言

18参数模型相对比16参数改动不是特别大，一个是公式的变动，另一个是广播星历对应参数的位置部分有变动。

公式的变动主要体现在下图，具体的还是参考中国卫星导航系统管理办公室提供的官方文件。

c70f71ae60b0412cb37d843d436cf9a9.png

二、实现步骤（具体每一步讲解会在代码中体现并带有注释）

ae0a9215ccfb4695987949967240ca75.png

93e0bd8a0cf94ad6a9f4e9faeacd737b.png

三、C语言源码

本代码仅适用于计算北斗卫星国内格式的文件，国外格式的文件需要做相对应的修正。

本代码可以直接读取整个大文件，并计算出所有北斗卫星的位置。

//**************************************** Message ***********************************//
//技术交流：folkore0118@163.com
//关注优快云博主：“彼稷”
//Author: 彼稷 
//All rights reserved	                               
//----------------------------------------------------------------------------------------
// Target Devices: 		DELL-G15
// Tool Versions:       DEVC++5.11
// File name:           Calculation_BDS_Position
// Last modified Date:  2023年10月18日20:00:00
// Last Version:        
// Descriptions:        广播星历&北斗卫星位置 &18参数 
//----------------------------------------------------------------------------------------
//****************************************************************************************//
#include <stdio.h>
#include <string.h>
#include <math.h>

#define	Mu			3.986004418e14
#define	Omega_e		7.292115e-5
#define	Pi			3.1415926535898
#define	Aref_MEO	27906100
#define	Aref_GEO	42162200

typedef	struct{
	char	PRN[10];
	int		year;
	int		month;
	int		day;
	int		hour;
	int		minute;
	double	second;
}TIME;
	
	/*--------------------------------------------------------------------------------
	get the BDS time
	*/

int GetsBDSTime(char PRN[10],int year,int month,int day,int hour,int minute,double second,int *WeekNum,double *SecondOfWeek){
	int	DayOfYear	=	0;
	int	DayOfMonth	=	0;
//	/*-------------------------------------
//	year
//	
	for(int DayOfBDS = 2006;DayOfBDS < year;DayOfBDS ++){
		if(((DayOfBDS % 4 == 0)&&(DayOfBDS % 100 != 0))||(DayOfBDS % 400 == 0)){
			DayOfYear	+=	366;
		}
		else{
			DayOfYear	+=	365;
		}
	}
//	/*--------------------------------------
//	month
//	
	for(int DayOfBDS = 1;DayOfBDS < month;DayOfBDS ++){
		if((DayOfBDS == 1)||(DayOfBDS == 3)||(DayOfBDS == 5)||(DayOfBDS == 7)||(DayOfBDS == 8)||(DayOfBDS == 10)||(DayOfBDS == 12)){
			DayOfMonth	+=	31;
		}	
		else if((DayOfBDS == 4)||(DayOfBDS == 6)||(DayOfBDS == 9)||(DayOfBDS == 11)){
			DayOfMonth	+=	30;
		} 
		else{
			if(((year % 4 == 0)&&(year % 100 != 0))||(year % 400 == 0)){
				DayOfMonth	+=	29;
			}
			else{
				DayOfMonth	+=	28;
			}
		}
	}
	int	VariableDay	=	0;
	int	ExtraDay	=	0;
	VariableDay		=	DayOfYear	+	DayOfMonth	+	day -	1;			//总天数 
	*WeekNum		=	VariableDay/7;										//总周数 
	ExtraDay		=	VariableDay - (7*(VariableDay/7));//VariableDay%7
	*SecondOfWeek	=	(double)((ExtraDay * 86400) + (hour * 3600) + (minute * 60) + second);
	//return DayOfBDS;
}


int	main(){
	
	TIME time;
	/*-------------------------------------------------------------------------------------
	parameters
	*/
	/*-----------------------------------------------
	timing
	*/
	double	a0	=	0;				//钟差 
	double	a1	=	0;				//钟速 
	double	a2	=	0;				//钟飘 
	double	toe	=	0;				//星历表参考历元 
	double	toc	=	0;				//参考时间 
	/*-----------------------------------------------
	Kepler orbit parameters
	*/
	double	i0		=	0;			//轨道倾角 
	double	e		=	0;			//轨道偏心率 
	double	OMEGA	=	0;			//Ω升交点赤经 
	double	omega	=	0;			//w近地点角距 
	double	M0		=	0;			//M0toe时刻的平近点角 
	double	DELTA_A	=	0;			//参考时刻长半轴较参考值的偏差 
	/*-----------------------------------------------
	perturbation parameters
	*/
	double	OMEGADOT=	0;			//升交点赤经变化率 
	double	IDOT	=	0;			//轨道倾角变化率 
	double	Cuc		=	0;			//纬度幅角的余弦调和项 
	double	Cic		=	0;			//轨道倾角的余弦调和项 
	double	Crc		=	0;			//轨道半径的余弦调和项 
	double	Cus		=	0;			//纬度幅角的正弦调和项 
	double	Cis		=	0;			//轨道倾角的正弦调和项 
	double	Crs		=	0;			//轨道半径的正弦调和项 
	double	DELTA_n0=	0;			//平均角速度摄动量 
	/*-----------------------------------------------
	extra parameters
	*/
	double	AODE	=	0;
	double	ADOT	=	0;
	double	data	=	0;
	double	BDTWeek	=	0;
	double	PrecisionIndex	=	0;
	double	HealthStatus	=	0;
	double	IGD		=	0;
	double	ISC		=	0;
	double	BDSLaunchTime	=	0;
	double	IODC	=	0;
	double	DELTA_n0DOT		=	0;
	double	SatelliteOrbit	=	0;
	/*----------------------------------------------------------------------------------------
	read file
	*/
//	FILE *pfile			=	fopen("文件路径（相对路径或绝对路径均可）","r");
//	FILE *pfileWrite	=	fopen("文件路径（相对路径或绝对路径均可）","w");
	FILE *pfile			=	fopen("ReadFile_18.txt","r");
	FILE *pfileWrite	=	fopen("WriteFile_18.txt","w");
	if(!pfile || !pfileWrite){
		printf("/-----------------------------------------------------------------------------------------/\n");
		printf("\topen the file error!\n");
	}
	else{
		printf("/-----------------------------------------------------------------------------------------/\n");
		printf("\topen the file success!\n");
	}
	printf("/-----------------------------------------------------------------------------------------/\n");
	/*---------------------------------------------
		first line
	*/
	
	while(1){
		fscanf(pfile,"%s %d %d %d %d %d %lf %lf %lf %lf\n",&time.PRN,&time.year,&time.month,&time.day,&time.hour,&time.minute,&time.second,&a0,&a1,&a2);
		if(feof(pfile)){
			printf("Reached end of file.\n");
			break;
		}
		else if(ferror(pfile)){
			printf("Error occurred while reading file.\n");
			continue;
		}
		else if((strcmp(time.PRN,"C01")==0||strcmp(time.PRN,"C02")==0||strcmp(time.PRN,"C03")==0||strcmp(time.PRN,"C04")==0||
			strcmp(time.PRN,"C05")==0||strcmp(time.PRN,"C59")==0||strcmp(time.PRN,"C60")==0||strcmp(time.PRN,"C11")==0||
			strcmp(time.PRN,"C12")==0||strcmp(time.PRN,"C14")==0||strcmp(time.PRN,"C19")==0||strcmp(time.PRN,"C20")==0||
			strcmp(time.PRN,"C21")==0||strcmp(time.PRN,"C22")==0||strcmp(time.PRN,"C23")==0||strcmp(time.PRN,"C24")==0||
			strcmp(time.PRN,"C25")==0||strcmp(time.PRN,"C26")==0||strcmp(time.PRN,"C27")==0||strcmp(time.PRN,"C28")==0||
			strcmp(time.PRN,"C29")==0||strcmp(time.PRN,"C30")==0||strcmp(time.PRN,"C32")==0||strcmp(time.PRN,"C33")==0||
			strcmp(time.PRN,"C34")==0||strcmp(time.PRN,"C35")==0||strcmp(time.PRN,"C36")==0||strcmp(time.PRN,"C37")==0||
			strcmp(time.PRN,"C41")==0||strcmp(time.PRN,"C42")==0||strcmp(time.PRN,"C43")==0||strcmp(time.PRN,"C44")==0||
			strcmp(time.PRN,"C45")==0||strcmp(time.PRN,"C46")==0||strcmp(time.PRN,"C06")==0||strcmp(time.PRN,"C07")==0||
			strcmp(time.PRN,"C08")==0||strcmp(time.PRN,"C09")==0||strcmp(time.PRN,"C10")==0||strcmp(time.PRN,"C13")==0||
			strcmp(time.PRN,"C16")==0||strcmp(time.PRN,"C38")==0||strcmp(time.PRN,"C39")==0||strcmp(time.PRN,"C40")==0
			)&&(time.year == 2023)){
		/*else if(strcmp(time.PRN,"C01")==0||strcmp(time.PRN,"C02")==0||strcmp(time.PRN,"C03")==0||strcmp(time.PRN,"C04")==0
				||strcmp(time.PRN,"C05")==0||strcmp(time.PRN,"C59")==0||strcmp(time.PRN,"C60")==0){*/
				printf("\t%s %d %d %d %d %d %.2lf %e\t%e\t%e\t\n",time.PRN,time.year,time.month,time.day,
				time.hour,time.minute,time.second,a0,a1,a2);
				/*----------------------------------------------
				second line
				*/
				fscanf(pfile,"%lf %lf %lf %lf\n",&AODE,&Crs,&DELTA_n0,&M0);
				printf("\t    %e     %e\t%e\t%e\t\n",AODE,Crs,DELTA_n0,M0);
				/*----------------------------------------------
				third line
				*/
				fscanf(pfile,"%lf %lf %lf %lf\n",&Cuc,&e,&Cus,&DELTA_A);
				printf("\t    %e     %e\t%e\t%e\t\n",Cuc,e,Cus,DELTA_A);
				/*----------------------------------------------
				fourth line
				*/	
				fscanf(pfile,"%lf %lf %lf %lf\n",&toe,&Cic,&OMEGA,&Cis);
				printf("\t    %e     %e\t%e\t%e\t\n",toe,Cic,OMEGA,Cis);
				/*----------------------------------------------
				fifth line
				*/
				fscanf(pfile,"%lf %lf %lf %lf\n",&i0,&Crc,&omega,&OMEGADOT);
				printf("\t    %e     %e\t%e\t%e\t\n",i0,Crc,omega,OMEGADOT);
				/*----------------------------------------------
				sixth
				*/
				fscanf(pfile,"%lf %lf %lf %lf\n",&IDOT,&data,&BDTWeek,&ADOT);
				printf("\t    %e     %e\t%e\t%e\t\n",IDOT,data,BDTWeek,ADOT);
				/*----------------------------------------------
				seventh
				*/
				fscanf(pfile,"%lf %lf %lf %lf\n",&PrecisionIndex,&HealthStatus,&IGD,&ISC);
				printf("\t    %e     %e\t%e\t%e\t\n",PrecisionIndex,HealthStatus,IGD,ISC);
				/*----------------------------------------------
				eighth
				*/
				fscanf(pfile,"%lf %lf %lf %lf\n",&BDSLaunchTime,&IODC,&DELTA_n0DOT,&SatelliteOrbit);
				printf("\t    %e     %e\t%e\t%e\t\n",BDSLaunchTime,IODC,DELTA_n0DOT,SatelliteOrbit);
				printf("/-----------------------------------------------------------------------------------------/\n");
				
				/*-------------------------------------------------------------------------------------
				calculate values
				*/
				int year,month,day,hour,minute,WeekNum;
				char PRN[10];
				double second,SecondOfWeek;
				GetsBDSTime(time.PRN,time.year,time.month,time.day,time.hour,time.minute,time.second,&WeekNum,&SecondOfWeek);
				
				printf("\t周内秒	= %.20lf\n",SecondOfWeek);
				/*---------------------------------------------
				一、计算归化时间 
				*/
				double	t	=	0.0;
				double	tk;
				//	double	delta_t;
				double trans_t	=	0;
				printf("\t输入参考时间t:(小于60分钟)\t====>0\n");
				//scanf("%lf",&t);
				//	delta_t	=	a0	+	a1 * ((t * 60 + SecondOfWeek) - toc)	+	a2 * pow(((t * 60 + SecondOfWeek) - toc),2);
				trans_t	=	t * 60	+	SecondOfWeek;
				tk	=	trans_t	-	toe	-	14;
				if(tk > 302400){
					tk	=	tk	-	604800;
				}
				else if(tk < -302400){
					tk	=	tk	+	604800;
				}
				printf("\t归化时间tk	= %.20lf\n",tk);
				/*---------------------------------------------
				二、计算半长轴 
				*/
				double	A0,Ak;
				if(strcmp(time.PRN,"C01")==0||strcmp(time.PRN,"C02")==0||strcmp(time.PRN,"C03")==0||strcmp(time.PRN,"C04")==0
				||strcmp(time.PRN,"C05")==0||strcmp(time.PRN,"C59")==0||strcmp(time.PRN,"C60")==0){		 //GEO
					A0	=	Aref_GEO	+	DELTA_A;
					printf("\tPRN	= %s,为GEO卫星\n",time.PRN);
				}
				else if(strcmp(time.PRN,"C06")==0||strcmp(time.PRN,"C07")==0||strcmp(time.PRN,"C08")==0||strcmp(time.PRN,"C09")==0
				||strcmp(time.PRN,"C10")==0||strcmp(time.PRN,"C13")==0||strcmp(time.PRN,"C16")==0||strcmp(time.PRN,"C38")==0||
				strcmp(time.PRN,"C39")==0||strcmp(time.PRN,"C40")==0){		 //IGSO
					A0	=	Aref_GEO	+	DELTA_A;
					printf("\tPRN = %s,为IGSO卫星\n",time.PRN);
				}
				else if(strcmp(time.PRN,"C11")==0||strcmp(time.PRN,"C12")==0||strcmp(time.PRN,"C14")==0||strcmp(time.PRN,"C19")==0||
				strcmp(time.PRN,"C20")==0||strcmp(time.PRN,"C21")==0||strcmp(time.PRN,"C22")==0||strcmp(time.PRN,"C23")==0||
				strcmp(time.PRN,"C24")==0||strcmp(time.PRN,"C25")==0||strcmp(time.PRN,"C26")==0||strcmp(time.PRN,"C27")==0||
				strcmp(time.PRN,"C28")==0||strcmp(time.PRN,"C29")==0||strcmp(time.PRN,"C30")==0||strcmp(time.PRN,"C32")==0||
				strcmp(time.PRN,"C33")==0||strcmp(time.PRN,"C34")==0||strcmp(time.PRN,"C35")==0||strcmp(time.PRN,"C36")==0||
				strcmp(time.PRN,"C37")==0||strcmp(time.PRN,"C41")==0||strcmp(time.PRN,"C42")==0||strcmp(time.PRN,"C43")==0||
				strcmp(time.PRN,"C44")==0||strcmp(time.PRN,"C45")==0||strcmp(time.PRN,"C46")==0){
					A0	=	Aref_MEO	+	DELTA_A;
					printf("\tPRN = %s,为MEO卫星\n",time.PRN);
				}
				Ak	=	A0	+	ADOT * tk;
				/*---------------------------------------------
				三、计算平均角速度n 
				*/
				
				double	n0,DELTA_nA,nA;
				n0			=	sqrt(Mu/pow(A0,3));
				DELTA_nA	=	DELTA_n0	+	DELTA_n0DOT * tk * 0.5;
				nA			=	n0	+	DELTA_nA;
				printf("\t参考时刻平均角速度n0	= %.20lf\n",n0);
				printf("\t平均角速度的偏差n		= %.20lf\n",DELTA_nA);
				printf("\t改正后的平均角速度n	= %.20lf\n",nA);
				/*---------------------------------------------
				四、计算平近点角Mk 
				*/
				double	Mk;
				Mk	=	M0	+	nA * tk;
				printf("\t平近点角Mk	= %.20lf\n",Mk);
				/*---------------------------------------------
				五、计算偏近点角（迭代） 
				*/
				double	Ek;
				int	i	=	0;
				Ek	=	Mk;
				do{
					Ek	=	Mk	+	e * sin(Ek);
					i++;
				}while(i < 7);
				printf("\t偏近点角Ek	= %.20lf\n",Ek);
				/*---------------------------------------------
				六、计算真近点角 
				*/
				double	Vk;
				Vk	=	atan(((sqrt(1 - e * e )*sin(Ek)))/(cos(Ek) - e));		//弧度 
				printf("\t真近点角Vk	= %.20lf\n",Vk);
				/*---------------------------------------------
				七、计算纬度幅值 
				*/
				double	Phik;
				Phik	=	Vk	+	omega;
				printf("\t升交距角φk	= %.20lf\n",Phik);
				/*---------------------------------------------
				八、计算摄动改正量 
				*/
				double	delta_u,delta_r,delta_i;
				delta_u	=	Cuc * cos(2 * Phik)	+	Cus * sin(2 * Phik); 
				delta_r	=	Crc * cos(2 * Phik)	+	Crs * sin(2 * Phik);
				delta_i	=	Cic * cos(2 * Phik)	+	Cis * sin(2 * Phik);
				printf("\t摄动改正量δu	= %.20lf\n",delta_u);
				printf("\t摄动改正量δr	= %.20lf\n",delta_r);
				printf("\t摄动改正量δi	= %.20lf\n",delta_i);
				/*---------------------------------------------
				九、计算摄动改正后的升交距角uk、卫星矢径rk、轨道倾角ik 
				*/
				double	uk,rk,ik;
				uk	=	Phik	+	delta_u;
				rk	=	Ak * (1 - e * cos(Ek))	+	delta_r;
				ik	=	i0	+	delta_i	+	IDOT * tk;	
				printf("\t摄动改正后的升交距角uk	= %.20lf\n",uk);
				printf("\t摄动改正后的卫星矢径rk	= %.20lf\n",rk);
				printf("\t摄动改正后的轨道倾角ik	= %.20lf\n",ik);
				/*---------------------------------------------
				十、计算轨道平面坐标系的坐标 
				*/
				double	Xk,Yk;
				Xk	=	rk * cos(uk);
				Yk	=	rk * sin(uk);
				printf("\t轨道平面坐标系X	= %.20lf\n",Xk);
				printf("\t轨道平面坐标系Y	= %.20lf\n",Yk);
			
				if(strcmp(time.PRN,"C01")==0||strcmp(time.PRN,"C02")==0||strcmp(time.PRN,"C03")==0||strcmp(time.PRN,"C04")==0
				||strcmp(time.PRN,"C05")==0||strcmp(time.PRN,"C59")==0||strcmp(time.PRN,"C60")==0){
					/*---------------------------------------------
					十、计算观测时刻的升交点经度 
					*/
					double	OMEGA_k;
					OMEGA_k	=	OMEGA	+	OMEGADOT * tk - Omega_e * toe;
					printf("\t升交点经度Ωk	= %.20lf\n",OMEGA_k);
					/*--------------------------------------------
					十一、计算自定义惯性系中的坐标 
					*/
					double	ApostOfX,ApostOfY,ApostOfZ; 
					double	AApostOfX,AApostOfY,AApostOfZ; 
					double	ResultX,ResultY,ResultZ,Resulta0;
					ApostOfX	=	Xk * cos(OMEGA_k)	-	Yk * cos(ik) * sin(OMEGA_k);
					ApostOfY	=	Xk * sin(OMEGA_k)	+	Yk * cos(ik) * cos(OMEGA_k);
					ApostOfZ	=	Yk * sin(ik);
					/*--------------------------------------------
					十一、计算自CGCS2000的坐标 
					*/
					double	RMatrixZ[3][3]	=	{{cos(Omega_e*tk),sin(Omega_e*tk),0},{-sin(Omega_e*tk),cos(Omega_e*tk),0},{0,0,1}};
					double	RMatrixX[3][3]	=	{{1,0,0},{0,cos(-Pi/180*5),sin(-Pi/180*5)},{0,-sin(-Pi/180*5),cos(-Pi/180*5)}};
					double	RMatrixGk[3][1]	=	{{ApostOfX},{ApostOfY},{ApostOfZ}};
					double	ZX[3][3];
					double	ZXGk[3][1];
					/*第一次矩阵*/
					for(int i = 0;i < 3;i++){
						for(int j = 0;j < 3;j++){
							ZX[i][j]	=	0;
							for(int k = 0;k < 3;k++){
								ZX[i][j]	+=	RMatrixZ[i][k] * RMatrixX[k][j];
							}	
						}	
					}
					printf("\t矩阵RMatrixZ和矩阵RMatrixX相乘的结果为：\n");
					for(int i = 0;i < 3;i++){
						for(int j = 0;j < 3;j++){
							printf("\t%.10lf\t",ZX[i][j]);
						}
						printf("\n");
					}
					/*第二次矩阵*/
					for(int i = 0;i < 3;i++){
						for(int j = 0;j < 1;j++){
							ZXGk[i][j]	=	0;
							for(int k = 0;k < 3;k++){
								ZXGk[i][j]	+=	ZX[i][k] * RMatrixGk[k][j];
							}	
						}	
					}
					printf("\t矩阵ZX和矩阵RMatrixGk相乘的结果为：\n");
					for(int i = 0;i < 3;i++){
						for(int j = 0;j < 1;j++){
							printf("\t%.10lf\t",ZXGk[i][j]);
						}
						printf("\n");
					}
					AApostOfX	=	ZXGk[0][0];
					AApostOfY	=	ZXGk[1][0];
					AApostOfZ	=	ZXGk[2][0];
					printf("\tCoordinateOfX	=	%.20lf\n",AApostOfX);
					printf("\tCoordinateOfY	=	%.20lf\n",AApostOfY);
					printf("\tCoordinateOfZ	=	%.20lf\n",AApostOfZ);
					ResultX		=	AApostOfX / 1000.0;
					ResultY		=	AApostOfY / 1000.0;
					ResultZ		=	AApostOfZ / 1000.0;
					Resulta0	=	a0 * 1000000.0;
					printf("\tResultX		=	%.6lf\n",ResultX);
					printf("\tResultY		=	%.6lf\n",ResultY);
					printf("\tResultZ		=	%.6lf\n",ResultZ);
					printf("\tResulta0		=	%.6lf\n",Resulta0);
						/*-----------------------------------------------------------------------------------------------
					将数据写入文件并保存 
					*/
					printf("/-----------------------------------------------------------------------------------------/\n");
					//fprintf(pfileWrite,"%s %d %d %d %d %.2lf %.2lf\t%.6lf\t%.6lf\t%.6lf\t%.6lf\n",time.PRN,time.year,time.month,time.day,
					//time.hour,t,time.second,ResultX,ResultY,ResultZ,Resulta0);
					fprintf(pfileWrite,"%s %.6lf %.6lf %.6lf %.6lf\n",time.PRN,ResultX,ResultY,ResultZ,Resulta0);
					printf("/-----------------------------------------------------------------------------------------/\n");
				}	
				else if(strcmp(time.PRN,"C11")==0||strcmp(time.PRN,"C12")==0||strcmp(time.PRN,"C14")==0||strcmp(time.PRN,"C19")==0||
				strcmp(time.PRN,"C20")==0||strcmp(time.PRN,"C21")==0||strcmp(time.PRN,"C22")==0||strcmp(time.PRN,"C23")==0||
				strcmp(time.PRN,"C24")==0||strcmp(time.PRN,"C25")==0||strcmp(time.PRN,"C26")==0||strcmp(time.PRN,"C27")==0||
				strcmp(time.PRN,"C28")==0||strcmp(time.PRN,"C29")==0||strcmp(time.PRN,"C30")==0||strcmp(time.PRN,"C32")==0||
				strcmp(time.PRN,"C33")==0||strcmp(time.PRN,"C34")==0||strcmp(time.PRN,"C35")==0||strcmp(time.PRN,"C36")==0||
				strcmp(time.PRN,"C37")==0||strcmp(time.PRN,"C41")==0||strcmp(time.PRN,"C42")==0||strcmp(time.PRN,"C43")==0||
				strcmp(time.PRN,"C44")==0||strcmp(time.PRN,"C45")==0||strcmp(time.PRN,"C46")==0||strcmp(time.PRN,"C06")==0||
				strcmp(time.PRN,"C07")==0||strcmp(time.PRN,"C08")==0||strcmp(time.PRN,"C09")==0||strcmp(time.PRN,"C10")==0||
				strcmp(time.PRN,"C13")==0||strcmp(time.PRN,"C16")==0||strcmp(time.PRN,"C38")==0||strcmp(time.PRN,"C39")==0||
				strcmp(time.PRN,"C40")==0){
					/*---------------------------------------------
					十、计算观测时刻的升交点经度 
					*/
					double	OMEGA_k;
					OMEGA_k	=	OMEGA	+	(OMEGADOT	-	Omega_e) * tk	-	Omega_e * toe;
					printf("\t升交点经度Ωk = %.20lf\n",OMEGA_k);
					printf("/-----------------------------------------------------------------------------------------/\n");
					/*---------------------------------------------
					十一、计算CGCS2000的坐标 	apostrophe--->撇 	coordinate--->坐标 
					*/
					double	ApostOfX,ApostOfY,ApostOfZ; 
					double	ResultX,ResultY,ResultZ,Resulta0;
					ApostOfX	=	Xk * cos(OMEGA_k)	-	Yk * cos(ik) * sin(OMEGA_k);
					ApostOfY	=	Xk * sin(OMEGA_k)	+	Yk * cos(ik) * cos(OMEGA_k);
					ApostOfZ	=	Yk * sin(ik);
					printf("\tCoordinateOfX	=	%.20lf\n",ApostOfX);
					printf("\tCoordinateOfY	=	%.20lf\n",ApostOfY);
					printf("\tCoordinateOfZ	=	%.20lf\n",ApostOfZ);
					ResultX		=	ApostOfX / 1000.0;
					ResultY		=	ApostOfY / 1000.0;
					ResultZ		=	ApostOfZ / 1000.0;
					Resulta0	=	a0 * 1000000.0;
					printf("\tResultX		=	%.6lf\n",ResultX);
					printf("\tResultY		=	%.6lf\n",ResultY);
					printf("\tResultZ		=	%.6lf\n",ResultZ);
					printf("\tResulta0	=	%.6lf\n",Resulta0);
					/*-----------------------------------------------------------------------------------------------
					将数据写入文件并保存 
					*/
					printf("/-----------------------------------------------------------------------------------------/\n");
					fprintf(pfileWrite,"%s %.6lf %.6lf %.6lf %.6lf\n",time.PRN,ResultX,ResultY,ResultZ,Resulta0);
					printf("/-----------------------------------------------------------------------------------------/\n");
				//continue; 
			}
		}
		else{
			continue;
		}
	}
	fclose(pfile);
	fclose(pfileWrite);
	return 0;
}

二、利用MATLAB，并根据18参数解算出来的卫星位置文件绘制三维图

一、代码

MATLAB读文件的位置直接修改为18参数写文件的路径，这样在C程序中直接修改卫星型号，MATLAB就能实时绘图。

%**************************************** Message ***********************************//
%//技术交流：folkore0118@163.com
%//关注优快云博主：“彼稷”
%//Author: 彼稷
%//All rights reserved	                               
%//----------------------------------------------------------------------------------------
%// Target Devices: 		DELL-G15
%// Tool Versions:       MATLAB R2021a
%// File name:           draw
%// Last modified Date:  2023年10月18日14:00:00
%// Last Version:        
%// Descriptions:        卫星位置&三维图
%//----------------------------------------------------------------------------------------

% 文件路径
filename = '18参数写文件的位置';

% 使用 importdata 函数读取文件
data = importdata(filename, ' ');

% 获取名称、X坐标、Y坐标和Z坐标
name = data.textdata;
x = data.data(:, 1);
y = data.data(:, 2);
z = data.data(:, 3);

% 绘制三维图
figure;
scatter3(x, y, z);
hold on;

% 添加球面
center = [0, 0, 0]; % 球心坐标
radius = 6371; % 球半径
[X, Y, Z] = sphere;
X = radius * X;
Y = radius * Y;
Z = radius * Z ;
surf(X, Y, Z, 'FaceColor', 'none', 'EdgeColor', 'red');
axis equal;
grid on;
xlabel('X');
ylabel('Y');
zlabel('Z');
title('SUM坐标图');

二、三维图

外围蓝色空心的是北斗卫星，中间红色的球面是地球。

根据图中卫星位置来看，在南极和北极上空均没有卫星分布；且IGSO卫星呈“8”字绕行，符合客观事实。

MEO卫星分布：

bd10bbbb0bd3408c9323cd95b80f4311.png

IGSO卫星分布：

52177ac2ccfa4811be045b381d12b515.png

GEO卫星分布：

c88be47b676342b384c3d58b901e444a.png

所有北斗卫星分布：

51cc637288e6478d9f689ad227ec7a79.png

 

三、最后

本文只是记录学习过程中的一些感想并复盘每阶段的任务，代码仍有优化的空间、内容有谬误之处请大家多多指正。

 

 

 

 

 

 

 

 

 

 

用广播星历计算卫星位置
04-09
用C++语言实现利用广播星历计算卫星位置
利用北斗BDS广播星历，进行卫星位置计算
qq_45568498的博客
12-19 9601
利用RINEX文件中的星历和卫星发射时刻计算北斗MEO/IGSO卫星和GEO卫星的位置
细说卫星导航：测距定位原理
最新发布
极客不孤独的博客
03-22 1516
本文介绍了卫星导航的测距定位原理，通过两个演示实验进行详细介绍。
基于广播星历的北斗定位解算原理（基于C语言和MATLAB实现）_卫星位置解算
2401_84978510的博客
05-16 1052
printf(“\t摄动改正量δu = %.20lf\n”,delta_u);printf(“\t摄动改正量δr = %.20lf\n”,delta_r);printf(“\t摄动改正量δi = %.20lf\n”,delta_i);printf(“\t摄动改正后的升交距角uk = %.20lf\n”,uk);printf(“\t摄动改正后的卫星矢径rk = %.20lf\n”,rk);printf(“\t输入参考时间t:(小于60分钟)\t====>”);
（武汉大学GPS作业）根据广播星历文件计算卫星位置
10-19
武汉大学GPS原理与应用课程结课作业，根据广播星历文件计算卫星位置，给有需要的小伙伴们参考^_^！ （项目工程环境Visual Studio 2010，在release文件夹中可找到独立运行的exe程序）
c 计算delta_生活丨由广播星历(导航电文)简单计算卫星位置
weixin_42358173的博客
11-29 1309
//写这个花了一天的时间 主要是那一块费了很久//于是专门把时间背景知识加进来做一个总结//经常感觉很口渴 思路也很闷 就像上面的停了水的喷灌器一、时间背景知识二、不同时间系统之间的转换利用导航电文计算卫星位置时 需要将目标计算时刻换算到GPS时间//由于导航电文中的是GPS周内的秒数//进行时间转换时：需要先将(通常是格里高利历)转换到，再将<儒略历下时刻>转换...
基于C语言 读取精密星历和广播星历并计算卫星坐标
05-12
本项目以C语言为工具，实现了读取精密星历和广播星历，并对卫星坐标进行计算的功能，这在GNSS信号处理、误差分析以及精度评估等方面具有重要意义。 首先，精密星历与广播星历是两种不同级别的星历数据。广播星历是...
基于不同星历对GPS解算结果的影响分析.pdf
09-08
本文深入探讨了广播星历和精密星历（包括最终星历、快速星历、超快星历）的精度差异，并通过实测数据解算分析了它们对GPS定位结果的影响。 首先，广播星历是由GPS地面控制系统提供并广播给全球用户的预估值，精度较...
matlab GPS卫星定位 利用广播星历计算卫星位置
09-28
5. MATLAB实现细节：在MATLAB中，可以编写脚本来读取和解析广播星历文件，然后用这些数据计算卫星的即时位置。利用“for”循环处理每一颗卫星的信息，结合伪距和多普勒频移修正，最终通过非线性优化算法（如...
广播星历和精密星历对解算某港区GPS控制网成果影响分析.pdf
09-08
本文主要探讨了广播星历和精密星历在GPS控制网解算中的影响。作者杜志彪通过对比两种星历的差异，利用某港区的GPS控制网外业观测数据，分析了它们对基线解算和网平差结果的差异。研究表明，在常规的短基线GPS控制网...
广播星历解算卫星位置
11-22
利用广播星历解算卫星位置。需先将星历文件拷贝至txt文档
GPS-广播星历计算卫星位置
09-13
本程序用C语言进行编程，提高了We和E的精度，实现了以下功能： (1)实现星历数据文件的导入并存储到数组中； (2)实现三个时刻坐标的同时计算并与验算数据比对； (3)实现了计算结果的屏幕输出和result文件输出； 计算结果精确，与提供的验算数据比对小数点后三位基本相同。
北斗四星定位解算程序
12-05
基于伪距的北斗导航的四星定位解算，并将解算结果转换为经纬度
C++实现由GPS广播星历计算卫星位置
05-08
武汉大学测绘学院GNSS原理及应用编程作业。包含数据文件，结果文件。
读取广播星历计算卫星位置
01-20
代码是在VS2010环境下编写的，在主函数中输入卫星的编号及年月日时分秒，获得卫星的位置信息
GNSS-基于广播星历计算卫星位置
m0_58307078的博客
12-20 3497
本文内容主要包括基于广播星历卫星位置计算原理、所需文件及其格式内容、数据结构、计算步骤、编程流程、算例及结果分析的阐述，以及笔者对编程过程中一些常见问题和注意事项的总结。
北斗信号服务器解算,北斗导航系统接收机定位解算设计与实现
weixin_39686230的博客
08-09 926
摘要：随着北斗导航系统的建设不断推进,其应用范围越来越广,因此北斗接收机需求也越来越大。不同的应用场景的接收机结构和侧重点有所不同,但是其中的定位解算模块都是其关键部分。本文主要对北斗接收机的整体结构和定位解算模块进行设计与开发。论文首先介绍了北斗系统的组成、信号结构和接收机定位原理,在此基础上研究了接收机伪距定位中常用的最小二乘解算方法和基于扩展卡尔曼滤波器的定位解算方法,分析了两种算法在静态动...
利用广播星历解算北斗卫星位置及精度分析
热门推荐
weixin_48630797的博客
09-28 1万+
网上已经有了比较多关于如何利用播发的广播星历来解算卫星位置的blog，此文章的目的是为了记录一下学习成果的同时回馈一下（毕竟看过很多blog了，但是从来没写过）。
评论 11
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

彼稷

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付：¥1
获取中
扫码支付

您的余额不足，请更换扫码支付或充值

打赏作者

实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值