武汉大学卫星导航原理课程作业——广播星历与卫星位置盘算
武汉大学卫星导航原理课程作业采用C++语言编译，实现rinex3.04格式的广播星历文件读取和卫星的位置解算。
[*]作业内容
[*]根据所给定的广播星历文件，盘算当天多颗GPS卫星（3颗及以上）在24h内每隔30s在瞬时地心地固系下的坐标。
[*]将上述盘算得到的GPS卫星位置与精密星历对比（15min间隔），评估广播星历盘算出的卫星位置的误差，绘制出卫星的轨道误差时序图。
[*]给定文件说明
[*]广播星历文件：BRDC00IGS_R_20213480000_01D_MN.rnx
[*]精密星历文件：WUM0MGXFIN_20213480000_01D_15M_ORB.SP3
[*]rinex 04格式说明文档。给出的广播星历文件格式为rinex 04，参数含义参照说明文档附录A5、A6。
一、广播星历文件的读取
1.首先读取所有内容，存储到lines容器里面，并记录下头文件竣事的位置，便于后续获取数据
vector<string> lines;
string end_of_header = " END OF HEADER";
int end_of_line = 0;//文件头结束的位置
void readfile(string filename)
{
int i = 0;
ifstream infile(filename);
if (!infile.is_open())
{
cerr << "open error!" << endl;
}
string line="0";
while (getline(infile, line))
{
if (line == end_of_header)
{
end_of_line = i;
}
lines.push_back(line);
i = i + 1;
}
infile.close();
}//读取广播星历文件 2.界说结构体并进行数据分配
class Data
{
private:
struct data_body
{
string sPRN;//卫星号
int TOC_Y;
int TOC_M;
int TOC_D;
int TOC_H;
int TOC_Min;
int TOC_Sec;//卫星钟的参考时刻，年月日...
double sa0;//钟差
double sa1;//钟速
double sa2;//钟漂
//广播轨道-1
double IODE;//星历发布时间
double Crs;
double deltan;
double M0;
//广播轨道-2
double Cuc;
double e;//轨道偏心角
double Cus;
double sqrtA;//长半轴平方根
//广播轨道-3
double TOE;//星历的参考时刻（GPS周内秒）
double Cic;
double OMEGA;//升交点赤经
double Cis;
//广播轨道-4
double i0;//参考时刻的轨道倾角
double Crc;
double omega;//近地点角距
double deltaomega;//
//广播轨道-5
double IDOT;
double L2code;
double GPSweek;//GPS周数
double L2flag;
//广播轨道-6
double aACC;//卫星精度
double sHEA;//卫星健康状态
double TGD;//TGD
double IODC;//钟的数据龄期
//广播轨道-7
double TTN;//电文发送时刻
double fit;//拟合区间
double none1;
double none2;//备用
};
struct position
{
int GPSweek;
double GPSsec;
double x;
double y;
double z;//地心地固坐标系下的坐标
};
public:
vector <data_body> data;
vector <position> data_P;
void getdata()
{
for (int j = end_of_line + 1; j <= size(lines) - 1; j = j + 8)
{
data_body data1;
istringstream iss(lines);
iss >> data1.sPRN;
iss >> data1.TOC_Y >> data1.TOC_M >> data1.TOC_D >> data1.TOC_H >> data1.TOC_Min >> data1.TOC_Sec;
iss>> data1.sa0 >> data1.sa1 >> data1.sa2;
istringstream iss1(lines);
iss1 >> data1.IODE >> data1.Crs >> data1.deltan >> data1.M0;
istringstream iss2(lines);
iss2 >> data1.Cuc >> data1.e >> data1.Cus >> data1.sqrtA;
istringstream iss3(lines);
iss3 >> data1.TOE >> data1.Cic >> data1.OMEGA >> data1.Cis;
istringstream iss4(lines);
iss4 >> data1.i0 >> data1.Crc >> data1.omega >> data1.deltaomega;
istringstream iss5(lines);
iss5 >> data1.IDOT >> data1.L2code >> data1.GPSweek >> data1.L2flag;
istringstream iss6(lines);
iss6 >> data1.aACC >> data1.sHEA >> data1.TGD >> data1.IODC;
istringstream iss7(lines);//逐行解析数据
vector<double> row;//存储数据的临时容器
double value=0;
while (iss7 >> value)
{
row.push_back(value);
}
while (row.size() < 4)
{
row.push_back(0);
}//数据不足4个时，补齐为0
data1.TTN = row;
data1.fit = row;
data1.none1 = row;
data1.none2 = row;
data.push_back(data1);//将解析完的数据添加到容器里面
}
} 二、卫星位置解算
1.盘算步骤：
https://i-blog.csdnimg.cn/direct/b7ec376c83b14ce5ac9a8204b1855928.png
https://i-blog.csdnimg.cn/direct/083bb72a7b374c15bf353697951b5621.png
2.时间转换和最佳历元选取
需要将世界时转换成GPS时，同时需要查找时间差最近的历元，以此历元的数据进行卫星位置解算，代码如下：
double get_GPST(int Y, int M, int D, int H, int min, double S)
{
double y = 0;
double m = 0;
if (M <= 2)
{
y = Y - 1;
m = M + 12;
}
else
{
y = Y;
m = M;
}
double UT = 0;//世界时
double MJD = 0;//简化儒略日
UT = H + min / 60.0 + S / 3600.0;
MJD = int(365.25*y) + int(30.6001*(m + 1)) + D + UT / 24 - 679019;
int GPSweek = 0;
double secofweek = 0;
GPSweek = int((MJD - 44244) / 7);//GPS周
secofweek = (MJD - 44244 - GPSweek * 7) * 86400;//GPS周秒
return secofweek;
}//计算GPS周内秒的函数
int select_epoch(double secofweek,string PRN)
{
int n = 0;
double min_s = 10000;
double Min;
for (int i = 0; i <= size(data) - 1; i++)
{
if (data.sPRN == PRN)
{
Min = fabs(secofweek - data.TOE);
if (Min <= min_s)
{
n = i;
min_s = Min;
}
}
}
return n;
}//选择时间差最小的历元计算
void caculate(int k,double secofweek)
{
double A = pow(data.sqrtA,2);
double n0 = sqrt(u / pow(A, 3));
double n = n0 + data.deltan;//计算卫星的平均角速度n
double tk = secofweek - data.TOE;//计算规划时间
if (tk > 302400)
{
tk = tk - 604800;
}
if (tk < -302400)
{
tk = tk + 604800;
}//考虑跨周的影响
double Mk = data.M0 + n * tk;//计算平近点角
double E0 = 0;
double E1 = Mk;
while (abs(E1-E0)>1e-6)
{
E0 = E1;
E1 = Mk + data.e*sin(E0);
}//计算偏近点角
double cosv = (cos(E1) - data.e) / (1 - data.e * cos(E1));
double sinv = sqrt(1 - pow(data.e, 2))*sin(E1) / (1 - data.e * cos(E1));
double Vk = atan2(sinv, cosv);
double w = Vk + data.omega;
double Ou = data.Cuc*cos(2 * w) + data.Cus*sin(2 * w);
double Or = data.Crc*cos(2 * w) + data.Crs*sin(2 * w);
double Oi = data.Cic*cos(2 * w) + data.Cis*sin(2 * w);
double uk = w + Ou;
double rk = A * (1 - data.e*cos(E1)) + Or;
double ik = data.i0 + Oi + data.IDOT*tk;
double xk1 = rk * cos(uk);
double yk1 = rk * sin(uk);
double Ok = data.OMEGA + (data.deltaomega - Oe)*tk - Oe * data.TOE;
double xk = xk1 * cos(Ok) - yk1 * cos(ik)*sin(Ok);
double yk = xk1 * sin(Ok) + yk1 * cos(ik)*cos(Ok);
double zk = yk1 * sin(ik);
position P;
P.GPSweek = data.GPSweek;
P.GPSsec = secofweek;
P.x = xk;
P.y = yk;
P.z = zk;
data_P.push_back(P);//把坐标数据添加到容器里面
} 三、结果文件保存
将盘算所得结果按照指定要求保存在文本文件中：
void savefile()
{
ofstream outfile("广播星历.txt");
outfile <<fixed<< setprecision(3);
for (int i = 0; i <= size(data_P) - 1; i++)
{
outfile << data_P.GPSweek << " " << data_P.GPSsec << " " << data_P.x << " " << data_P.y << " " << data_P.z << endl;
}
}//将广播星历的数据写入文件 武汉大学测绘学院在读本科生一枚，本次的内容就分享到这里了，渴望各人多多关注，后面会连续给各人分享更多的优质内容。本次作业的所有测试用例和源代码都可以免费下载，渴望各人点点关注支持一下。
【免费】武汉大学卫星导航原理课程作业-广播星历与卫星位置盘算.zip资源-CSDN文库
免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。
页:
[1]