基于广播星历的北斗卫星轨道计算（16参数和18参数）

已于 2025-03-17 16:55:56 修改
阅读量3.9k 收藏 57
点赞数 11
分类专栏： 卫星定位 文章标签： 算法 开发语言 c语言
于 2023-10-20 12:41:35 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/m0_62942334/article/details/133935512
版权

最近在用C语言写基于广播星历的北斗卫星轨道计算，采用轨道根数为16和18两种方式，特在此记下，方便之后学习的同学。

一、16参数卫星轨道计算

我们首先来看一下计算公式。

 分为MEO/IGSO卫星和GEO卫星

 我使用的数据文件的名字为 BRDM00DLR_S_20231520000_01D_MN.rnx

16参数位置如下，为RINEX3.04标准（广播星历格式），按下表读取我们需要的参数。

采用的Dev C++编写的C代码，粘贴在下方

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<string.h>
#define MAX_DATA_SIZE 300000

const double GM = 3.986005e+14,we = 7.292115e-05,pi = 3.1415926;

typedef struct
{
	int year;  //年
	int month; //月
	int day;   //日
	int hour;  //时
	int minute;//分
	double second;//秒
}Time;

typedef struct
{
	char name[10];
}Parameter;

double UTC2BDST(int year, int month, int day, int hour, int minute, double second, double *secondOfweek)
{
	double JD, WN;
	if (month <= 2)
	{
    	year -= 1;
    	month += 12;
	}
	JD = 365.25 * year + (int)(30.6001 * (month + 1)) + day + 1720981.5 + hour / 24.0 + minute / 1440 + second / 86400 ;  // 计算儒略日
	WN = (int)((JD - 2453736.5) / 7) ; // WN:GPS_week number 目标时刻的GPS周
	*secondOfweek = (int(JD - 2453736.5) - (7.0 * WN)) * 24 * 3600.0 ;  // tGPS:目标时刻的GPS秒 减去14秒为BDT
}

double UTC2BDS(int year, int month, int day, int hour, int minute, double second, double *secondOfweek)
{ 
   int DayofYear = 0;
   int DayofMonth = 0;
   for (int i = 2006; i < year; i++)  
   {
      if ((i % 4 == 0 && i % 100 != 0) || i % 400 == 0)
                DayofYear += 366;
      else
               DayofYear += 365;
   }
   for (int i = 1; i < month; i++)   
   {
         if (i == 1 || i == 3 || i == 5 || i == 7 || i == 8 || i == 10 || i ==12)
                DayofMonth += 31;
          else if (i == 4 || i == 6 || i == 9 || i == 11)
                DayofMonth += 30;
          else
          {
            if ((year % 4 == 0 && year % 100 != 0) || year % 400 == 0)
                DayofMonth += 29;
            else
                DayofMonth += 28;
          }
    }
   int Day;
   Day = DayofMonth + day + DayofYear-1;//2006年1月1日起始历元
   *secondOfweek = Day % 7 * 86400 + hour * 3600 + minute * 60 + second;

   return *secondOfweek;
}

