（武汉大学GPS作业）根据广播星历文件计算卫星位置

#include "stdafx.h" #include "afx.h" #include "math.h" #define GM 3.986005e14 //基本广播星历块 struct EPHEMERISBLOCK //每小时一个卫星对应一个基本星历块 { //PRN号 int PRN; double a0,a1,a2;//时间改正数 //六个轨道参数 double IODE,Crs,Deltan,M0;// ORBIT - 1 double Cuc,e,Cus,SqrtA;// ORBIT - 2 double Toe,Cic,OMEGA,Cis;// ORBIT - 3 double i0,Crc,omega,OMEGAdot;// ORBIT - 4 double IDOT,GpsWeekNumber,L2C,L2P;// ORBIT - 5 double SatAccuracy,SatHealth,TGD,IODC;// ORBIT - 6 }; struct GPSTIME { int weekno; double weekSecond; }; int Calendar2GpsTime (int nYear, int nMounth, int nDay, int nHour, int nMinute, double dSecond,double &WeekSecond); int ReadBrodcastEphemeris(CString strEpheNam, int &EphemerisBlockNum); EPHEMERISBLOCK *m_pGpsEphemeris = NULL; //读广播星历文件，数据存储与上面定义的指针中 //参数：strEpheNam表示广播星历文件的完整路径 // EphemerisBlockNum 返回读取到的星历块个数 int ReadBrodcastEphemeris(CString strEpheNam, int &EphemerisBlockNum) { int HeadLineNum = 0; int WeekNo; double WeekSecond; //打开文件 CStdioFile pfEph; BOOL IsEn = pfEph.Open(strEpheNam, CFile::modeRead); if(!IsEn) return 0; //读入头文件 CString strLine; while(IsEn) { IsEn = pfEph.ReadString(strLine); HeadLineNum++; int index = strLine.Find("END OF HEADER"); if( -1 != index ) break; } //计算星历块数 int AllNum = 0; while(IsEn) { IsEn = pfEph.ReadString(strLine); AllNum++; } //临时读入星历块 EphemerisBlockNum = (AllNum + 1) / 8; m_pGpsEphemeris = new EPHEMERISBLOCK[EphemerisBlockNum]; GPSTIME *pGpsTime = new GPSTIME[EphemerisBlockNum]; if(!m_pGpsEphemeris || !pGpsTime) return 0; //将文件指针调整到数据位置 pfEph.SeekToBegin(); for(int i=0; i<HeadLineNum; i++) IsEn = pfEph.ReadString(strLine); //定义读取的参数 int mPrn;//卫星号PRNo int year,month,day,hour,minute;//卫星钟参考时刻 double msecond; 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;//卫星精度、卫星健康、电离层群迟改正数 for(int i=0; i<EphemerisBlockNum; i++) { //读取卫星PRN号，星历参考时间 IsEn = pfEph.ReadString(strLine); strLine.Replace('D', 'e'); sscanf(strLine,"%d %d %d %d %d %d %lf %lf %lf %lf",&mPrn,&year,&month,&day,&hour,&minute,&msecond, &a0, &a1, &a2); year += 2000; WeekNo = Calendar2GpsTime(year, month, day, hour, minute, msecond, WeekSecond); pGpsTime[i].weekno = WeekNo; pGpsTime[i].weekSecond = WeekSecond; IsEn = pfEph.ReadString(strLine); strLine.Replace('D', 'e'); sscanf(strLine,"%lf %lf %lf %lf",&IODE, &Crs, &DeltN, &M0); //读 Cuc,e,Cus,sqrtA IsEn = pfEph.ReadString(strLine); strLine.Replace('D', 'e'); sscanf(strLine,"%lf %lf %lf %lf",&Cuc, &e, &Cus, &sqrtA); //Toe,Cic,OMEGA,Cis; IsEn = pfEph.ReadString(strLine); strLine.Replace('D', 'e'); sscanf(strLine,"%lf %lf %lf %lf",&Toe, &Cic, &OMEGA, &Cis); //i0,Crc,w,OMEGADOT IsEn = pfEph.ReadString(strLine); strLine.Replace('D', 'e'); sscanf(strLine,"%lf %lf %lf %lf",&i0, &Crc, &omega, &OMEGADOT); //IDOT,L2Cod,GPSWeek,L2PCod IsEn = pfEph.ReadString(strLine); strLine.Replace('D', 'e'); sscanf(strLine,"%lf %lf %lf %lf", &IDOT, &L2C, &GPSWeek, &L2P); //AccuracyofSat,HealthofSat,TGD,IODC IsEn = pfEph.ReadString(strLine); strLine.Replace('D', 'e'); sscanf(strLine,"%lf %lf %lf %lf",&AccuracyofSat, &HealthofSat, &TGD, &IODC); // IsEn = pfEph.ReadString(strLine); //赋值 m_pGpsEphemeris[i].PRN = mPrn; m_pGpsEphemeris[i].a0 = a0; m_pGpsEphemeris[i].a1 = a1; m_pGpsEphemeris[i].a2 = a2; //&IODE, &Crs, &DeltN, &M0 m_pGpsEphemeris[i].IODE = IODE; m_pGpsEphemeris[i].Crs = Crs; m_pGpsEphemeris[i].Deltan = DeltN; m_pGpsEphemeris[i].M0 = M0; //&Cuc, &e, &Cus, &sqrtA m_pGpsEphemeris[i].Cuc = Cuc; m_pGpsEphemeris[i].e = e; m_pGpsEphemeris[i].Cus = Cus; m_pGpsEphemeris[i].SqrtA = sqrtA; //Toe,Cic,OMEGA,Cis; m_pGpsEphemeris[i].Toe = Toe; m_pGpsEphemeris[i].Cic = Cic; m_pGpsEphemeris[i].OMEGA = OMEGA; m_pGpsEphemeris[i].Cis = Cis; //i0,Crc,omega,OMEGADOT m_pGpsEphemeris[i].i0 = i0; m_pGpsEphemeris[i].Crc = Crc; m_pGpsEphemeris[i].omega = omega; m_pGpsEphemeris[i].OMEGAdot = OMEGADOT; //iDOT,L2Cod,GPSWeek,L2PCod m_pGpsEphemeris[i].IDOT = IDOT; m_pGpsEphemeris[i].L2C = L2C; m_pGpsEphemeris[i].L2P = L2P; m_pGpsEphemeris[i].GpsWeekNumber = GPSWeek; //AccuracyofSat,HealthofSat,TGD,IODC m_pGpsEphemeris[i].SatAccuracy = AccuracyofSat; m_pGpsEphemeris[i].SatHealth = HealthofSat; m_pGpsEphemeris[i].TGD = TGD; m_pGpsEphemeris[i].IODC = IODC; } pfEph.Close(); if(pGpsTime) {delete []pGpsTime; pGpsTime = NULL; } return 1; } //从年 月 日 转换为GPS 周秒 int Calendar2GpsTime (int nYear, int nMounth, int nDay, int nHour, int nMinute, double dSecond,double &WeekSecond) { int DayofMonth = 0; int DayofYear = 0; int weekno = 0; int dayofweek; int m; if (nYear < 1980 || nMounth < 1 || nMounth > 12 || nDay < 1 || nDay > 31) return -1; //计算从1980年到当前的前一年的天数 for( m = 1980 ; m < nYear ; m++ ) { if ( (m%4 == 0 && m%100 != 0) || (m%400 == 0) ) { DayofYear += 366; } else DayofYear += 365; } //计算当前一年内从元月到当前前一月的天数 for( m = 1;m < nMounth; m++) { if(m==1 || m==3 || m==5 || m==7 || m==8 || m==10 || m==12) DayofMonth += 31; else if (m==4 || m==6 || m==9 || m==11) DayofMonth += 30; else if (m ==2) { if ( (nYear%4 == 0 && nYear%100 != 0) || (nYear%400 == 0) ) DayofMonth += 29; else DayofMonth += 28; } } DayofMonth = DayofMonth + nDay - 6;//加上当月天数/减去1980年元月的6日 weekno = (DayofYear + DayofMonth) / 7;//计算GPS周 dayofweek = (DayofYear + DayofMonth) % 7; //计算GPS 周秒时间 WeekSecond = dayofweek*86400 + nHour*3600 + nMinute*60 + dSecond; return weekno; } int _tmain(int argc, _TCHAR* argv[]) { CString strEpheNam; strEpheNam="P0201880.18N"; //文件名 int * EphemerisBlockNum = new int[100]; ReadBrodcastEphemeris(strEpheNam,*EphemerisBlockNum); int nYear,nMonth,nDay,nHour,nMinute,number; double dSecond,WeekSecond; printf("\n***************************** xxx班 xxx 广播星历文件计算卫星位置小程序*****************************\n\n");//在此处写上你自己的班级姓名 printf("请输入要计算的时刻的年份："); scanf("%d",&nYear); printf("\n月："); scanf("%d",&nMonth); printf("\n日："); scanf("%d",&nDay); printf("\n时："); scanf("%d",&nHour); printf("\n分："); scanf("%d",&nMinute); printf("\n秒："); scanf("%lf",&dSecond); printf("\n输入星历块编号（1-11的整数)："); scanf("%d",&number); Calendar2GpsTime(nYear,nMonth,nDay, nHour, nMinute, dSecond,WeekSecond); double t; t = WeekSecond; double n0; n0=sqrt(GM)/pow(m_pGpsEphemeris[number-1].SqrtA,3); double n; n=n0+m_pGpsEphemeris[number-1].Deltan; //计算卫星在观测时刻的平均角速度n double M; M=m_pGpsEphemeris[number-1].M0 + n*(t - m_pGpsEphemeris[number-1].Toe); //计算平近角点M double E; E=M; for(int i=0;i<10;i++) { E=M+m_pGpsEphemeris[number-1].e*sin(E); //迭代的方法计算偏近点角E } double f; f=atan((sqrt(1 - m_pGpsEphemeris[number-1].e*m_pGpsEphemeris[number-1].e)*sin(E)) / (cos(E) - m_pGpsEphemeris[number-1].e)); //计算真近点角f double u1; u1=m_pGpsEphemeris[number-1].omega+f; //计