C语言计算GPS卫星位置

“GPS广播星历计算卫星位置和速度”及“GPS伪距定位”计算试验1．试验内容及上交成果1.1 试验内容应用C语言按预定格式（函数、输入输出变量之名称、类型）编写“GPS广播星历计算卫星位置和速度”函数SatPos_Vel( )、“GPS伪距定位”函数Positioning( )。

将此两个函数组成文件F2.cpp，并包含于文件GPS_Positioning.cpp中。

编译、连接并运行文件GPS_Positioning.cpp，逐一时刻读取广播星历（Ephemeris.dat）、观测时间及伪距、卫星号（Observation.dat）信息，计算WGS84坐标系中观测时刻相应的卫星位置、速度以及载体位置，结果保存于文件Position.dat中。

1.2 上交成果磁盘文件F2.cpp、Position.dat，并存于“学号作者中文姓名”目录中。

2．函数说明2.1 星历文件读取函数void EFileReading(Efile)功能：读取星历文件，给星历数据结构体Ephemeris赋值。

输入变量：EFile 字符串，文件名。

2.2 观测数据读取函数int ObsReading(fp_Obs,Time,Rho,Mark)功能：从文件Observation.dat中读取某一时刻的伪距、卫星号。

读取成功函数值返回“1”，失败返回“－1”（读错，或至文件尾）。

输入变量：fp_Obs 文件指针；输出变量：Time double，时间（秒）；Rho double[12]，伪距（米）；Mark int[12]，卫星号，“－1”表示此通道无卫星、无伪距。

2.3 最小二乘估计函数int LeastSquareEstimation(Y,A,P,m,n, X)功能：最小二乘方法求解观测方程Y=AX+ε，其中观测值方差阵的逆阵为P（也称为权阵），得未知参数X。

成功返回“1”，失败返回“－1”（亏秩）。

输入变量：Y double[m]，观测方程自由项（米）；A double[m×n]，系数阵（无量纲），按第1行第1、2……n元素，第2行……顺序存放；P double[m]，权矩阵对角线元素，0或1（无量纲）；m int，观测值个数；n int，未知参数个数；输出变量：X double[n]，未知参数（米）。

坐标转换源代码--GPS定位程序(C++)GPS数据处理中为了满足不同的需要，处理的数据要进行坐标转换，得到在不同坐标系统下的结果，下面是笛卡尔坐标系，大地坐标系，站心地平坐标系（线型和极坐标形式）之间的转换源代码：头文件：#ifndef _COORDCOVERT_H#define _COORDCOVERT_H#include "stdlib.h"//WGS-84椭球体参数const double a=6378137.0;//长半轴const double flattening=1/298.257223563;//扁率const double delta=0.0000001;typedef struct tagCRDCARTESIAN{double x;double y;double z;}CRDCARTESIAN;typedef CRDCARTESIAN *PCRDCARTESIAN;//笛卡尔坐标系typedef struct tagCRDGEODETIC{double longitude;double latitude;double height;}CRDGEODETIC;typedef CRDGEODETIC *PCRDGEODETIC;//大地坐标系typedef struct tagCRDTOPOCENTRIC{double northing;double easting;double upping;}CRDTOPOCENTRIC;typedef CRDTOPOCENTRIC *PCRDTOPOCENTRIC; //站心地平坐标系（线坐标形式）typedef struct tagCRDTOPOCENTRICPOLAR{ double range;double azimuth;double elevation;}CRDTOPOCENTRICPOLAR;typedef CRDTOPOCENTRICPOLAR *PCRDTOPOCENTRICPOLAR;//站心地平坐标系（极坐标形式）//由笛卡尔坐标转换为大地坐标void CartesianToGeodetic (PCRDGEODETIC pcg, PCRDCARTESIAN pcc, double dSemiMajorAxis, double dFlattening);//pcg：指向所转换出的大地坐标的指针；//pcc：指向待转换的笛卡尔坐标的指针；//dSemiMajorAxis：参考椭球的长半轴；//dFlattening：参考椭球的扁率。

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+001.410000000000E+02-1.721875000000E+01 4.502687555010E-09 1.413760604187E+00-7.990747690201E-07 7.598234573379E-03 1.118145883083E-05 5.153709835052E+034.536000000000E+05-1.303851604462E-08-1.095067942661E-01 1.527369022369E-079.571235745530E-01 1.640000000000E+02-2.656176299285E+00-8.0374********E-09-5.193073455211E-10 1.000000000000E+00 1.389000000000E+03 0.000000000000E+002.000000000000E+00 0.000000000000E+00-1.024*********E-08 1.410000000000E+024.464490000000E+05 4.000000000000E+00程序设计部分：#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){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 mY = 24133486.931401 mZ = 6032585.919385 m。