void FigureOut(FILE *file,FILE *FileWrite)
{
	int data_count = 0;
	int EphemerisData[MAX_DATA_SIZE];
	char ReadPath[50] = "E:\\Buffer\\data_16.txt";
	char WritePath[50] = "E:\\Buffer\\result_16.txt";
	file = fopen(ReadPath,"r");
	FileWrite = fopen(WritePath,"w");
	if(!file||!FileWrite)
	{
		printf("open the file error!\n");
	}
	
	Time time;//年月日时分秒
	Parameter para;
    double a0, a1, a2;//钟差 钟速 钟漂 
	double IODE, Crs, Delta_n, M0;//星历数据龄期  轨道半径的正弦调和项 平均角速度摄动量 Toe时刻的平近点角
	double Cuc, e, Cus, sqrt_a;//纬度幅角的余弦调和项 轨道偏心率 纬度幅角的正弦调和项 轨道长半径的平方根
	double toe, Cic, OMEGA_0, Cis;//星历表参考历元 轨道倾角的余弦调和项 升交点赤经 轨道倾角的正弦调和项
	double i0, Crc, w, OMEGA_DOT;//轨道倾角 轨道半径的余弦调和项 近地点角距 升交点赤经变化率 
	double IDOT, cflgl2, weekno, pflgl2;//轨道倾角变化率 电文来源 BDS的整周计数 长半轴变化率
	double svacc, svhlth, tgd, IODC;//卫星精度 卫星自主健康标识 TGD1星上设备时延差 TGD2星上设备时延差
	double xxx1,xxx2;
	
	char name[5];
	double secondOfweek1, secondOfweek2, secondOfweek;
	double n, n0, t_observe, t, tk , Delta_t, toc, Mk, Ek, Vk, fai_k;
	double _6u, _6r, _6i, uk, rk, ik, Xk, Yk, E0, E1;
	double OMEGA_k, X_obs_k, Y_obs_k, Z_obs_k, X__obs_k, Y__obs_k, Z__obs_k;
	
		while(1)
		{	
			if(feof(file))
			{
				printf("文件读取结束");
				break;
			}
		
		fscanf(file,"%s",&name);
		//printf("%s\n",name);
		
		fscanf(file, " %d %d %d %d %d %lf %lf %lf %lf",&time.year,&time.month,
		&time.day,&time.hour,&time.minute,&time.second,&a0, &a1, &a2);
		//printf(" %d %d %d %d %d %.20lf %.20lf %.20lf %.20lf \n", time.year, time.month,
		//time.day, time.hour, time.minute, time.second, a0, a1, a2);
			
		fscanf(file, " %lf %lf %lf %lf ", &IODE,&Crs,&Delta_n,&M0);
		//printf(" %.10lf \t %.10lf \t %.10lf \t %.10lf \n", IODE, Crs, Delta_n, M0);

		fscanf(file, " %lf %lf %lf %lf ", &Cuc, &e, &Cus, &sqrt_a);
		//printf(" %.10lf \t %.10lf \t %.10lf \t %.10lf \n", Cuc, e, Cus, sqrt_a);

		fscanf(file, " %lf %lf %lf %lf ", &toe, &Cic, &OMEGA_0, &Cis);
		//printf(" %.10lf \t %.10lf \t %.10lf \t %.10lf \n", toe, Cic, OMEGA_0, Cis);

		fscanf(file, " %lf %lf %lf %lf ", &i0, &Crc, &w, &OMEGA_DOT);
		//printf(" %.10lf \t %.10lf \t %.10lf \t %.10lf \n", i0, Crc, w, OMEGA_DOT);

		fscanf(file, " %lf %lf %lf %lf ", &IDOT, &cflgl2, &weekno, &pflgl2);
		//printf(" %.20lf \t %.20lf \t %.20lf \t %.20lf \n", IDOT, cflgl2, weekno, pflgl2);
			
		fscanf(file, " %lf %lf %lf %lf ", &svacc, &svhlth, &tgd, &IODC);
		//printf(" %.10lf \t %.10lf \t %.10lf \t %.10lf \n", svacc, svhlth, tgd, IODC);
		
		fscanf(file, " %lf %lf ", &xxx1, &xxx2);
		//printf(" %.10lf \t %.10lf \n", xxx1,xxx2);
					
		/*计算卫星运行的平均角速度n*/
		n0 = sqrt(GM)/pow(sqrt_a,3);
		n = n0 + Delta_n;
		
		/*计算归化时间t_k*/
		t_observe = UTC2BDST(time.year,time.month,time.day,time.hour,time.minute,time.second,&secondOfweek1); 
		toc = UTC2BDST(time.year,time.month,time.day,time.hour,time.minute,time.second,&secondOfweek2);
		Delta_t = a0 + a1*(t_observe-toc) + a2*(t_observe-toc)*(t_observe-toc);
		t = t_observe - Delta_t;
		//t = UTC2BDS(time.year,time.month,time.day,time.hour,time.minute ,time.second,&secondOfweek);
		tk = t - toe - 14;
		if(tk>302400)
			tk -= 604800;
		else
		if(tk<-302400)
			tk += 604800;
		printf("归化时间tk = %.20lf\n",tk);
		
		/*计算观测时刻卫星平近点角Mk*/
		Mk = M0 + n*tk;
		printf("平近点角Mk = %.20lf\n",Mk);
		
		/*计算偏近点角Ek*/
		  Ek = Mk;
		  Ek = Mk + e*sin(Ek);
		  for(int i=0;i<2;i++)
		 {
			Ek = Mk + e*sin(Ek);
		 }
		
		//Ek=Mk;
		//E1=Mk+e*sin(Ek);
		//while ( abs(E1-Ek)>1e-10 )
		//{
		//	Ek=E1;
    	//	E1=Mk+e*sin(E1);	
		//}
		//Ek=E1;
		printf("偏近点角Ek = %.20lf\n",Ek);

		/*计算真近点角Vk*/
		//Vk = atan(sqrt(1-e*e)*sin(Ek)/(cos(Ek)-e));
		//Vk=2*atan2(sqrt((1+e)/(1-e))*tan(Ek/2),1);
		Vk = atan2(sqrt(1-e*e)*sin(Ek),(cos(Ek)-e));
		printf("真近点角Vk = %.20lf\n",Vk);
		
		/*计算升交距角fai_k*/
		fai_k = Vk + w;
		printf("升交距角fai_k = %.20lf\n",fai_k);
		
		/*计算摄动改正量*/
		_6u = Cuc*cos(2*fai_k) + Cus*sin(2*fai_k);
		_6r = Crc*cos(2*fai_k) + Crs*sin(2*fai_k);
		_6i = Cic*cos(2*fai_k) + Cis*sin(2*fai_k);
		printf("升摄动改正量 _6u = %.20lf\n",_6u);
		printf("升摄动改正量 _6r = %.20lf\n",_6r);
		printf("升摄动改正量 _6i = %.20lf\n",_6i);
		
		/*计算摄动改正后的 升交距角uk 卫星矢径rk 轨道倾角ik*/
		uk = fai_k + _6u;
		rk = pow(sqrt_a,2)*(1-e*cos(Ek)) + _6r;
		ik = i0 + _6i + IDOT*tk;
		printf("升交距角uk = %.20lf\n",uk);
		printf("卫星矢径rk = %.20lf\n",rk);
		printf("轨道倾角ik = %.20lf\n",ik);
		
		/*计算卫星在轨道平面坐标系的坐标*/
		Xk = rk*cos(uk);
		Yk = rk*sin(uk); 
		printf("Xk = %.20lf\n",Xk);
		printf("Yk = %.20lf\n",Yk);
	
		if(strcmp(name,"C01")==0||strcmp(name,"C02")==0||strcmp(name,"C03")==0||strcmp(name,"C04")==0||strcmp(name,"C05")==0)
		{
			/*GEO*/
			OMEGA_k = OMEGA_0 + OMEGA_DOT*tk - we*toe;
			printf("升交点经度OMEGA_k = %.20lf\n",OMEGA_k);			
		
			X__obs_k = Xk*cos(OMEGA_k)-Yk*cos(ik)*sin(OMEGA_k);
			Y__obs_k = Xk*sin(OMEGA_k)+Yk*cos(ik)*cos(OMEGA_k);
			Z__obs_k = Yk*sin(ik);
		
			X_obs_k = cos(we*tk)*X__obs_k + sin(we*tk)*cos(-5.0/180*pi)*Y__obs_k + sin(we*tk)*sin(-5.0/180*pi)*Z__obs_k;
			Y_obs_k = -1*sin(we*tk)*X__obs_k + cos(we*tk)*cos(-5.0/180*pi)*Y__obs_k + cos(we*tk)*sin(-5.0/180*pi)*Z__obs_k;
			Z_obs_k = -1*sin(-5.0/180*pi)*Y__obs_k + cos(-5.0/180*pi)*Z__obs_k;
					
		}
		else	
		{
			/*MEO/IGSO*/
			/*计算观测时刻升交点经度OMEGA_k*/
			OMEGA_k = OMEGA_0 + (OMEGA_DOT-we)*tk - we*toe;
			printf("升交点经度OMEGA_k = %.20lf\n",OMEGA_k);
	
			/*计算卫星在CGCS2000中的坐标*/
			X_obs_k = Xk*cos(OMEGA_k)-Yk*cos(ik)*sin(OMEGA_k);
			Y_obs_k = Xk*sin(OMEGA_k)+Yk*cos(ik)*cos(OMEGA_k);
			Z_obs_k = Yk*sin(ik);
		}
		printf("%s\n",name);
		printf("X_obs_k=%.6lf Y_obs_k=%.6lf Z_obs_k=%.6lf a0=%.6lf\n\n",X_obs_k/1000.0,Y_obs_k/1000.0,Z_obs_k/1000.0,a0*1e+06);
		
		fprintf(FileWrite,"%s %.6lf %.6lf %.6lf %.6lf\n",name,X_obs_k/1000.0,Y_obs_k/1000.0,Z_obs_k/1000.0,a0*1e+06);
		
	}
	
	fclose(file);
	fclose(FileWrite);
}

int main()
{
	FILE *file;
	FILE *FileWrite;
	
	FigureOut(file,FileWrite);
	
	return 0;
}

 该解算有三个需要注意的地方。

第一，计算周内秒 t 有两个方法。一是程序中的UTC2BDS函数，用参考时刻的年月日减去BDST的起始历元2006年1月1日0时0分0秒，从而得出周内秒。二是使用儒略日计算，见程序中的UTC2BDST函数。因为解算方法是在GPS的基础上，BDST = GPST -14s，所以 tk = t - toe -14。

第二，计算偏近点角时使用正确的迭代方法，我在网上找了多种迭代方法，经过尝试发现代码中这种方法最好。

第三，计算真近点角时用的是arctan2。因为arctan的取值范围是-pi/2到pi/2，不满足真近点角0到2pi的取值范围。

关于卫星类型，在 星座状态 (csno-tarc.cn) 网站查询卫星信息，看是GEO，MEO还是IGSO。

二、18参数卫星轨道计算

2.1国内IGMAS RINEX3.03标准

首先我们先来看计算公式

18参数RINEX3.03位置如图所示

 C语言代码如下。

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<string.h> 
#define MAX_DATA_SIZE 300000

const double GM = 3.986005e+14, we = 7.292115e-05,
			 pi = 3.1415926, miu = 3.986004418e+14,
			 Aref_MEO = 27906100, Aref_IGSOGEO = 42162200;

