GPS卫星位置的计算（C++程序计算）

GPS卫星位置的计算（C++程序计算）
GPS卫星位置的计算
利用C++编写了一段能计算单一瞬时卫星坐标的程序，在运行程序之前，需做部分准备工作：（1）在F盘下建立一名为“单一卫星广播星历”的txt文件。

（2）从“广播星历.txt”文件中拷贝从卫星PRN号开始的8行数据到“单一卫星广播星历.txt”中
（3）在编辑选项中，将全部的“D”替换为“E”。

下面为我所选取的一个广播星历：
18 06 8 25 6 0 0.0-2.472363412380E-04-1.023*********E-12
0.000000000000E+00
1.410000000000E+02-1.721875000000E+01
4.502687555010E-09 1.413760604187E+00
-7.990747690201E-07 7.598234573379E-03 1.118145883083E-05 5.153709835052E+03
4.536000000000E+05-1.303851604462E-08-
1.095067942661E-01 1.527369022369E-07
9.571235745530E-01 1.640000000000E+02-2.656176299285E+00-8.0374********E-09
-5.193073455211E-10 1.000000000000E+00
1.389000000000E+03 0.000000000000E+00
2.000000000000E+00 0.000000000000E+00-1.024*********E-08 1.410000000000E+02
4.464490000000E+05 4.000000000000E+00
程序设计部分：
#include
#include
int main()
{
int i = 0;
double n[50], n0, nn, t, tk, Mk, Ek, Vk, Yk, Gu, Gr, Gi, uk, rk, ik, xk, yk, zk, X, Y, Z, Lk, UT, yy, mm, JD, gpsz;
FILE *fp;
fp = fopen("F:\\单一卫星广播星历.txt", "r");
if (fp == NULL)
{
printf ("文件打开失败!\n");
return 0;
}
while (! feof (fp))
{
fscanf(fp, "%lf", &n[i]);
i++;
}
n0 = (sqrt(3986005E+8))/pow(n[17], 3);
nn = n0 + n[12];/*计算卫星运行的平均角速度*/
UT = n[4] + (n[5] / 60) + (n[7] / 3600);/*民用日的时分秒化为实数时*/
if (n[1] >= 80)/*广播星历中年只有后两位，化为4位，参考1980年1月6日0点*/ {
if (n[1] == 80 && n[2] == 1 && n[3] < 6)
{
n[1] = n[1] + 2000;
}
n[1] = n[1] + 1900;
}
else
{
n[1] = n[1] + 2000;
}
if (n[2] <= 2)
{
yy = n[1] - 1;
mm = n[2] + 12;
}
if (n[2] > 2)
{
yy = n[1];
mm = n[2];
}
JD = (int)(365.25 * yy) + (int)(30.6001 * (mm + 1)) + n[3] + (UT / 24) + 1720981.5;/*化为儒略日*/
gpsz = (int)((JD - 2444244.5) / 7);/*计算GPS周*/
t = (JD - 2444244.5 - 7 * gpsz) * 24 * 3600;/*得出GPS秒*/
tk = t - n[18];/*tk1为中间值，用以判断tk与正负302400的关系，然后返回到tk上*/ while (tk > 302400 || tk < -302400) {
if (tk > 302400)
{
tk = tk - 604800;
}
else
{
tk = tk + 604800;
}
}/*计算归化观测时间*/
Mk = n[13] + nn * tk;/*观测时刻的卫星平近点角*/
Ek = Mk;
Ek = Mk + n[15] * sin(Ek);
Ek = Mk + n[15] * sin(Ek);/*迭代两次计算观测时刻的偏近点角*/ Vk = atan(sqrt(1 - n[15] * n[15]) * sin(Ek)) / (cos(Ek) - n[15]);/*真近点角*/
Yk = Vk + n[24];/*升交距角*/
Gu = n[14] * cos(2 * Yk) + n[16] * sin(2 * Yk);
Gr = n[23] * cos(2 * Yk) + n[11] * sin(2 * Yk);
Gi = n[19] * cos(2 * Yk) + n[21] * sin(2 * Yk);/*摄动改正项*/
uk = Yk + Gu;
rk = n[17] * n[17] * (1 - n[15] * cos(Ek)) + Gr;
ik = n[22] + Gi + n[26] * tk;/*经摄动改正后的升交距角、卫星矢径、轨道倾角*/
xk = rk * cos(uk);
yk = rk * sin(uk);
zk = 0;/*卫星在轨道坐标系的坐标*/
Lk = n[20] + (n[25] - 7.29211515E-5) * tk - 7.29211515E-5 * n[18];/*观测时刻t的升交点经度*/
X = xk * cos(Lk) - yk * cos(ik) * sin(Lk);
Y = xk * sin(Lk) + yk * cos(ik) * cos(Lk);
Z = yk * sin(ik);/*卫星在WGS-84坐标系的坐标*/
printf("该卫星在WGS-84坐标系中的坐标为：\nX = %lf m\nY = %lf m\nZ = %lf m\n", X, Y, Z);
fclose(fp);
return 0;
}
计算结果：
该卫星在WGS-84坐标系中的坐标为：
X = 9223153.692525 m
Y = 24133486.931401 m
Z = 6032585.919385 m。