根据导航文件计算卫星坐标的程序

#include "iostream.h" #include "stdio.h" #include "afx.h" #include "math.h" void readdata();//读取导航文件数据 void dealdata();//处理导航文件数据 typedef struct Navedata //定义导航数据存储结构 { int PRN;//卫星的PRn号 int year,month,day,hour,minute; //卫星钟参考时刻 double second; double a0,a1,a2; //卫星钟飘参数 double IODE,Crs,DeltN,M0; //数据星历发布时间，在轨道径向方向上周期改正正弦的振幅 double Cuc,e,Cus,sqrtA; //轨道延迹方向上周期改正余弦振幅 、扁心率、轨道延迹方向上周期改正正弦振幅 、长半轴平方根 double Toe,Cic,OMEGA,Cis; //星历参考时刻、轨道倾角周期改正余弦项振幅、参考时刻升交点赤径主项、轨道倾角周期改正正弦项振幅 double i0,Crc,omega,OMEGADOT; //参考时间轨道倾角、在轨道径向方向上周期改正余余弦的振幅、近地点角距、升交点赤径在赤道平面中的长期变化 double IDOT,L2C,GPSWeek,L2P; ////轨道倾角变化率、gps周等 double AccuracyofSat,HealthofSat,TGD,IODC; //卫星精度、卫星健康、电离层群迟改正数 double eight1,eight2,eight3,eight4;//导航数据第八行的参数（计算中不用） }; typedef struct xyz { int prn; double x,y,z; }; #define allstarnum 9 //导航文件中需要计算的卫星的个数 struct Navedata stars[allstarnum];//定义每个卫星数据的结构数组 struct xyz positions[allstarnum]; //每个卫星坐标的结构数组 void readdata()//读取数据文件 { CStdioFile fpt; CString string; BOOL isend=fpt.Open("50191100.10n",CFile::modeRead);//打开导航文件 if(!isend) printf("file open failed!"); int headlinenum=0; while (isend)//读取头文件 { isend=fpt.ReadString(string); headlinenum++; int headend=string.Find("END OF HEADER"); if (-1!=headend) break; } fpt.SeekToBegin();//把指针指到文件的开头 for (int i=0;i<headlinenum;i++)//把头文件部分去除 isend=fpt.ReadString(string); for (i=0;i<allstarnum;i++)//读取数据部分 { //读取数据部分的第一行prn,year等 isend=fpt.ReadString(string); string.Replace('D','e'); sscanf(string,"%d %d %d %d %d %d %lf %lf %lf %lf", &stars[i].PRN,&stars[i].year,&stars[i].month,&stars[i].day,&stars[i].hour,&stars[i].minute,&stars[i].second, &stars[i].a0, &stars[i].a1, &stars[i].a2); //读取IDOT,Crs,deltn，m0； isend=fpt.ReadString(string); string.Replace('D','e'); sscanf(string,"%lf %lf %lf %lf",&stars[i].IDOT,&stars[i].Crs,&stars[i].DeltN,&stars[i].M0); //读取cuc，e，cus，sqrta isend=fpt.ReadString(string); string.Replace('D','e'); sscanf(string,"%lf %lf %lf %lf",&stars[i].Cuc, &stars[i].e, &stars[i].Cus, &stars[i].sqrtA); //读取toe，cic，omega，cis isend=fpt.ReadString(string); string.Replace('D','e'); sscanf(string,"%lf %lf %lf %lf",&stars[i].Toe, &stars[i].Cic, &stars[i].OMEGA, &stars[i].Cis); //读取i0，crc，omega，omegadot isend=fpt.ReadString(string); string.Replace('D','e'); sscanf(string,"%lf %lf %lf %lf",&stars[i].i0, &stars[i].Crc, &stars[i].omega, &stars[i].OMEGADOT); //读取idot，l2c，gpsweek，l2p isend=fpt.ReadString(string); string.Replace('D','e'); sscanf(string,"%lf %lf %lf %lf", &stars[i].IDOT, &stars[i].L2C, &stars[i].GPSWeek, &stars[i].L2P); //读取accracyofsat等 isend=fpt.ReadString(string); string.Replace('D','e'); sscanf(string,"%lf %lf %lf %lf",&stars[i].AccuracyofSat, &stars[i].HealthofSat, &stars[i].TGD, &stars[i].IODC); //读取eight等 isend=fpt.ReadString(string); string.Replace('D','e'); sscanf(string,"%lf %lf %lf %lf",&stars[i].eight1,&stars[i].eight2,&stars[i].eight3,&stars[i].eight4); } fpt.Close();//关闭文件 } void dealdata()//gps卫星位置的计算函数 { for (int i=0;i<allstarnum;i++) { double udot,f,n0,n,M,t,E,deltu,deltr,delti,u,r,s; double xdot,ydot,zOMEGA,GAST,L; n0=3.9860047E+14/stars[i].sqrtA;//计算参考时刻的平均角速度 n=n0+stars[i].DeltN;//观测时刻卫星的平均角速度 t=5*86400 + 2*3600+30*60;//观测时刻一周内的秒数 M=stars[i].M0+n*(t-stars[i].Toe);//计算观测瞬间卫星的平近点角 double E0=M; while (1)// 用迭代法计算偏近点角 { E=M+stars[i].e*sin(E0); if (fabs(E-E0)<1e-12) break; else E0=E; } f=atan(sqrt(1-pow(stars[i].e,2))*sin(E)/(cos(E)-stars[i].e));//计算真近点角 udot=stars[i].omega+f;//计算升交距角 deltu=stars[i].Cuc*cos(2*udot)+stars[i].Cus*sin(2*udot);//计算摄动改正项 deltr=stars[i].Crc*cos(2*udot)+stars[i].Crs*sin(2*udot); delti=stars[i].Cic*cos(2*udot)+stars[i].Cis*sin(2*udot); u=udot+deltu;//对上面的进行摄动改正 r=pow(stars[i].sqrtA,2)*(1-stars[i].e*cos(E))+deltr; s=stars[i].i0+delti+stars[i].IDOT*(t-stars[i].Toe); xdot=r*cos(u);//计算卫星在轨道面坐标系中的位置 ydot=r*sin(u); zOMEGA=stars[i].OMEGA+stars[i].OMEGADOT*(t-stars[i].Toe);//计算观测瞬间升交点的经度 double we=7.292115e-5; GAST=stars[i].GPSWeek+we*t; L=zOMEGA-GAST; positions[i].x=xdot*cos(L)-ydot*cos(s)*sin(L);//计算卫星在瞬时地球坐标系中的位置 positions[i].y=xdot*sin(L)+ydot*cos(s)*cos(L); positions[i].z=ydot*sin(s); positions[i].prn=stars[i].PRN; } } void main() { readdata(); dealdata(); for (int i=0;i<allstarnum;i++)//输出每个卫星的prn和位置坐标 { printf("卫星PRN为 %d 在瞬时地球坐标系中的位置：x=%lf,y=%lf,z=%lf \n",positions[i].prn,positions[i].x,positions[i].y,positions[i].z); } }