typedef struct
{
	int year;  //年
	int month; //月
	int day;   //日
	int hour;  //时
	int minute;//分
	double second;//秒
}Time;

typedef struct
{
	char name[10];
    double a0, a1, a2;//钟差 钟速 钟漂 
	double IODE, Crs, Delta_n0, M0;//星历参数版本号  轨道半径的正弦调和改正项的振幅 参考时刻卫星平均角速度与计算值之差 Toe时刻的平近点角
	double Cuc, e, Cus, Delta_A;//纬度幅角的余弦调和改正项的振幅 轨道偏心率 纬度幅角的正弦调和改正项的振幅 轨道长半径的平方根
	double toe, Cic, OMEGA_0, Cis;//星历表参考历元 轨道倾角的余弦调和项 升交点赤经 轨道倾角的正弦调和项
	double i0, Crc, w, OMEGA_DOT;//轨道倾角 轨道半径的余弦调和项 近地点角距 升交点赤经变化率 
	double IDOT, Data, weekno, A_DOT;//轨道倾角变化率 电文来源 BDS的整周计数 长半轴变化率
	double svacc, HS_SatH1, TGD_TGD1, ISC_TGD2;//卫星精度 卫星自主健康标识 TGD1星上设备时延差 TGD2星上设备时延差
	double BDS_Second, IODC_AODC, Delta_n0_DOT, SatType;//信息的发射时间（ BDT 的周内秒）钟差参数版本号/时钟数据龄期
														//参考时刻卫星平均角速度与值之差的变化率  留空 卫星轨道类型
	double secondOfweek1, secondOfweek2, secondOfweek;
	double n, n0, t_observe, t, tk , Delta_t, toc, Mk, Vk, phi_k, A0, Ak, Delta_nA, nA;
	double _6u, _6r, _6i, uk, rk, ik, Xk, Yk, Ek, E1;                                       
	double OMEGA_k, X_obs_k, Y_obs_k, Z_obs_k, X__obs_k, Y__obs_k, Z__obs_k;	
}Parameter;

double UTC2BDST(int year, int month, int day, int hour, int minute, double second, double *secondOfweek)
{
	double JD, WN;
	if (month <= 2)
	{
    	year -= 1;
    	month += 12;
	}
	JD = 365.25 * year + int(30.6001 * (month + 1)) + day + 1720981.5 + hour / 24.0 + minute / 1440 + second / 86400 ;  // 计算儒略日
	WN = int((JD - 2453736.5) / 7) ; // WN:GPS_week number 目标时刻的GPS周
	*secondOfweek = (int(JD - 2453736.5) - (7.0 * WN)) * 24 * 3600.0 ;  // tGPS:目标时刻的GPS秒 减去14秒为BDT

    return *secondOfweek;
}

double UTC2BDS(int year, int month, int day, int hour, int minute, double second, double *secondOfweek)
{ 
   int DayofYear = 0;
   int DayofMonth = 0;
   for (int i = 2006; i < year; i++)  
   {
      if ((i % 4 == 0 && i % 100 != 0) || i % 400 == 0)
                DayofYear += 366;
      else
               DayofYear += 365;
   }
   for (int i = 1; i < month; i++)   
   {
         if (i == 1 || i == 3 || i == 5 || i == 7 || i == 8 || i == 10 || i ==12)
                DayofMonth += 31;
          else if (i == 4 || i == 6 || i == 9 || i == 11)
                DayofMonth += 30;
          else
          {
            if ((year % 4 == 0 && year % 100 != 0) || year % 400 == 0)
                DayofMonth += 29;
            else
                DayofMonth += 28;
          }
    }
   int Day;
   Day = DayofMonth + day + DayofYear-1;//2006年1月1日起始历元
   *secondOfweek = Day % 7 * 86400 + hour * 3600 + minute * 60 + second;

   return *secondOfweek;
}