C语言计算GPS卫星位置GPS（全球定位系统）是一种利用卫星定位来确定全球地理位置的技术。

GPS系统由一组卫星和地面接收器组成，可以帮助我们在地球上的任何位置确定自己的准确位置。

在C语言中，我们可以利用一些数学和物理公式来计算GPS卫星位置。

首先，我们需要了解GPS系统的原理。

GPS卫星运行在地球轨道上，同时向地球发送信号。

地面接收器接收到卫星发出的信号，并通过计算信号的距离和时间来确定自身的位置。

计算GPS卫星位置的关键是测量信号的传播时间。

当接收器接收到来自至少4个不同卫星的信号时，我们可以使用以下公式来计算GPS卫星的位置：速度=距离/时间由于信号的传播速度恒定（大约为光速），我们可以知道信号的传播时间等于距离与速度的比。

在C语言中，我们可以使用以下代码计算GPS卫星的位置：```c#include <stdio.h>#include <math.h>typedef structdouble x;double y;double z;} Point3D;Point3D calculateSatellitePosition(double distance, double latitude, double longitude, double altitude)Point3D position;position.x = (altitude + distance) * cos(E) * cos(longitude);position.y = (altitude + distance) * cos(E) * sin(longitude);return position;int maidouble latitude = 37.7749; // 纬度（假设）double longitude = -122.4194; // 经度（假设）double altitude = 0.0; // 海拔（假设）Point3D satellitePosition =calculateSatellitePosition(distance, latitude, longitude, altitude);printf("GPS卫星位置：(x=%.2f, y=%.2f, z=%.2f)\n", satellitePosition.x, satellitePosition.y, satellitePosition.z);return 0;```在这段代码中，我们首先定义了光速的常量。

c语言 gps课程设计一、教学目标本课程的目标是让学生掌握C语言在GPS领域的应用。

通过本课程的学习，学生将能够：1.理解GPS系统的基本原理和工作流程。

2.掌握C语言的基本语法和编程技巧。

3.能够使用C语言进行GPS数据的采集、解析和处理。

4.能够利用C语言实现简单的GPS导航功能。

二、教学内容本课程的教学内容主要包括以下几个部分：1.GPS系统的基本原理和工作流程。

2.C语言的基本语法和编程技巧。

3.GPS数据的采集、解析和处理。

4.GPS导航功能的实现。

三、教学方法为了达到本课程的教学目标，我们将采用多种教学方法，包括：1.讲授法：用于讲解GPS系统的基本原理和工作流程，以及C语言的基本语法和编程技巧。

2.案例分析法：通过分析具体的GPS应用案例，让学生了解GPS技术的实际应用。

3.实验法：通过实验让学生亲手操作，掌握GPS数据的采集、解析和处理方法，以及实现简单的导航功能。

四、教学资源为了支持本课程的教学内容和教学方法的实施，我们将准备以下教学资源：1.教材：选用合适的C语言和GPS相关教材，为学生提供理论学习的参考。

2.参考书：提供相关的参考书籍，拓展学生的知识面。

3.多媒体资料：制作课件和视频资料，直观地展示GPS系统和C语言的相关内容。

4.实验设备：准备GPS接收器、计算机等实验设备，让学生进行实际操作。

五、教学评估为了全面、客观地评估学生在C语言GPS课程中的学习成果，我们将采取以下评估方式：1.平时表现：通过学生在课堂上的参与度、提问回答、小组讨论等表现来评估其学习态度和理解程度。

2.作业：布置与课程内容相关的编程作业，评估学生对C语言编程和GPS应用的掌握情况。

3.考试：定期进行理论知识考试和编程实践考试，全面评估学生的知识掌握和实际应用能力。

4.项目报告：要求学生完成一个GPS相关的编程项目，通过项目报告评估学生的综合应用能力。

六、教学安排本课程的教学安排如下：1.教学进度：按照教材和大纲进行，确保覆盖所有重要知识点。

GPS编码格式及C语⾔解码有关磁偏⾓和地图定位的问题：地图的⽅向：上北、下南、左西、右东是⼤多数地图的⽅向，但这可不是通⽤原则，如果地图上有⽅向标，可以通过⽅向标了解到这些。

