GPS卫星位置的计算(C++程序计算)

1、GPS卫星位置的计算利用C+编写了一段能计算单一瞬时卫星坐标的程序，在运行程序之前，需做部分准备工作：（1） 在F盘下建立一名为“单一卫星广播星历”的txt文件。（2） 从“广播星历.txt”文件中拷贝从卫星PRN号开始的8行数据到“单一卫星广播星历.txt”中（3） 在编辑选项中，将全部的“D”替换为“E”。下面为我所选取的一个广播星历：18 06 8 25 6 0 0.0-2.0E-04-1.5E-12 0.0E+00 1.0E+02-1.0E+01 4.0E-09 1.7E+00 -7.1E-07 7.9E-03 1.3E-05 5.2E+03 4.0E+05-1.2E-08-1.1E-

2、01 1.9E-07 9.0E-01 1.0E+02-2.5E+00-8.9E-09 -5.1E-10 1.0E+00 1.0E+03 0.0E+00 2.0E+00 0.0E+00-1.7E-08 1.0E+024.0E+05 4.0E+00程序设计部分：#include#includeint main()int i = 0;double n50, 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:

3、单一卫星广播星历.txt, r);if (fp = NULL)printf (文件打开失败!n);return 0;while (! feof (fp)fscanf(fp, %lf, &ni);i+;n0 = (sqrt(E+8)/pow(n17, 3);nn = n0 + n12;/*计算卫星运行的平均角速度*/UT = n4 + (n5 / 60) + (n7 / 3600);/*民用日的时分秒化为实数时*/if (n1 = 80)/*广播星历中年只有后两位，化为4位，参考1980年1月6日0点*/if (n1 = 80 & n2 = 1 & n3 6)n1 = n1 + 2000; n1

4、 = n1 + 1900;elsen1 = n1 + 2000;if (n2 2)yy = n1;mm = n2;JD = (int)(365.25 * yy) + (int)(30.6001 * (mm + 1) + n3 + (UT / 24) + .5;/*化为儒略日*/gpsz = (int)(JD - .5) / 7);/*计算GPS周*/t = (JD - .5 - 7 * gpsz) * 24 * 3600;/*得出GPS秒*/tk = t - n18;/*tk1为中间值，用以判断tk与正负的关系，然后返回到tk上*/while (tk | tk )tk = tk - ;else

5、tk = tk + ;/*计算归化观测时间*/Mk = n13 + nn * tk;/*观测时刻的卫星平近点角*/Ek = Mk;Ek = Mk + n15 * sin(Ek);Ek = Mk + n15 * sin(Ek);/*迭代两次计算观测时刻的偏近点角*/Vk = atan(sqrt(1 - n15 * n15) * sin(Ek) / (cos(Ek) - n15);/*真近点角*/Yk = Vk + n24;/*升交距角*/Gu = n14 * cos(2 * Yk) + n16 * sin(2 * Yk);Gr = n23 * cos(2 * Yk) + n11 * sin(2

6、* Yk);Gi = n19 * cos(2 * Yk) + n21 * sin(2 * Yk);/*摄动改正项*/uk = Yk + Gu;rk = n17 * n17 * (1 - n15 * cos(Ek) + Gr;ik = n22 + Gi + n26 * tk;/*经摄动改正后的升交距角、卫星矢径、轨道倾角*/xk = rk * cos(uk);yk = rk * sin(uk);zk = 0;/*卫星在轨道坐标系的坐标*/Lk = n20 + (n25 - 7.E-5) * tk - 7.E-5 * n18;/*观测时刻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 = %