void FigureOut(FILE *file,FILE *FileWrite)
{
	int data_count = 0;
	int EphemerisData[MAX_DATA_SIZE];
	file = fopen("E:\\Buffer\\data18_IGMAS.txt","r");
	FileWrite = fopen("E:\\Buffer\\result18_IGMAS.txt","w");
	if(!file||!FileWrite)
	{
		printf("open the file error!\n");
	}
	
	Time time;//年月日时分秒
	Parameter para;
	
		while(1) 
		{	
			if(feof(file))
			{
				printf("文件读取结束");
				break;
			}
		
		fscanf(file,"%s",&para.name);
		//printf("%s\n",name);
		
		fscanf(file, " %d %d %d %d %d %lf %lf %lf %lf",&time.year, &time.month, &time.day, &time.hour, &time.minute, &time.second, &para.a0, &para.a1, &para.a2);
		//printf(" %d %d %d %d %d %.20lf %.20lf %.20lf %.20lf \n", time.year, time.month,time.day, time.hour, time.minute, time.second, para.a0, para.a1, para.a2);
			
		if( (strcmp(para.name,"C19")==0||strcmp(para.name,"C20")==0||
		strcmp(para.name,"C21")==0||strcmp(para.name,"C22")==0||strcmp(para.name,"C23")==0||strcmp(para.name,"C24")==0||
		strcmp(para.name,"C25")==0||strcmp(para.name,"C26")==0||strcmp(para.name,"C27")==0||strcmp(para.name,"C28")==0||
		strcmp(para.name,"C29")==0||strcmp(para.name,"C30")==0||strcmp(para.name,"C32")==0||strcmp(para.name,"C33")==0||
		strcmp(para.name,"C34")==0||strcmp(para.name,"C35")==0||strcmp(para.name,"C36")==0||strcmp(para.name,"C37")==0||
		strcmp(para.name,"C38")==0||strcmp(para.name,"C39")==0||strcmp(para.name,"C40")==0||strcmp(para.name,"C41")==0||
		strcmp(para.name,"C42")==0||strcmp(para.name,"C43")==0||strcmp(para.name,"C44")==0||strcmp(para.name,"C45")==0||
		strcmp(para.name,"C46")==0)&&(time.year==2023) )
        //if( (strcmp(para.name,"C06")==0||strcmp(para.name,"C07")==0||strcmp(para.name,"C08")==0||strcmp(para.name,"C09")==0||
        //strcmp(para.name,"C10")==0||strcmp(para.name,"C13")==0||strcmp(para.name,"C16")==0||strcmp(para.name,"C31")==0||
        //strcmp(para.name,"C38")==0||strcmp(para.name,"C39")==0||strcmp(para.name,"C40")==0||strcmp(para.name,"C56")==0)&&(time.year==2023) )
		{
		fscanf(file, " %lf %lf %lf %lf ", &para.IODE,&para.Crs,&para.Delta_n0,&para.M0);
		//printf(" %.10lf \t %.10lf \t %.10lf \t %.10lf \n", para.IODE, para.Crs, para.Delta_n0, para.M0);

		fscanf(file, " %lf %lf %lf %lf ", &para.Cuc, &para.e, &para.Cus, &para.Delta_A);
		//printf(" %.10lf \t %.10lf \t %.10lf \t %.10lf \n", para.Cuc, para.e, para.Cus, para.Delta_A);

		fscanf(file, " %lf %lf %lf %lf ", &para.toe, &para.Cic, &para.OMEGA_0, &para.Cis);
		//printf(" %.10lf \t %.10lf \t %.10lf \t %.10lf \n", para.toe, para.Cic, para.OMEGA_0, para.Cis);

		fscanf(file, " %lf %lf %lf %lf ", &para.i0, &para.Crc, &para.w, &para.OMEGA_DOT);
		//printf(" %.10lf \t %.10lf \t %.10lf \t %.10lf \n", para.i0, para.Crc, para.w, para.OMEGA_DOT);

		fscanf(file, " %lf %lf %lf %lf ", &para.IDOT, &para.Data, &para.weekno, &para.A_DOT);
		//printf(" %.20lf \t %.20lf \t %.20lf \t %.20lf \n", para.IDOT, para.Data, para.weekno, para.A_DOT);
			
		fscanf(file, " %lf %lf %lf %lf ", &para.svacc, &para.HS_SatH1, &para.TGD_TGD1, &para.ISC_TGD2);
		//printf(" %.10lf \t %.10lf \t %.10lf \t %.10lf \n", para.svacc, para.HS_SatH1, para.IGD_TGD1, para.ISC_TGD2);
		
		fscanf(file, " %lf %lf %lf %lf ", &para.BDS_Second, &para.IODC_AODC, &para.Delta_n0_DOT, &para.SatType);
		//printf(" %.10lf \t %.10lf \n", para.BDS_Second, para.IODC_AODC, para.Delta_n0_DOT, para.SatType);
		
		/*计算归化时间t_k*/
		//t_observe = UTC2BDST(time.year,time.month,time.day,time.hour,time.minute,time.second,&para.secondOfweek1); 
		//toc = UTC2BDST(time.year,time.month,time.day,time.hour,time.minute,time.second,&para.secondOfweek2);
		//para.Delta_t = para.a0 + para.a1*(para.t_observe-para.toc) + para.a2*(para.t_observe-para.toc)*(para.t_observe-para.toc);
		//t = t_observe - Delta_t;
		para.t = UTC2BDS(time.year,time.month,time.day,time.hour,time.minute ,time.second,&para.secondOfweek);
		para.tk = para.t - para.toe - 14;
		if(para.tk>302400)
        {
            para.tk -= 604800;
        }
		else
		if(para.tk<-302400)
        {
            para.tk += 604800;
        }

		printf("归化时间tk = %.20lf\n",para.tk);

        /*计算参考时刻的长半轴*/
        if(strcmp(para.name,"C11")==0||strcmp(para.name,"C12")==0||strcmp(para.name,"C14")==0||strcmp(para.name,"C19")==0||strcmp(para.name,"C20")==0||strcmp(para.name,"C21")==0
        ||strcmp(para.name,"C22")==0||strcmp(para.name,"C23")==0||strcmp(para.name,"C24")==0||strcmp(para.name,"C25")==0||strcmp(para.name,"C26")==0||strcmp(para.name,"C27")==0
        ||strcmp(para.name,"C28")==0||strcmp(para.name,"C29")==0||strcmp(para.name,"C30")==0||strcmp(para.name,"C32")==0||strcmp(para.name,"C33")==0||strcmp(para.name,"C34")==0
        ||strcmp(para.name,"C35")==0||strcmp(para.name,"C36")==0||strcmp(para.name,"C37")==0||strcmp(para.name,"C41")==0||strcmp(para.name,"C42")==0||strcmp(para.name,"C43")==0
        ||strcmp(para.name,"C44")==0||strcmp(para.name,"C45")==0||strcmp(para.name,"C46")==0||strcmp(para.name,"C57")==0||strcmp(para.name,"C58")==0)
        {para.A0 = Aref_MEO +  para.Delta_A;}
        else
        {para.A0 = Aref_IGSOGEO +  para.Delta_A;}
		printf("参考时刻的长半轴A0 = %.20lf\n",para.A0);
		 
		/*计算长半轴Ak*/
		para.Ak = para.A0 +  para.A_DOT*para.tk;
		printf("长半轴Ak = %.20lf\n",para.Ak);
		
		/*计算参考时刻的卫星平均角速度n0*/
		para.n0 = sqrt(miu/pow(para.A0,3));
		printf("参考时刻的卫星平均角速度n0 = %.20lf\n",para.n0);
		
		/*计算卫星平均角速度的偏差Delta_nA*/
		para.Delta_nA = para.Delta_n0 + (1.0/2)*para.Delta_n0_DOT*para.tk;
		printf("卫星平均角速度的偏差Delta_nA = %.20lf\n",para.Delta_nA);
		
		/*计算改正后的卫星平均角速度nA*/
		para.nA = para.n0 + para.Delta_nA;
		printf("改正后的卫星平均角速度nA = %.20lf\n",para.nA);
		
		/*计算平近点角Mk*/
		para.Mk = para.M0 + para.nA*para.tk;
		printf("平近点角Mk = %.20lf\n",para.Mk);
		
		/*计算偏近点角Ek*/
		para.Ek = para.Mk;
		para.Ek = para.Mk + para.e*sin(para.Ek);
		for(int i=0;i<2;i++)
		{
			para.Ek = para.Mk + para.e*sin(para.Ek);
		}
		
		//int count = 0;
		//para.Ek=para.Mk;
		//para.E1=para.Mk+para.e*sin(para.Ek);
		//while ( abs(para.E1-para.Ek)>1e-10 )
		//{
		//	para.Ek=para.E1;
    	//	para.E1=para.Mk+para.e*sin(para.E1);
    	//	count++;
		//	if(count>1e+06)
		//	{
		//		break;	
		//	}	
		//}
		//para.Ek=para.E1;
		printf("偏近点角Ek = %.20lf\n",para.Ek);

		/*计算真近点角Vk*/
		//para.Vk = atan(sqrt(1-para.e*para.e)*sin(para.Ek)/(cos(para.Ek)-para.e));
		//para.Vk=2*atan2(sqrt((1+para.e)/(1-para.e))*tan(para.Ek/2),1);
		para.Vk = atan2(sqrt(1-para.e*para.e)*sin(para.Ek),(cos(para.Ek)-para.e));
		printf("真近点角Vk = %.20lf\n",para.Vk);
		
		/*计算纬度幅角fai_k*/
		para.phi_k = para.Vk + para.w;
		printf("升交距角fai_k = %.20lf\n",para.phi_k);
		
		/*计算摄动改正量*/
		para._6u = para.Cus*sin(2*para.phi_k) + para.Cuc*cos(2*para.phi_k);
		para._6r = para.Crs*sin(2*para.phi_k) + para.Crc*cos(2*para.phi_k);
		para._6i = para.Cis*sin(2*para.phi_k) + para.Cic*cos(2*para.phi_k);
		printf("纬度幅角改正量 _6u = %.20lf\n",para._6u);
		printf("径向距离改正量 _6r = %.20lf\n",para._6r);
		printf("轨道倾角改正量 _6i = %.20lf\n",para._6i);
		
		/*计算摄动改正后的 纬度幅角uk 径向距离rk 轨道倾角ik*/
		para.uk = para.phi_k + para._6u;
		para.rk = para.Ak*(1-para.e*cos(para.Ek)) + para._6r;
		para.ik = para.i0 + para.IDOT*para.tk + para._6i;
		printf("升交距角uk = %.20lf\n",para.uk);
		printf("卫星矢径rk = %.20lf\n",para.rk);
		printf("轨道倾角ik = %.20lf\n",para.ik);
		
		/*计算卫星在轨道平面坐标系的坐标*/
		para.Xk = para.rk*cos(para.uk);
		para.Yk = para.rk*sin(para.uk); 
		printf("Xk = %.20lf\n",para.Xk);
		printf("Yk = %.20lf\n",para.Yk);
	
		if(strcmp(para.name,"C01")==0||strcmp(para.name,"C02")==0||strcmp(para.name,"C03")==0||strcmp(para.name,"C04")==0||strcmp(para.name,"C05")==0)
		{
			/*GEO*/
			para.OMEGA_k = para.OMEGA_0 + para.OMEGA_DOT*para.tk - we*para.toe;
		
			para.X__obs_k = para.Xk*cos(para.OMEGA_k)-para.Yk*cos(para.ik)*sin(para.OMEGA_k);
			para.Y__obs_k = para.Xk*sin(para.OMEGA_k)+para.Yk*cos(para.ik)*cos(para.OMEGA_k);
			para.Z__obs_k = para.Yk*sin(para.ik);
		
			para.X_obs_k = cos(we*para.tk)*para.X__obs_k + sin(we*para.tk)*cos(-5.0/180*pi)*para.Y__obs_k + sin(we*para.tk)*sin(-5.0/180*pi)*para.Z__obs_k;
			para.Y_obs_k = -1*sin(we*para.tk)*para.X__obs_k + cos(we*para.tk)*cos(-5.0/180*pi)*para.Y__obs_k + cos(we*para.tk)*sin(-5.0/180*pi)*para.Z__obs_k;
			para.Z_obs_k = -1*sin(-5.0/180*pi)*para.Y__obs_k + cos(-5.0/180*pi)*para.Z__obs_k;
					
		}
		else	
		{
			/*MEO/IGSO*/
			/*计算改正后的升交点经度OMEGA_k*/
			para.OMEGA_k = para.OMEGA_0 + (para.OMEGA_DOT-we)*para.tk - we*para.toe;
			printf("升交点经度OMEGA_k = %.20lf\n",para.OMEGA_k);
	
			/*计算卫星在CGCS2000中的坐标*/
			para.X_obs_k = para.Xk*cos(para.OMEGA_k)-para.Yk*cos(para.ik)*sin(para.OMEGA_k);
			para.Y_obs_k = para.Xk*sin(para.OMEGA_k)+para.Yk*cos(para.ik)*cos(para.OMEGA_k);
			para.Z_obs_k = para.Yk*sin(para.ik);
		}
		printf("%s\n",para.name);
		printf("X_obs_k=%.6lf Y_obs_k=%.6lf Z_obs_k=%.6lf a0=%.6lf\n\n",para.X_obs_k/1000.0,para.Y_obs_k/1000.0,para.Z_obs_k/1000.0,para.a0*1e+06);

        //fprintf(FileWrite,"%d %d %d %d %d %.2lf",time.year,time.month,time.day,time.hour,time.minute,time.second);
		fprintf(FileWrite,"%s %.6lf %.6lf %.6lf %.6lf\n",para.name,para.X_obs_k/1000.0,para.Y_obs_k/1000.0,para.Z_obs_k/1000.0,para.a0*1e+06);			
		}
		 
	}
	fclose(file);
	fclose(FileWrite);
}