地磁极是接近南极和北极的，但并不和南极、北极重合，⼀个约在北纬72°、 西经96°处；⼀个约在南纬70°、东经150°处。

磁北极距地理北极⼤约相差1500km。

现在磁北极的位置在加拿⼤北⽅，在⼀天中磁北极的位置也是不停的变动，它的轨迹⼤致为⼀椭圆形，磁北极平均每天向北以40m。

⽬前磁北极正在逐渐离开加拿⼤，⼤约于2005年进⼊俄罗斯境内。

东经25度地区，磁偏⾓在1-2度之间；北纬25度以上地区，磁偏⾓⼤于2度；若在西经低纬度地区，磁偏⾓是5-20度；西经45度以上，磁偏⾓为25-50度。

在我国，正常情况下，磁偏⾓最⼤可达6度，⼀般情况为2-3度。

在我国除部分磁⼒异常的地⽅外，⼀般磁偏⾓都是西偏。

磁偏⾓还是不断有规律变化的，地图上的磁偏⾓只是测图时的磁偏⾓（磁北⽐真北偏左，加上磁偏⾓；磁北⽐真北偏右，减去磁偏⾓；在我国⼀般是加上）。

使⽤地图本⾝所注的磁偏⾓要注意出版年限，地图太⽼误差较⼤。

