GPS卫星位置的计算(C 程序计算)
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
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);/*民用日的时分秒化为实数时*/
程序设计部分：
#include<stdio.h> #include<math.h>
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) {
return 0; }
计算结果：
该卫星在 WGS-84 坐标系中的坐标为： X = 9223153.692525 m Y = 24133486.931401 m Z = 6032585.919385 m
4.536000000000E+05-1.303851604462E-08-1.095067942661E-01 1.527369022369E-07 9.571235745530E-01 1.640000000000E+02-2.656176299285E+00-8.037477650349E-09 -5.193073455211E-10 1.000000000000E+00 1.389000000000E+03 0.000000000000E+00 2.000000000000E+00 0.000000000000E+00-1.024454832077E-08 1.410000000000E+02 4.464490000000E+05 4.000000000000E+00
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.023181539495E-12 0.000000000000E+00
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];/*升交距角*/
1.410000000000E+0来自百度文库-1.721875000000E+01 4.502687555010E-09 1.413760604187E+00 -7.990747690201E-07 7.598234573379E-03 1.118145883083E-05 5.153709835052E+03
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; } 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);/*民用日的时分秒化为实数时*/
程序设计部分：
#include<stdio.h> #include<math.h>
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) {
return 0; }
计算结果：
该卫星在 WGS-84 坐标系中的坐标为： X = 9223153.692525 m Y = 24133486.931401 m Z = 6032585.919385 m
4.536000000000E+05-1.303851604462E-08-1.095067942661E-01 1.527369022369E-07 9.571235745530E-01 1.640000000000E+02-2.656176299285E+00-8.037477650349E-09 -5.193073455211E-10 1.000000000000E+00 1.389000000000E+03 0.000000000000E+00 2.000000000000E+00 0.000000000000E+00-1.024454832077E-08 1.410000000000E+02 4.464490000000E+05 4.000000000000E+00
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.023181539495E-12 0.000000000000E+00
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];/*升交距角*/
1.410000000000E+0来自百度文库-1.721875000000E+01 4.502687555010E-09 1.413760604187E+00 -7.990747690201E-07 7.598234573379E-03 1.118145883083E-05 5.153709835052E+03
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);