int main()
{
	FILE *file;
	FILE *FileWrite;
	
	FigureOut(file,FileWrite);

	return 0;
}

 2.2国外RINEX4.00标准

国外标准分为CNAV1与CNAV2。笔者写了CNAV1的代码供参考。

国外标准的公式与国内IGMAS标准有一个区别，要将A0 = Aref + DeltaA换为A0 = sqrtA的平方。

CNAV1参数位置如下，按下图提取我们所需参数。

 C语言代码如下。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#define MAX_DATA_SIZE 300000

const double GM = 3.986005e+14, we = 7.292115e-05,
			 pi = 3.1415926, miu = 3.986004418e+14,
             Aref_MEO = 27906100, Aref_IGSOGEO = 42162200;

typedef struct
{
	int year;  //年
	int month; //月
	int day;   //日
	int hour;  //时
	int minute;//分
	double second;//秒
}Time;

typedef struct
{
    char name[10],name1[20];
    double a0, a1, a2;//钟差 钟速 钟漂
    double A_DOT, Crs, Delta_n0, M0;//星历参数版本号  轨道半径的正弦调和改正项的振幅 参考时刻卫星平均角速度与计算值之差 Toe时刻的平近点角
    double Cuc, e, Cus, Delta_A;//纬度幅角的余弦调和改正项的振幅 轨道偏心率 纬度幅角的正弦调和改正项的振幅 轨道长半径的平方根
    double toe, Cic, OMEGA_0, Cis;//星历表参考历元 轨道倾角的余弦调和项 升交点赤经 轨道倾角的正弦调和项
    double i0, Crc, w, OMEGA_DOT;//轨道倾角 轨道半径的余弦调和项 近地点角距 升交点赤经变化率
    double IDOT, n0_DOT, SatType, t_op;//轨道倾角变化率 电文来源 BDS的整周计数 长半轴变化率
    double SlSAl_oe, SISAI_ocb, SISAI_oc1, SISAI_oc2;//卫星精度 卫星自主健康标识 TGD1星上设备时延差 TGD2星上设备时延差
    double ISC_B1Cd, Spare_x1, TGD_B1Cp, TGD_B2ap;//信息的发射时间（ BDT 的周内秒）钟差参数版本号/时钟数据龄期
    //参考时刻卫星平均角速度与值之差的变化率  留空 卫星轨道类型
    double SISMAI , Health, B1C, IODC;
    double t_tm, Spare_x2, IODE, Note;

    double secondOfweek1, secondOfweek2, secondOfweek;
    double n, n0, t_observe, t, tk , Delta_t, toc, Mk, Vk, phi_k, A0, Ak, sqrt_A, nA;
    double _6u, _6r, _6i, uk, rk, ik, Xk, Yk, Ek, E1, Delta_nA;
    double OMEGA_k, X_obs_k, Y_obs_k, Z_obs_k, X__obs_k, Y__obs_k, Z__obs_k;
}Parameter;
double UTC2BDST(int year, int month, int day, int hour, int minute, double second, double *secondOfweek) {
    double JD, WN;
    if (month <= 2) {
        year -= 1;
        month += 12;
    }
    JD = 365.25 * year + int(30.6001 * (month + 1)) + day + 1720981.5 + hour / 24.0 + minute / 1440 +
         second / 86400;  // 计算儒略日
    WN = int((JD - 2453736.5) / 7); // WN:GPS_week number 目标时刻的GPS周
    *secondOfweek = (int(JD - 2453736.5) - (7.0 * WN)) * 24 * 3600.0;  // tGPS:目标时刻的GPS秒 减去14秒为BDT

    return *secondOfweek;
}