GPS 0183协议GGA、GLL、GSA、GSV、RMC、VTG解释$GPGGA例：$GPGGA,092204.999,4250.5589,S,14718.5084,E,1,04,24.4,19.7,M,,,,0000*1F字段0：$GPGGA，语句ID，表明该语句为Global Positioning System Fix Data（GGA）GPS定位信息字段1：UTC 时间，hhmmss.sss，时分秒格式字段2：纬度ddmm.mmmm，度分格式（前导位数不⾜则补0）字段3：纬度N（北纬）或S（南纬）字段4：经度dddmm.mmmm，度分格式（前导位数不⾜则补0）字段5：经度E（东经）或W（西经）字段6：GPS状态，0=未定位，1=⾮差分定位，2=差分定位，3=⽆效PPS，6=正在估算字段7：正在使⽤的卫星数量（00 - 12）（前导位数不⾜则补0）字段8：HDOP⽔平精度因⼦（0.5 - 99.9）字段9：海拔⾼度（-9999.9 - 99999.9）字段10：地球椭球⾯相对⼤地⽔准⾯的⾼度字段11：差分时间（从最近⼀次接收到差分信号开始的秒数，如果不是差分定位将为空）字段12：差分站ID号0000 - 1023（前导位数不⾜则补0，如果不是差分定位将为空）字段13：校验值$GPGLL例：$GPGLL,4250.5589,S,14718.5084,E,092204.999,A*2D字段0：$GPGLL，语句ID，表明该语句为Geographic Position（GLL）地理定位信息字段1：纬度ddmm.mmmm，度分格式（前导位数不⾜则补0）字段2：纬度N（北纬）或S（南纬）字段3：经度dddmm.mmmm，度分格式（前导位数不⾜则补0）字段4：经度E（东经）或W（西经）字段5：UTC时间，hhmmss.sss格式字段6：状态，A=定位，V=未定位字段7：校验值$GPGSA例：$GPGSA,A,3,01,20,19,13,,,,,,,,,40.4,24.4,32.2*0A字段0：$GPGSA，语句ID，表明该语句为GPS DOP and Active Satellites（GSA）当前卫星信息字段1：定位模式，A=⾃动⼿动2D/3D，M=⼿动2D/3D字段2：定位类型，1=未定位，2=2D定位，3=3D定位字段3：PRN码（伪随机噪声码），第1信道正在使⽤的卫星PRN码编号（00）（前导位数不⾜则补0）字段4：PRN码（伪随机噪声码），第2信道正在使⽤的卫星PRN码编号（00）（前导位数不⾜则补0）字段5：PRN码（伪随机噪声码），第3信道正在使⽤的卫星PRN码编号（00）（前导位数不⾜则补0）字段6：PRN码（伪随机噪声码），第4信道正在使⽤的卫星PRN码编号（00）（前导位数不⾜则补0）字段7：PRN码（伪随机噪声码），第5信道正在使⽤的卫星PRN码编号（00）（前导位数不⾜则补0）字段8：PRN码（伪随机噪声码），第6信道正在使⽤的卫星PRN码编号（00）（前导位数不⾜则补0）字段9：PRN码（伪随机噪声码），第7信道正在使⽤的卫星PRN码编号（00）（前导位数不⾜则补0）字段10：PRN码（伪随机噪声码），第8信道正在使⽤的卫星PRN码编号（00）（前导位数不⾜则补0）字段11：PRN码（伪随机噪声码），第9信道正在使⽤的卫星PRN码编号（00）（前导位数不⾜则补0）字段12：PRN码（伪随机噪声码），第10信道正在使⽤的卫星PRN码编号（00）（前导位数不⾜则补0）字段13：PRN码（伪随机噪声码），第11信道正在使⽤的卫星PRN码编号（00）（前导位数不⾜则补0）字段14：PRN码（伪随机噪声码），第12信道正在使⽤的卫星PRN码编号（00）（前导位数不⾜则补0）字段15：PDOP综合位置精度因⼦（0.5 - 99.9） 字段16：HDOP⽔平精度因⼦（0.5 - 99.9）字段17：VDOP垂直精度因⼦（0.5 - 99.9）字段18：校验值$GPGSV例：$GPGSV,3,1,10,20,78,331,45,01,59,235,47,22,41,069,,13,32,252,45*70字段0：$GPGSV，语句ID，表明该语句为GPS Satellites in View（GSV）可见卫星信息字段1：本次GSV语句的总数⽬（1 - 3）字段2：本条GSV语句是本次GSV语句的第⼏条（1 - 3）字段3：当前可见卫星总数（00 - 12）（前导位数不⾜则补0）字段4：PRN 码（伪随机噪声码）（01 - 32）（前导位数不⾜则补0）字段5：卫星仰⾓（00 - 90）度（前导位数不⾜则补0）字段6：卫星⽅位⾓（00 - 359）度（前导位数不⾜则补0）字段7：信噪⽐（00－99）dbHz字段8：PRN 码（伪随机噪声码）（01 - 32）（前导位数不⾜则补0）字段9：卫星仰⾓（00 - 90）度（前导位数不⾜则补0）字段10：卫星⽅位⾓（00 - 359）度（前导位数不⾜则补0）字段11：信噪⽐（00－99）dbHz字段12：PRN 码（伪随机噪声码）（01 - 32）（前导位数不⾜则补0）字段13：卫星仰⾓（00 - 90）度（前导位数不⾜则补0）字段14：卫星⽅位⾓（00 - 359）度（前导位数不⾜则补0）字段15：信噪⽐（00－99）dbHz字段16：校验值$GPRMC例：$GPRMC,024813.640,A,3158.4608,N,11848.3737,E,10.05,324.27,150706,,,A*50字段0：$GPRMC，语句ID，表明该语句为Recommended Minimum Specific GPS/TRANSIT Data（RMC）推荐最⼩定位信息 字段1：UTC时间，hhmmss.sss格式字段2：状态，A=定位，V=未定位字段3：纬度ddmm.mmmm，度分格式（前导位数不⾜则补0）字段4：纬度N（北纬）或S（南纬）字段5：经度dddmm.mmmm，度分格式（前导位数不⾜则补0）字段6：经度E（东经）或W（西经）字段7：速度，节，Knots字段8：⽅位⾓，度字段9：UTC⽇期，DDMMYY格式字段10：磁偏⾓，（000 - 180）度（前导位数不⾜则补0）字段11：磁偏⾓⽅向，E=东W=西字段16：校验值$GPVTG例：$GPVTG,89.68,T,,M,0.00,N,0.0,K*5F字段0：$GPVTG，语句ID，表明该语句为Track Made Good and Ground Speed（VTG）地⾯速度信息字段1：运动⾓度，000 - 359，（前导位数不⾜则补0）字段2：T=真北参照系字段3：运动⾓度，000 - 359，（前导位数不⾜则补0）字段4：M=磁北参照系字段5：⽔平运动速度（0.00）（前导位数不⾜则补0）字段6：N=节，Knots字段7：⽔平运动速度（0.00）（前导位数不⾜则补0）字段8：K=公⾥/时，km/h字段9：校验值GPS编码格式及C语⾔解码GPS接收机只要处于⼯作状态就会源源不断地把接收并计算出的GPS导航定位信息通过串⼝传送到计算机中。

