基于广播星历的北斗定位解算原理（基于C语言和MATLAB实现） |
您所在的位置：网站首页 › 卫星轨道的机动原理图 › 基于广播星历的北斗定位解算原理（基于C语言和MATLAB实现） |
基于广播星历的北斗定位解算原理（基于C语言和MATLAB实现）
|
一、利用广播星历解算卫星位置
网上已经有很多解算北斗卫星位置的资料和算法，此文章的目的只是记录自己的学习过程，并且积累一些写文章博客的经验。本文先用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，可以参考北斗卫星导航系统办公室发表的文件。如下图所示
16/18参数在广播星历中各自所对应的位置在网上均可找到，此处不做赘述。 二、实现步骤（具体每一步讲解会在代码中体现并带有注释）
值得注意的是，toc转到BDT是从2006年1月1日（周日）作为起始点；因为BDT和GPST相差14秒，所以计算归化时间tk时，需要人为减去14秒，调整到和GPST同步。计算tk时，t相当于是对卫星位置做预测，且北斗卫星每隔一个小时更新一次星历，故t的取值小于60分钟。 二、18参数算法一、前言 18参数模型相对比16参数改动不是特别大，一个是公式的变动，另一个是广播星历对应参数的位置部分有变动。 公式的变动主要体现在下图，具体的还是参考中国卫星导航系统管理办公室提供的官方文件。
二、实现步骤（具体每一步讲解会在代码中体现并带有注释）
三、C语言源码 本代码仅适用于计算北斗卫星国内格式的文件，国外格式的文件需要做相对应的修正。 本代码可以直接读取整个大文件，并计算出所有北斗卫星的位置。 //**************************************** Message ***********************************// //技术交流：[email protected] //关注CSDN博主：“彼稷” //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 #include #include #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 ***********************************// %//技术交流：[email protected] %//关注CSDN博主：“彼稷” %//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卫星分布：
IGSO卫星分布：
GEO卫星分布：
所有北斗卫星分布：
三、最后 本文只是记录学习过程中的一些感想并复盘每阶段的任务，代码仍有优化的空间、内容有谬误之处请大家多多指正。
|
【本文地址】
公司简介
联系我们
今日新闻 |
点击排行 |
|
推荐新闻 |
图片新闻 |
|
专题文章 |