double UTC2BDS(int year, int month, int day, int hour, int minute, double second, double *secondOfweek)
{ 
   int DayofYear = 0;
   int DayofMonth = 0;
   for (int i = 2006; i < year; i++)  
   {
      if ((i % 4 == 0 && i % 100 != 0) || i % 400 == 0)
                DayofYear += 366;
      else
               DayofYear += 365;
   }
   for (int i = 1; i < month; i++)   
   {
         if (i == 1 || i == 3 || i == 5 || i == 7 || i == 8 || i == 10 || i ==12)
                DayofMonth += 31;
          else if (i == 4 || i == 6 || i == 9 || i == 11)
                DayofMonth += 30;
            else
            {
            if ((year % 4 == 0 && year % 100 != 0) || year % 400 == 0)
                DayofMonth += 29;
            else
                DayofMonth += 28;
          }
    }
   int Day;
   Day = DayofMonth + day + DayofYear-1;//2006年1月1日起始历元
   *secondOfweek = Day % 7 * 86400 + hour * 3600 + minute * 60 + second;

   return *secondOfweek;
}

void FigureOut(FILE *file,FILE *FileWrite)
{
	int data_count = 0;
    char ReadPath[50] = "E:\\Buffer\\data18_RINEX4.00.txt";
    char WritePath[50] = "E:\\Buffer\\result18_RINEX4.00.txt";
	file = fopen(ReadPath,"r");
	FileWrite = fopen(WritePath,"w");
	if(!file||!FileWrite)
	{
		printf("open the file error!\n");
	}
	
	Time time;//年月日时分秒
    Parameter para;
	
		while(1) 
		{
			if(feof(file))
			{
				printf("文件读取结束");
				break;
			}

        fgets(para.name1,15,file);
        //printf("%s\n",para.name);

        if( strcmp(para.name1,"> EPH C10 CNV1") == 0 || strcmp(para.name1,"> EPH C16 CNV1") == 0 || strcmp(para.name1,"> EPH C26 CNV1") == 0 || strcmp(para.name1,"> EPH C29 CNV1") == 0 ||
        strcmp(para.name1,"> EPH C30 CNV1") == 0 || strcmp(para.name1,"> EPH C35 CNV1") == 0 || strcmp(para.name1,"> EPH C36 CNV1") == 0 || strcmp(para.name1,"> EPH C38 CNV1") == 0 ||
        strcmp(para.name1,"> EPH C39 CNV1") == 0 || strcmp(para.name1,"> EPH C40 D1") == 0 )
        {
            fscanf(file, "%s %d %d %d %d %d %lf %lf %lf %lf",&para.name,&time.year,&time.month,
                   &time.day,&time.hour,&time.minute,&time.second,&para.a0, &para.a1, &para.a2);
            //printf(" %d %d %d %d %d %.20lf %.20lf %.20lf %.20lf \n", time.year, time.month,
            //time.day, time.hour, time.minute, time.second, para.a0, para.a1, para.a2);
            if(time.year == 2023 && time.month == 6 && time.day == 1 && time.hour == 0 && time.minute == 0 && time.second == 0.0)
            {
                fscanf(file, " %lf %lf %lf %lf ", &para.A_DOT,&para.Crs,&para.Delta_n0,&para.M0);
                //printf(" %.10lf \t %.10lf \t %.10lf \t %.10lf \n", para.A_DOT, para.Crs, para.Delta_n0, para.M0);

                fscanf(file, " %lf %lf %lf %lf ", &para.Cuc, &para.e, &para.Cus, &para.sqrt_A);
                //printf(" %.10lf \t %.10lf \t %.10lf \t %.10lf \n", para.Cuc, para.e, para.Cus, para.sqrt_A);

                fscanf(file, " %lf %lf %lf %lf ", &para.toe, &para.Cic, &para.OMEGA_0, &para.Cis);
                //printf(" %.10lf \t %.10lf \t %.10lf \t %.10lf \n", para.toe, para.Cic, para.OMEGA_0, para.Cis);

                fscanf(file, " %lf %lf %lf %lf ", &para.i0, &para.Crc, &para.w, &para.OMEGA_DOT);
                //printf(" %.10lf \t %.10lf \t %.10lf \t %.10lf \n", para.i0, para.Crc, para.w, para.OMEGA_DOT);

                fscanf(file, " %lf %lf %lf %lf ", &para.IDOT, &para.n0_DOT, &para.SatType, &para.t_op);
                //printf(" %.20lf \t %.20lf \t %.20lf \t %.20lf \n", para.IDOT, para.n0_DOT, para.SatType, para.t_op);

                fscanf(file, " %lf %lf %lf %lf ", &para.SlSAl_oe, &para.SISAI_ocb, &para.SISAI_oc1, &para.SISAI_oc2);
                //printf(" %.10lf \t %.10lf \t %.10lf \t %.10lf \n", para.SlSAl_oe, para.SISAI_ocb, para.SISAI_oc1, para.SISAI_oc2);

                fscanf(file, " %lf %lf %lf %lf ", &para.ISC_B1Cd, &para.TGD_B1Cp, &para.TGD_B2ap);
                //printf(" %.10lf \t %.10lf \n", para.ISC_B1Cd, para.TGD_B1Cp, para.TGD_B2ap);

                fscanf(file, " %lf %lf %lf %lf ", &para.SISMAI , &para.Health, &para.B1C, &para.IODC);
                //printf(" %.10lf \t %.10lf \n", para.SISMAI, para.Health, para.B1C, para.TGD_B2ap);

                fscanf(file, " %lf %lf %lf %lf ", &para.t_tm,  &para.Note);
                //printf(" %.10lf \t %.10lf \n", para.t_tm, para.Note);

                /*计算归化时间t_k*/
                //t_observe = UTC2BDST(time.year,time.month,time.day,time.hour,time.minute,time.second,&secondOfweek1);
                //toc = UTC2BDST(time.year,time.month,time.day,time.hour,time.minute,time.second,&secondOfweek2);
                //Delta_t = a0 + a1*(t_observe-toc) + a2*(t_observe-toc)*(t_observe-toc);
                //t = t_observe - Delta_t;
                para.t = UTC2BDS(time.year,time.month,time.day,time.hour,time.minute ,time.second,&para.secondOfweek);
                para.tk = para.t - para.toe - 14;
                if(para.tk>302400)
                    para.tk -= 604800;
                else
                if(para.tk<-302400)
                    para.tk += 604800;
                printf("归化时间tk = %.20lf\n",para.tk);

                /*计算参考时刻的长半轴*/
                if(strcmp(para.name,"C11")==0||strcmp(para.name,"C12")==0||strcmp(para.name,"C14")==0||strcmp(para.name,"C19")==0||strcmp(para.name,"C20")==0||strcmp(para.name,"C21")==0
                   ||strcmp(para.name,"C22")==0||strcmp(para.name,"C23")==0||strcmp(para.name,"C24")==0||strcmp(para.name,"C25")==0||strcmp(para.name,"C26")==0||strcmp(para.name,"C27")==0
                   ||strcmp(para.name,"C28")==0||strcmp(para.name,"C29")==0||strcmp(para.name,"C30")==0||strcmp(para.name,"C32")==0||strcmp(para.name,"C33")==0||strcmp(para.name,"C34")==0
                   ||strcmp(para.name,"C35")==0||strcmp(para.name,"C36")==0||strcmp(para.name,"C37")==0||strcmp(para.name,"C41")==0||strcmp(para.name,"C42")==0||strcmp(para.name,"C43")==0
                   ||strcmp(para.name,"C44")==0||strcmp(para.name,"C45")==0||strcmp(para.name,"C46")==0||strcmp(para.name,"C57")==0||strcmp(para.name,"C58")==0)
                {para.A0 = Aref_MEO +  para.Delta_A;}
                else
                {para.A0 = Aref_IGSOGEO +  para.Delta_A;}
                printf("参考时刻的长半轴A0 = %.20lf\n",para.A0);

                /*计算长半轴Ak*/
                para.Ak = para.A0 +  para.A_DOT*para.tk;
                printf("长半轴Ak = %.20lf\n",para.Ak);

                /*计算参考时刻的卫星平均角速度n0*/
                para.n0 = sqrt(miu/pow(para.A0,3));
                printf("参考时刻的卫星平均角速度n0 = %.20lf\n",para.n0);

                /*计算卫星平均角速度的偏差Delta_nA*/
                para.Delta_nA = para.Delta_n0 + (1.0/2)*para.n0_DOT*para.tk;
                printf("卫星平均角速度的偏差Delta_nA = %.20lf\n",para.Delta_nA);

                /*计算改正后的卫星平均角速度nA*/
                para.nA = para.n0 + para.Delta_nA;
                printf("改正后的卫星平均角速度nA = %.20lf\n",para.nA);

                /*计算平近点角Mk*/
                para.Mk = para.M0 + para.nA*para.tk;
                printf("平近点角Mk = %.20lf\n",para.Mk);

                /*计算偏近点角Ek*/
                para.Ek = para.Mk;
                para.Ek = para.Mk + para.e*sin(para.Ek);
                for(int i=0;i<2;i++)
                {
                    para.Ek = para.Mk + para.e*sin(para.Ek);
                }

                //int count = 0;
                //Ek=Mk;
                //E1=Mk+e*sin(Ek);
                //while ( abs(E1-Ek)>1e-10 )
                //{
                //	Ek=E1;
                //	E1=Mk+e*sin(E1);
                //	count++;
                //	if(count>1e+06)
                //	{
                //		break;
                //	}
                //}
                //Ek=E1;
                printf("偏近点角Ek = %.20lf\n",para.Ek);

                /*计算真近点角Vk*/
                //Vk = atan(sqrt(1-e*e)*sin(Ek)/(cos(Ek)-e));
                //Vk=2*atan2(sqrt((1+e)/(1-e))*tan(Ek/2),1);
                para.Vk = atan2(sqrt(1-para.e*para.e)*sin(para.Ek),(cos(para.Ek)-para.e));
                printf("真近点角Vk = %.20lf\n",para.Vk);

                /*计算纬度幅角fai_k*/
                para.phi_k = para.Vk + para.w;
                printf("升交距角fai_k = %.20lf\n",para.phi_k);

                /*计算摄动改正量*/
                para._6u = para.Cus*sin(2*para.phi_k) + para.Cuc*cos(2*para.phi_k);
                para._6r = para.Crs*sin(2*para.phi_k) + para.Crc*cos(2*para.phi_k);
                para._6i = para.Cis*sin(2*para.phi_k) + para.Cic*cos(2*para.phi_k);
                printf("纬度幅角改正量 _6u = %.20lf\n",para._6u);
                printf("径向距离改正量 _6r = %.20lf\n",para._6r);
                printf("轨道倾角改正量 _6i = %.20lf\n",para._6i);

                /*计算摄动改正后的 纬度幅角uk 径向距离rk 轨道倾角ik*/
                para.uk = para.phi_k + para._6u;
                para.rk = para.Ak*(1-para.e*cos(para.Ek)) + para._6r;
                para.ik = para.i0 + para.IDOT*para.tk + para._6i;
                printf("升交距角uk = %.20lf\n",para.uk);
                printf("卫星矢径rk = %.20lf\n",para.rk);
                printf("轨道倾角ik = %.20lf\n",para.ik);

                /*计算卫星在轨道平面坐标系的坐标*/
                para.Xk = para.rk*cos(para.uk);
                para.Yk = para.rk*sin(para.uk);
                printf("Xk = %.20lf\n",para.Xk);
                printf("Yk = %.20lf\n",para.Yk);

                if(strcmp(para.name,"C01")==0||strcmp(para.name,"C02")==0||strcmp(para.name,"C03")==0||strcmp(para.name,"C04")==0||strcmp(para.name,"C05")==0)
                {
                    /*GEO*/
                    para.OMEGA_k = para.OMEGA_0 + para.OMEGA_DOT*para.tk - we*para.toe;

                    para.X__obs_k = para.Xk*cos(para.OMEGA_k)-para.Yk*cos(para.ik)*sin(para.OMEGA_k);
                    para.Y__obs_k = para.Xk*sin(para.OMEGA_k)+para.Yk*cos(para.ik)*cos(para.OMEGA_k);
                    para.Z__obs_k = para.Yk*sin(para.ik);

                    para.X_obs_k = cos(we*para.tk)*para.X__obs_k + sin(we*para.tk)*cos(-5.0/180*pi)*para.Y__obs_k + sin(we*para.tk)*sin(-5.0/180*pi)*para.Z__obs_k;
                    para.Y_obs_k = -1*sin(we*para.tk)*para.X__obs_k + cos(we*para.tk)*cos(-5.0/180*pi)*para.Y__obs_k + cos(we*para.tk)*sin(-5.0/180*pi)*para.Z__obs_k;
                    para.Z_obs_k = -1*sin(-5.0/180*pi)*para.Y__obs_k + cos(-5.0/180*pi)*para.Z__obs_k;

                }
                else
                {
                    /*MEO/IGSO*/
                    /*计算改正后的升交点经度OMEGA_k*/
                    para.OMEGA_k = para.OMEGA_0 + (para.OMEGA_DOT-we)*para.tk - we*para.toe;
                    printf("升交点经度OMEGA_k = %.20lf\n",para.OMEGA_k);

                    /*计算卫星在CGCS2000中的坐标*/
                    para.X_obs_k = para.Xk*cos(para.OMEGA_k)-para.Yk*cos(para.ik)*sin(para.OMEGA_k);
                    para.Y_obs_k = para.Xk*sin(para.OMEGA_k)+para.Yk*cos(para.ik)*cos(para.OMEGA_k);
                    para.Z_obs_k = para.Yk*sin(para.ik);
                }
                printf("%s\n",para.name);
                printf("X_obs_k=%.6lf Y_obs_k=%.6lf Z_obs_k=%.6lf a0=%.6lf\n\n",para.X_obs_k/1000.0,para.Y_obs_k/1000.0,para.Z_obs_k/1000.0,para.a0*1e+06);

                //fprintf(FileWrite,"%.6lf %.6lf %.6lf %.6lf\n",para.X_obs_k/1000.0,para.Y_obs_k/1000.0,para.Z_obs_k/1000.0,para.a0);
                fprintf(FileWrite,"%.12lf  %.12lf  %.12lf  %.12lf  %.12lf  %.12lf  %.12lf  %.12lf  %.12lf\n",para.X_obs_k/1000.0,para.Y_obs_k/1000.0,para.Z_obs_k/1000.0,para.a0,para.a1,para.a2,para.sqrt_A,para.Ek,para.e);
            }
        }

	}
	fclose(file);
	fclose(FileWrite);
}