#include <stdio.h>#include <math.h>#include <stdlib.h>#define bGM84 3.986005e14#define bOMEGAE84 7.2921151467e-5void main(){long double roota=0.515365263176E+04; //轨道长半轴的平方根（根号a）long double toe=0.720000000000E+04; //观测时刻toelong double m0=-0.290282040486E+00; //参考时刻toe的平近点角long double e=0.678421219345E-02; //轨道偏心率elong double delta_n=0.451411660250E-08;//卫星的摄动改正数△nlong double smallomega=-0.258419417299E+01;//近地点角距ωlong double cus=0.912137329578E-05;//纬度幅角正弦调和项改正的振幅（弧度）long double cuc=0.189989805222E-06;//纬度幅角余弦调和项改正的振幅（弧度）long double crs=0.406250000000E+01;//轨道半径的余弦调和项改正的振幅（m）long double crc=0.201875000000E+03;//轨道半径的正弦调和项改正的振幅（m）long double cis=0.949949026108E-07;//轨道倾角的余弦调和项改正的振幅（弧度）long double cic=0.130385160446E-07;//轨道倾角的正弦调和项改正的振幅（弧度）long double idot=-0.253939149013E-09;//轨道倾角变化率Ilong double i0=0.958512160302E+00; //轨道倾角（弧度）long double bigomega0=-0.137835982556E+01;//升交点赤经long double earthrate=bOMEGAE84; //地球自转的速率welong double bigomegadot=-0.856928551657e-08;long double t=0.720000000000E+04; //加入卫星钟差改正的归化时间long double A;long double n0=0,n,tk;long double mk,ek,tak,ik,omegak,phik,uk,rk;long double corr_u,corr_r,corr_i;long double xpk,ypk,xk,yk,zk;int i;printf("输入的数据：\n");printf("√a=%e \n",roota);printf("toe=%e \n",toe);printf("e=%e \n",e);printf("i0=%e \n",i0);printf("ω=%e \n",smallomega);printf("△n=%e \n",delta_n);printf("Ω0=%e \n",bigomega0);printf("I=%e \n",idot);printf("Cuc=%e \n",cuc);printf("Cus=%e \n",cus);printf("Crc=%e \n",crc);printf("Crs=%e \n",crs);printf("Cic=%e \n",cic);printf("Cis=%e \n",cis);printf("\n\n输出的结果为：\n",e);A=roota*roota;n0=sqrt(bGM84/(A*A*A));//平均角速度n0printf("n0=%.10lf \n",n0);tk=t-toe;//相对于参考时刻toe的归化时间tkprintf("tk=%.10lf \n",tk);n=n0+delta_n;//加摄动改正后的卫星平均角速度printf(" n=%.10lf \n",n);mk=m0+n*tk;//卫星平近点角printf("mk=%.10lf \n",mk);ek=mk;for(i=0;i<10;i++) ek=mk+e*sin(ek);//利用迭代法求偏近点角ekprintf("ek=%.10lf \n",ek);tak=atan2(sqrt(1.0-e*e)*sin(ek)/cos(ek)-e);//真近点角Vk的计算printf("Vk=%.10lf \n",tak);phik=tak+smallomega;//升交距角φk的计算printf("φk=%.10lf \n",phik);corr_u=cus*sin(2.0*phik)+cuc*cos(2.0*phik);//升交距角u的摄动改正δu printf("δu=%.10lf \n",corr_u);corr_r=crs*sin(2.0*phik)+crc*cos(2.0*phik);//卫星矢量r的摄动改正δr printf("δr=%.10lf \n",corr_r);corr_i=cis*sin(2.0*phik)+cic*cos(2.0*phik);//轨道倾角i的摄动改正δi printf("δi=%.10lf \n",corr_i);uk=phik+corr_u;//升交距角uprintf("uk=%.10lf \n",uk);rk=A*(1.0-e*cos(ek))+corr_r;//卫星矢量rprintf("rk=%.10lf \n",rk);ik=i0+idot*tk+corr_i;//轨道倾角iprintf("ik=%.10lf \n",ik);xpk=rk*cos(uk);//卫星在轨道平面坐标系的坐标ypk=rk*sin(uk);printf("xpk=%.10lf \n",xpk);printf("ypk=%.10lf \n",ypk);omegak=bigomega0+(bigomegadot-earthrate)*tk-earthrate*toe;//升交点经度Ωk的计算printf("Ωk=%.10lf \n\n",omegak);xk=xpk*cos(omegak)-ypk*sin(omegak)*cos(ik);//地心固定坐标系的直角坐标yk=xpk*sin(omegak)+ypk*cos(omegak)*cos(ik);zk=ypk*sin(ik);printf("Xk=%.4lf \n",xk);printf("Yk=%.4lf \n",yk);printf("Zk=%.4lf \n",zk);}。