int main()
{
	FILE *file;
	FILE *FileWrite;
	
	FigureOut(file,FileWrite);

	return 0;
}

如有错误，欢迎批评指正。

利用广播星历解算北斗卫星位置及精度分析
weixin_48630797的博客
09-28 1万+
网上已经有了比较多关于如何利用播发的广播星历来解算卫星位置的blog，此文章的目的是为了记录一下学习成果的同时回馈一下（毕竟看过很多blog了，但是从来没写过）。
基于广播星历的北斗定位解算原理（基于C语言和MATLAB实现）_卫星位置解算
04-20 1331
/ FILE *pfileWrite = fopen(“文件路径（相对路径或绝对路径均可）”,“w”);// FILE *pfile = fopen(“文件路径（相对路径或绝对路径均可）”,“r”);printf(“\t摄动改正量δu = %.20lf\n”,delta_u);printf(“\t摄动改正量δr = %.20lf\n”,delta_r);printf(“\t输入参考时间t:(小于60分钟)\t====>”);printf(“\t升交距角φk = %.20lf\n”,Phik);
卫星轨道计算
03-11
利用开普勒六参数进行卫星轨道计算的详细说明
16参数卫星轨道位置计算（基于Matlab）
m0_73984031的博客
02-16 2542
1.存储方式不同：单元格数组一般存储混合数据类型的表格或其他非均匀数据集，而数值数组以规则的网格形式存储数据，所有元素都按照相同的规则排列，转换有利于下标来索引。2.运算上区别：单元格数组不能直接进行数学运算，数值数组支持各种数学运算，例如加法、乘法、矩阵运算等。使用vpa函数来设置有效数字的位数。
【一文读懂】卫星轨道的轨道参数（六根数）和位置速度矢量转换及其在终端距离、相对速度和多普勒频移计算中的应用
最新发布
u011376987的博客
02-09 3283
本文首先介绍了卫星轨道的六根数参数：半长轴 a、偏心率 e、轨道倾角 i、升交点赤经 Ω、近地点幅角 ω 以及真近点角 ν，并探讨了如何利用这些参数计算卫星与地面终端之间的距离、相对速度及由此产生的多普勒频移。利用开普勒运动公式计算当前瞬间的轨道半径；根据真近点角 ν 在围心坐标系中计算位置矢量和速度矢量；利用旋转矩阵（依次旋转 ω、i、Ω）将围心坐标系的状态矢量转换到地心惯性坐标系；最终结合地面终端信息，计算卫星与终端的欧氏距离、沿视线的相对速度及多普勒频移。
利用北斗BDS广播星历，进行卫星位置计算
qq_45568498的博客
12-19 8891
利用RINEX文件中的星历和卫星发射时刻计算北斗MEO/IGSO卫星和GEO卫星的位置
（二）GNSS定位中的卫星轨道位置计算
qq_40185784的博客
03-21 8263
计算卫星位置
用广播星历计算卫星位置
04-09
用C++语言实现利用广播星历计算卫星位置
读取广播星历计算卫星位置
01-20
广播星历是由全球导航卫星系统（如GPS、GLONASS、Galileo或北斗）向地球广播的一系列参数，用于计算地面接收机的精确位置。本项目以VS2010作为开发环境，编写了一个程序，能够根据用户输入的卫星编号和时间，计算出...
C# 计算GPS卫星位置（使用广播星历）
06-25
首先，广播星历是GPS卫星向地面发送的一种预计算数据，包含了卫星的轨道参数，如偏心率、升交点径、平均运动和历元等。这些参数是通过数学模型预测的，可帮助用户估算卫星在特定时间的精确位置。要解析广播星历数据...
基于Matlab的北斗卫星位置计算及精度分析.pdf
06-27
通过Matlab编程实现北斗卫星位置计算的关键步骤通常包括利用广播星历提供的开普勒轨道参数和轨道摄动改正项参数，来计算卫星在特定时刻的瞬时坐标。接着，将计算得到的卫星位置与同时刻通过精密星历获得的卫星位置...
卫星轨道及定位计算流程
07-07
详细介绍根据导航星历的卫星位置计算，及定位算法流程，非常实用
（武汉大学GPS作业）根据广播星历文件计算卫星位置
10-19
武汉大学GPS原理与应用课程结课作业，根据广播星历文件计算卫星位置，给有需要的小伙伴们参考^_^！ （项目工程环境Visual Studio 2010，在release文件夹中可找到独立运行的exe程序）
卫星轨道计算.pdf
11-17
本资源主要介绍卫星轨道，主要包括：卫星运动特性、卫星的空间定位、卫星覆盖计算、轨道摄动、轨道对通信系统性能的影响、卫星发射
北斗卫星广播星历计算卫星位置
weixin_47168560的博客
10-12 3856
通过星历计算北斗igso，meo，geo卫星位置方法
利用广播星历计算北斗二号卫星位置
qq_43381998的博客
12-05 2262
北斗二号广播星历由6个开普勒轨道根数、3个摄动线性改正项和6个周期改正项系数 等一共 15个描述轨道特征的参数和1个星历参考时间组成，其广播星历参数的更新时间是 1h。从每周六/每周日子夜零点起算，每周的这时候都会归零。升交点赤经、近地点角距、轨道倾角，三者被称为定向参数。这三者唯一确定了卫星在三维空间中的方向。有3个参数可以描述卫星的位置：真近点角ν、偏近点角E和平近点角M。三者被称为位置参数。
如何利用北斗卫星导航系统实现高精度的精密定轨和广播星历轨道精度评估？
11-01
为了掌握北斗卫星导航系统的精密定轨和广播星历轨道精度评估方法，建议参考这篇深入解析北斗卫星导航系统精密定轨与广播星历精度的资料：《北斗卫星导航系统精密定轨与广播星历精度深度解析》。这本资料详细探讨了北斗系统在轨道精度方面的技术特点和分析手段，特别是基于IGS跟踪网数据的应用，对于提高定位准确性具有重要的指导意义。 参考资源链接：[北斗卫星导航系统精密定轨与广播星历精度深度解析](https://wenku.csdn.net/doc/6n7rn1ztpy?spm=1055.2569.3001.10343) 具体来说，精密定轨涉及到对北斗GEO卫星的PANDA软件分析，通过载波相位观测数据进行精确度分析，从而发现卫星切向分量的恒定偏差问题。解决这一问题的关键在于综合利用时间比对和平滑伪距观测数据，以此来获得更高的轨道精度。在操作过程中，你需要注意数据采集的精确度以及数据处理的技术细节，以确保结果的可靠性。 另外，为了评估广播星历的轨道精度，可以通过对IGS、M-GEX、Furgo、武汉大学BETS和ILRS等跟踪网提供的观测数据进行分析。这些数据能够帮助研究者得到切向和法向的轨道精度，从而对北斗系统的导航定位能力进行科学的评估。 在了解了精密定轨和广播星历轨道精度评估的方法后，你可以更好地掌握北斗卫星导航系统的性能特点。对于追求更高精度定位服务的用户来说，这些知识尤为重要。如果你希望对北斗卫星导航系统的性能进行深入研究，本资料将是你不可多得的学习资源，它将带你深入了解卫星导航系统在精密定轨和轨道精度评估方面的核心技术和实践应用。 参考资源链接：[北斗卫星导航系统精密定轨与广播星历精度深度解析](https://wenku.csdn.net/doc/6n7rn1ztpy?spm=1055.2569.3001.10343)
评论 14
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

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

余额充值