卫星位置计算程序
从广播星历计算卫星位置
从广播星历计算卫星位置： 1. 计算卫星运动的平均角速度n首先根据广播星历中给出的参数A 计算出参考时刻TOE 的平均角速度0n ：30)(A GM n =，式中，GM 为万有引力常数G 与地球总质量M 之积。
然后根据广播星历中给出的摄动参数n ∆计算观测时刻卫星的平均角速度n ：n n n ∆+=0。
2. 计算观测瞬间卫星的平近点角M :)(0TOE t n M M -+=式中，0M 为参考时刻TOE 时的平近点角，由广播星历给出。
3. 计算偏近点角E ：E e M E sin +=解上述方程可用迭代法或微分方程改正法。
4. 计算真近点角f⎪⎪⎩⎪⎪⎨⎧--=--=E e E e f Ee e Ef cos 1sin 1sin cos 1cos cos 2式中，e 为卫星轨道的偏心率，由广播星历给出。
5. 计算升交距角u '：f u +='ω式中，ω为近地点角距，由广播星历给出。
6. 计算摄动改正项i r u δδδ,,：广播星历中给出了下列6个摄动参数：is ic rs rc us uc C C C C C C ,,,,,，据此可以求出由于2J 项而引起的升交距角u 的摄动改正项u δ、卫星矢径r 的摄动改正项u δ和卫星轨道倾角i 的摄动改正项i δ。
计算公式如下：⎪⎩⎪⎨⎧'+'='+'='+'=u C u C u C u C u C u C is ic irs rc r us uc u 2sin 2cos 2sin 2cos 2sin 2cos δδδ 7. 计算0,,i r u ''进行摄动改正⎪⎪⎩⎪⎪⎨⎧-++=+-=+'=+'=)()cos 1(0TOE t dt di i i E a r r u u i r r u δδδδ 式中：a 为卫星轨道的长半径，2)(A a =，0i 为TOE 时刻的轨道倾角，由广播星历中的开普勒六参数给出，dtdi为i 的变化率，由广播星历中的摄动九参数给出。
GPS卫星坐标计算
GPS卫星坐标计算GPS（全球定位系统）是一种通过地球上的卫星提供位置和时间信息的导航系统。
GPS卫星坐标计算是指根据接收到的卫星信号来确定观测站位于球面上的位置。
GPS系统是由一组位于中轨道上的卫星组成，它们每天绕地球运行两次，以提供全球的覆盖范围。
每个卫星都携带有高精度的原子钟，用来产生精确的时间信号。
GPS接收机位于地面上，它接收到来自多颗卫星的信号，并测量信号的到达时间和卫星位置。
经过计算，接收机可以确定自身的空间坐标。
计算GPS卫星坐标的过程可以分为以下几个步骤：1.接收卫星信号：GPS接收机通过天线接收到来自多颗卫星的信号。
2.测量信号到达时间：接收机测量每个信号的到达时间，这需要精确的时钟。
由于GPS接收机一般没有原子钟那样的高精度时钟，所以需要利用接收到的卫星信号来校准本地时钟。
3.计算卫星位置：GPS接收机需要知道每颗卫星在接收时间点的准确位置。
每颗卫星通过广播自身的位置和时间信息，接收机可以根据接收到的信号来计算卫星的位置。
4.求解距离：接收机通过测量信号到达时间和卫星位置计算出距离。
由于信号的传播速度是已知的大约是光速，我们可以根据距离和到达时间计算出信号的传播时间。
5.根据接收到的信号来计算自身的位置。
接收机通过多个卫星信号的距离来确定自身的位置，这涉及到多种解算方法，例如最小二乘估计等。
接收机需要至少接收到四颗卫星的信号来解算自身的位置。
这些步骤涉及到大量的数学和物理计算，例如测量时间、测量距离、计算坐标等。
为了提高计算的精度，还需要考虑一些因素，例如信号传播时的大气延迟等。
总的来说，GPS卫星坐标计算是一项复杂而精确的工程，涉及到多个步骤和数学模型。
随着技术的不断进步，GPS定位的精度和可靠性也在不断提高，为导航、地球科学等领域的应用提供了重要的支持。
卫星导航定位算法与程序设计_第4课_伪距单点定位数学模型
5% 50% 95% 分布
5cm 10cm 20cm 50cm 1m 2m 5m 10m 20m 50m 100m
引用
定位精度与用户需求的关系
精密 工程 监控
地壳 形变 工程 监控
精密 大地 定位
工程与大 地定位
地理 信息 更新
地理 信息 更新
精密 交通 监控
交通 近海 自引导 监控 交通 导航
控制
卫星导航定位算法与程序设计
主讲： 刘晖 副教授
武汉大学卫星导航定位技术研究中心
课程讲授次序中的位置
软件设计相关 编程语言相关 定位算法相关
•建立起软件设计的宏观概念（战略层次）; •巩固软件设计方法的知识; •训练从设计到实现的工作流程;
•掌握软件编程平台的使用（战术层次）; •培养良好的编程习惯 •编程技巧的训练
GDOP – Geometry Dilution of Precision
PDOP – Position Dilution of Precision
TDOP – Time Dilution of Precision
HDOP – Horizontal Dilution of Precision
内容
GNSS定位技术回顾
定位技术的分类及发展 定位技术及其精度 定位精度与用户需求的关系 当前技术热点
伪距单点定位数学模型 伪距差分定位数学模型
概述
GNSS可提供全天候、高精度、高可用性、 高时效的三维空间定位
GNSS：3G+1C
3G：GPS、GLONASS、Galileo 1C：COMPASS
P ( X s X )2 (Ys Y )2 (Zs Z )2 atmos dt dT
卫星的位置计算.
• 计算思路
– 以地心为原点，建立轨道直角坐标系； – 计算卫星在轨道直角坐标系下的坐标； – 然后将轨道直角坐标系作一系列旋转，求出卫星在天球坐标 系下的坐标； – 将天球坐标系转换至地球坐标系。
无摄运动下卫星位置的计算
轨道直角坐标系
x0 r cos s f s y0 r sin s f s
u s f s
x0 r cos u y r sin u 0 z0 0
z0 0
轨道坐标系转换为天球坐标系
• 第一步：绕X’轴顺转角度i，以使Z’轴和Z轴重合 • 第二步：绕Z轴顺转角度Ω ，以使X’轴和X轴重合
x x0 y R R i y 1 0 3 z z0 cos sin 0 R3 sin cos 0 0 1 0 0 0 1 R1 i 0 cos i sin i 0 sin i cos i
天球坐标系转换至地球坐标系
• 起始子午面与过春分点的子午面的夹角实质为春分点的格 林尼治恒星时GAST。将瞬时天球坐标系绕Z轴逆转角度 GAST，即与瞬时地球坐标系重合。
X x Y R GAST y 3 Z z
cos GAST sin GAST 0 R3 GAST ቤተ መጻሕፍቲ ባይዱ sin GAST cos GAST 0 0 0 1
卫星导航定位算法与程序设计_第13课_卫星位置钟差计算
卫星位置的计算（2）
• 任意时刻t卫星位置的计算
– 原理：插值法 – 方法：拉格朗日插值法、切比雪夫多项式、三次样 条内插、三角多项式内插等 已知函数y f ( x)的n个结点x0 , x1 ,..., xn 及其对应的
函数值y0 , y1 ,..., yn 对于插值区间内的任一点x，其函数 值为 x xi f ( x) ( ) yk k 0 i 0 xk xi ik
xk rk cos uk yk rk sin uk
18
计算卫星位置(5/5)
14 计算改正后的升交点经度
Lk 0 earth t toe
15计算在地固坐标系下的位置
xk X Y RZ ( k ) Rx (ik ) yk Z 0
GM
e 7.2921151467 105 rad s
地球自转角速度
卫星星历
toe , A , e, M 0 , , i0 , n, i IDOT , Cus , Cuc , Crs , Crc , Cis , Cic toc , ClkBias, ClkDrift , ClkDriftRate
星历计算中常用常量和参数 星历参数详解
广播星历参数
精密星历参数
根据广播星历计算卫星位置等参数 根据精密星历计算卫星位置等参数 上机实习
星历计算中常用常数和符合
3.1415926535898
c 2.99792458 108 m s
真空中的光速
2
3.9860047 1014 m3 s
9
广播星历
A , e, M 0 , , i0 , toe , 轨道根数 n, i IDOT , Cus , Cuc , Crs , Crc , Cis , Cic 轨道摄动量
GPS卫星坐标计算
GPS卫星坐标计算
GPS系统由全球定位系统(GPS)组成，包括24颗运行在近地轨道上的
人造卫星，地面控制台和GPS接收器。
这些卫星以精确的轨道方式固定的
环绕着地球，它们通过无线电波将时间和位置信息传输到地面的GPS接收器。
具体步骤如下：
1.接收卫星信号：GPS接收器会接收到至少4颗卫星发出的信号。
这
些信号包括卫星的位置信息、时间戳和卫星信号的延迟。
2.计算信号传播时间差：接收器通过比较接收到的卫星信号和接收器
内部的原子钟产生的时间信号之间的差异，计算出信号传播的时间差。
3.确定接收器与卫星的距离：通过信号传播时间差和光速
（299,792,458米/秒），可以计算出接收器与卫星之间的距离。
公式为：距离=时间差x光速。
4.计算接收器的位置：通过接收到的至少4个卫星的距离信息，可以
计算出接收器相对于卫星的位置。
每个卫星会提供一个球面坐标，通过这
些球面坐标的交点，可以确定接收器的位置。
5.校准接收器的时间：接收器内部的原子钟会有一定的误差，因此需
要通过接收到的卫星信号的时间戳来校准接收器的时间。
6.确定地球的形状和尺寸：GPS系统还会考虑地球的形状和尺寸，以
便更精确地确定接收器的位置。
通过考虑地球的椭球形状、重力场和大气
层对卫星信号的影响，可以提高GPS定位的精确度。
总结起来，GPS卫星坐标计算的过程涉及接收卫星信号、计算信号传播时间差、确定接收器与卫星的距离、计算接收器的位置、校准接收器的时间以及考虑地球的形状和尺寸等步骤。
通过这些计算，可以精确测量地球上其中一点的位置坐标。
轨道卫星运动位置计算
1）计算平均角速度n
已知卫星轨道长半轴a，利用计算平均角速度。
2）计算平近点角M和偏近点角E
已知卫星过近地点时刻τ和卫星轨道离心率e，利用平近点交M和时间t的关系式计算平近点角M。
利用开普勒方程计算偏近点角E。
3）计算卫星向径的模r
利用式计算。
4）计算卫星真近点角f
利用计算。
5）计算卫星在轨道平面直角坐标系中的坐标（x’,y’）
利用式计算，其中r已计算出，如下
注意此处可以通过偏近点角E和椭圆参数直接计算坐标。
6）卫星在天球坐标系中的位置
由轨道倾角i，升交点赤经Ω和近升角距ω三个轨道参数，可以计算出卫星在天球坐标系中的位置。
此处R下标1、2、3对应的为x’、y’、z’的坐标轴，负号表示顺时针旋转。
7）卫星在瞬时地球坐标系中的位置
上面所求的坐标或速度，一般为惯性坐标系J2000.0。
要实现天球坐标系到地球坐标系的转换，应该首先考虑岁差和章动的影响先转换到瞬时真天球坐标系中。
但在实际应用中，如GPS导航电文的轨道根数提供的轨道根数，所求的上述结果已对应于瞬时天球坐标系，因而只需进行Z轴旋转GAST(t)，就转换到瞬时地球坐标系。
8）卫星在协议地球坐标系中的位置
若考虑极移的影响，有
以上就完成了卫星位置的计算。
GPS卫星坐标计算
卫星定位技术与方法--根据广播星历参数计算卫星坐标作业报告指导教师：熊永良班级：测绘二班学生姓名：段海东学生学号： 2 0 0 8 0 7 8 3作业日期：2010 年12月08 日目录一. 已知数据 (2)二. 计算步骤 (2)1.平均角速度 (mean angular speed) (2)2.规化时刻(normal time) (3)3.平近点角(mean anomaly) (3)4.偏近点角(eccentric anomaly) (3)5.真近点角(true anomaly) (3)6.升交距角(argument of ascending node) (3)7. 轨道向径（Orbital radius） (3)8. 扰动改正（Perturbed correction） (4)10.卫星在升交点轨道直角坐标系中的坐标 (4)11. 升交点经度（Longitude of ascending node） (5)三. 源程序 (5)四．程序运行结果 (14)七.作业体会 (15)根据广播星历参数计算卫星坐标一. 已知数据： 根据以下的广播星历参数计算UTC2004年1月30日8点0分00秒—20分00秒，每隔一分钟的PRN7的卫星坐标。
Compute the coordinate of PRN7 with interval of 1 minute.Navigation data:卫星导航文件格式：二. 计算步骤：The steps for satellite coordinates1.平均角速度 (mean angular speed)：∆n 由广播星历获得, GM=3.986005e+14 nn n ∆+=030a GM n =2.规化时刻(normal time)：t0已知（由广播星历获得），t 为GPS 周秒3.平近点角(mean anomaly)：M0已知（由广播星历获得）4.偏近点角(eccentric anomaly)：迭代求解：初始值取E=M ，以弧度为单位5.真近点角(true anomaly)：6.升交距角(argument of ascending node)：ω近地点角距(argument of perigee)7. 轨道向径（Orbital radius ）：k k t n M M ⋅+=0)cos 1(k k E e a r ⋅-⋅=e E E e V k kk --=cos sin 1arctan 2k k k E e M E sin ⋅+=0t t t k -=ωφ+=k 0V8. 扰动改正（Perturbed correction ）：• 升交角距（Argument of ascending node ）•• 轨道向径（Orbital radius ）• 轨道顷角（Orbital inclination ）是升交角距 (the argument of ascending node)9. 改正后升交角距、轨道向径、轨道倾角改正后升交角距（Corrected argument of ascending nod ）改正后的轨道向径（Corrected orbital radius)改正后的轨道倾角（Corrected orbital inclination ）10.卫星在升交点轨道直角坐标系中的坐标：如下图所示0φ00i 2sin 2cos φφδiS iC C C +＝00r 2sin 2cos φφδrS rC C C +＝00u 2sin 2cos φφδuS uC C C +＝rk k )cosE e 1(a r δ+⋅-⋅=u0k u δφ+=ki 0k t )IDOT (i i ++=δ11. 升交点经度（Longitude of ascending node ）：如下图所示12. 在地固坐标系中的卫星位置（Expressed in spheric coordinate system ）三. 源程序：using System;using System.Collections.Generic;using ponentModel;using System.Data;using System.Drawing;using System.Text;using System.Windows.Forms;using System.IO;using System.Text.RegularExpressions;kk k k k k sinu r y cosu r x ==s/rad 102921151467.7t t )(5-e 0e k e 0t ⨯=--Ω+Ω=ωωωλekk k k k k k k k k i y Z i y x Y i y x sin cos cos sin sin cos cos k k k k k =+=-=X λλλλnamespace Test{public partial class Form1 : Form{public Form1(){InitializeComponent();listView1.Columns.Add("序号", 40);listView1.Columns.Add("星历内容", 130);listView1.Columns.Add("导航数据", 130);listView1.GridLines = true; //显示表格线listView1.View = View.Details;//显示表格细节listView1.HeaderStyle = ColumnHeaderStyle.Clickable;//对表头进行设置listView2.Columns.Add("时间", 60);listView2.Columns.Add("x坐标", 150);listView2.Columns.Add("y坐标", 150);listView2.Columns.Add("z坐标", 150);listView2.GridLines = true; //显示表格线listView2.View = View.Details;//显示表格细节listView2.HeaderStyle = ColumnHeaderStyle.Clickable;//对表头进行设置}private void button1_Click(object sender, EventArgs e){//读取相对路径string str1 =AppDomain.CurrentDomain.SetupInformation.ApplicationBase;string filename = str1 + "navigation data.txt ";//读取卫星广播星历文件StreamReader myreader = new StreamReader(filename, Encoding.Default);string myinfo = myreader.ReadToEnd();myreader.Close();//把卫星广播星历里的D改为e，以便后续计算string mystring = myinfo;string myinfor = mystring;myinfor = mystring.Replace("D", "e");//把数据分开读入到一个数组中存储string[] split = new string[] { " " };string[] arrs = myinfor.Split(split,StringSplitOptions.RemoveEmptyEntries);double[] M = new double[arrs.Length];//星历代码string[] N = new string[arrs.Length];N[0] = "PRN"; N[1] = "Yer"; N[2] = "Mon"; N[3] = "day"; N[4] = "H"; N[5] = "M"; N[6] = "sec"; N[7] = "a0"; N[8] = "a1"; N[9] = "a2"; N[10] = "IODE"; N[11] = "Crs"; N[12] = "delta-n"; N[13] = "M0";N[14] = "Cuc"; N[15] = "e"; N[16] = "Cus"; N[17] = "sqrt(a)"; N[18] = "t0e"; N[19] = "Cic"; N[20] = "omega0"; N[21] = "Cis"; N[22] = "i0"; N[23] = "Crc"; N[24] = "omega"; N[25] = "omega-spot";N[26] = "IDOT"; N[27] = "Codes on L2 channel"; N[28] = "GPS Week"; N[29] = "L2 P data flag";N[30] = "SV accuracy"; N[31] = "SV health"; N[32] = "TGD"; N[33] = "IODC Issue of Data";N[34] = "Transmission time of message";//把卫星广播星历读入到数组中for (int i = 0; i < arrs.Length; i++){ListViewItem li = new ListViewItem();li.Text = (i + 1).ToString();li.SubItems.Add(N[i]);li.SubItems.Add(arrs[i]);listView1.Items.Add(li);M[i] = double.Parse(arrs[i]);}double t0e=460800.00;double t;double[] XK = new double[21];double[] YK = new double[21];double[] ZK = new double[21];for (int l = 0; l < 21;l++ ){//计算平均角速度double GM = 3986004.418e008;double n, n0;n0 = Math.Sqrt(GM / (Math.Pow(M[17], 6)));n = M[12] + n0;//规划时刻t=t0e+l*60;double tk = t-t0e;double Mk = M[13] + n * tk;// 迭代计算平近点角的计算double Ek, Ek1;Ek = Mk;Ek1 = Mk + M[15] * Math.Sin(Ek);do{Ek = Ek1;Ek1 = Mk + M[15] * Math.Sin(Ek);}while (Math.Abs(Ek1 - Ek) > 1e-15);//计算真近点角double Vk = Math.Atan(((Math.Sqrt(1 - M[15] * M[15]) * Math.Sin(Ek))) / (Math.Cos(Ek) - M[15]));//反正切值的象限处理double A1 = Math.Atan(((Math.Sqrt(1 - M[15] * M[15]) * Math.Sin(Ek))));double A2 = Math.Cos(Ek) - M[15];if (Vk < 0){if (A1 < 0 && A2 > 0){Vk += 2 * Math.PI;}if (A1 > 0 && A2 < 0){Vk += Math.PI;}}else{if (A1 <= 0 && A2 <= 0)Vk = Vk + Math.PI;}//计算升交角距double fk = Vk + M[24];//摄动改正值的计算double du = M[14] * Math.Cos(2 * fk) + M[16] * Math.Sin(2 * fk);double dr = M[23] * Math.Cos(2 * fk) + M[11] * Math.Sin(2 * fk);double di = M[19] * Math.Cos(2 * fk) + M[21] * Math.Sin(2 * fk);//摄动改正double uk = fk + du;double rk = Math.Pow(M[17], 2) * (1 - M[15] * Math.Cos(Ek))+ dr;double ik = M[22] + di + M[26] * tk;//计算卫星在升交点轨道直角坐标系的坐标double xk = rk * Math.Cos(uk);double yk = rk * Math.Sin(uk);//计算升交点经度double we = 7.2921151467e-05;double jdt = M[20] + (M[25] - we) * tk - we * M[18];//卫星在地固坐标系中的空间直角坐标XK[l] = xk * Math.Cos(jdt) - yk * Math.Cos(ik) * Math.Sin(jdt);YK[l] = xk * Math.Sin(jdt) + yk * Math.Cos(ik) * Math.Cos(jdt);ZK[l] = yk * Math.Sin(ik);}//输出卫星在地固坐标系中的空间直角坐标for (int i = 0; i < 21; i++){ListViewItem list = new ListViewItem();list.Text = "第" + i.ToString() + "分钟";list.SubItems.Add(XK[i].ToString());list.SubItems.Add(YK[i].ToString());list.SubItems.Add(ZK[i].ToString());listView2.Items.Add(list);}}private void button2_Click(object sender, EventArgs e) {//读取相对路径string str1 =AppDomain.CurrentDomain.SetupInformation.ApplicationBase;string filename = str1 + "navigation data.txt ";//读取卫星广播星历文件StreamReader myreader = new StreamReader(filename, Encoding.Default);string myinfo = myreader.ReadToEnd();myreader.Close();//把卫星广播星历里的D改为e，以便后续计算string mystring = myinfo;string myinfor = mystring;myinfor = mystring.Replace("D", "e");//把数据分开读入到一个数组中存储string[] split = new string[] { " " };string[] arrs = myinfor.Split(split, StringSplitOptions.RemoveEmptyEntries);//星历代码string[] N = new string[arrs.Length];N[0] = "PRN"; N[1] = "Yer"; N[2] = "Mon"; N[3] = "day"; N[4] = "H"; N[5] = "M"; N[6] = "sec"; N[7] = "a0"; N[8] = "a1"; N[9] = "a2"; N[10] = "IODE"; N[11] = "Crs"; N[12] = "delta-n"; N[13] = "M0";N[14] = "Cuc"; N[15] = "e"; N[16] = "Cus"; N[17] = "sqrt(a)"; N[18] = "t0e"; N[19] = "Cic"; N[20] = "omega0"; N[21] = "Cis"; N[22] = "i0"; N[23] = "Crc"; N[24] = "omega"; N[25] = "omega-spot";N[26] = "IDOT"; N[27] = "Codes on L2 channel"; N[28] = "GPS Week"; N[29] = "L2 P data flag";N[30] = "SV accuracy"; N[31] = "SV health"; N[32] = "TGD"; N[33] = "IODC Issue of Data";N[34] = "Transmission time of message";//把卫星广播星历文件内容输出到listView进行查看for (int i = 0; i < arrs.Length; i++){ListViewItem li = new ListViewItem();li.Text = (i + 1).ToString();li.SubItems.Add(N[i]);li.SubItems.Add(arrs[i]);listView1.Items.Add(li);}}}}四．程序运行结果：星历内容的读取、显示：卫星坐标计算结果：运行界面：作业体会：此次作业，收获颇多。
GPS卫星WGS-84坐标计算
������������ = ������0 + ������������ + i������������ (9) 计算卫星在轨道平面上的坐标 ������������ = ������������ ������������������������������ ������������ = ������������ ������������������������������ (10) 计算观测时刻 t 的升交点经度������������ ������������ = Ω0 + Ω − ������������ ������������ − ������������ ������������������ (11) 计算卫星在 WGS-84 坐标系中的坐标 ������������ = ������������ ������������������������������ − ������������ ������������������������������ ������������������������������ ������������ = ������������ ������������������������������ + ������������ ������������������������������ ������������������������������ ������������ = ������������ ������������������������������
下面是用户位置Байду номын сангаас算
估计位置坐标表示估计位置与真实位置的偏移量
设
GPS 卫星的广播星历提供 16 个星历参数,其中包括 1 个参考时刻、6 个相应参考时刻的开普 勒轨道参数和 9 个轨道摄动修正参数。 用广播星历参数计算任一时刻 t 的卫星位置的步骤如 下: (1) 计算卫星运行的平均角速度 n 卫星运行的平均角速度������ = ������0 + ∆������，式中������0 =
GPS卫星位置计算(C语言)
#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);}。
北斗导航卫星在轨位置计算数学模型
北斗导航卫星在轨位置计算数学模型下载提示：该文档是本店铺精心编制而成的，希望大家下载后，能够帮助大家解决实际问题。
文档下载后可定制修改，请根据实际需要进行调整和使用，谢谢！本店铺为大家提供各种类型的实用资料，如教育随笔、日记赏析、句子摘抄、古诗大全、经典美文、话题作文、工作总结、词语解析、文案摘录、其他资料等等，想了解不同资料格式和写法，敬请关注！Download tips: This document is carefully compiled by this editor. I hope that after you download it, it can help you solve practical problems. The document can be customized and modified after downloading, please adjust and use it according to actual needs, thank you! In addition, this shop provides you with various types of practical materials, such as educational essays, diary appreciation, sentence excerpts, ancient poems, classic articles, topic composition, work summary, word parsing, copy excerpts, other materials and so on, want to know different data formats and writing methods, please pay attention!引言北斗导航系统是我国自主研发的卫星导航定位系统，为了确保系统的准确性和可靠性，需要对其在轨位置进行精确计算。
GPS卫星定位坐标计算及程序设计
Ai X i li 0
(3-5)
对式（3-5）求解，便得到接收机地心坐标的唯一
解
X i Ai1li
4.程序设计
• 1、GPS时间转换程序 • 2、利用广播星历计算卫星坐标程序 • 3、地面点近似坐标计算程序
5.实例计算和精度分析
• 以2009年5月7日南京工业大学江浦校区控 制网20号控制点观测数据为例，来说明如 何利用该程序计算卫星坐标和地面点的近 似坐标。该数据利用华测GPS接收机观测， 观测时间为2小时。
• 3.新儒略日（Modified Julian Day-MJD）：从儒略 日中减去2400000.5天来得到，给出的是从1858年11 月17日子夜开始的天数。特点是数值比儒略日小。
• 4.年积日（Day Of Year-DOY）：从当前1月1日开始 的天数。
• 5.GPS时（GPS Time）：以1980年1月6日子夜为起点， 用周数和周内秒数来表示，为GPS系统内部计时法。
2.3GPS卫星的信号
• 导航电文 导航电文是包含有关卫星的星历、卫星工作状态 时间系统、卫星钟运行状态、轨道摄动改正、大 气折射改正和C/A码捕获P码等导航信息的数据码 （或D码），是利用GPS进行定位的数据基础。 导航电文的内容包括遥测码（TLW）、转换码 （HOW）、第一数据块、第二数据块和第三数据块 5部分。
RINEX数据格式
目前，RINEX格式已成为各厂商、学校、研究单 位在编制软件时采用的标准输入格式。RINEX格式 是纯ASCII码文本文件，共包含4个文件：
（1）观测数据文件：ssssdddf.yyo （2）导航文件：ssssdddf.yyn （3）气象数据文件：ssssdddf.yym （4）GLONASS数据文件：ssssdddf.yyg 其中：ssss——4个字母的测站名；
卫星导航定位算法PVT
卫星导航定位算法PVT1.PVT算法基本原理PVT算法的基本原理是通过接收到的卫星导航信号，测量接收机与每颗卫星之间的距离差异，并利用这些距离信息来计算接收机的位置、速度和时间。
2.卫星信号接收和检测PVT算法首先需要接收到至少四颗卫星的导航信号，每颗卫星将发送时间和位置等信息。
接收机通过接收和解调信号，获得卫星发送的导航数据，包括卫星的时钟偏差、频率偏差和星历信息等。
3.伪距测量PVT算法通过测量接收机与每颗卫星之间的时间差，即伪距来计算位置。
伪距的计算包含两个步骤：远程时钟校正和卫星-接收机距离计算。
-远程时钟校正：卫星发送的时间信息可能会受到发送时钟的误差影响。
接收机通过接收到的导航数据中的时间信息对其进行校正，使得接收到的时间信息与接收机自身的时间保持一致。
-距离计算：通过测量到达接收机的信号的传播时间，可以计算出卫星与接收机之间的距离。
由于其传播速度等于光速，可以通过接收到的时间差来计算卫星和接收机之间的距离。
4.定位解算PVT算法利用接收机到至少四颗卫星的距离信息，通过三角定位法计算出接收机的位置。
-三角测量法：基于距离的三角定位方法是最基本的定位算法之一、通过测量到多个卫星的距离，并将这些距离和卫星的位置信息转换为几何方程，可以利用三角法求解接收机的位置坐标。
-加权最小二乘法：由于测量误差和不同卫星的位置分布不均匀，可以采用加权最小二乘法对定位结果进行优化。
5.速度计算通过接收到的多个卫星的距离变化率，可以计算接收机的速度。
速度计算的基本原理是利用两次位置的差值和时间差来估计速度。
6.时间维护PVT算法不仅可以计算位置和速度，还可以计算接收机的时间。
通过接收到的卫星的时钟信息，可以校正接收机的时间偏差，并实时更新接收机的时间。
总结：卫星导航定位算法PVT通过接收到的卫星导航信号，并通过测量卫星-接收机之间的距离差异，实现对接收机位置、速度和时间的计算。
PVT 算法的主要步骤包括信号接收和检测、伪距测量、定位解算、速度计算和时间维护。
BDS卫星位置和速度
Created by Hui
X k xk cos k xk sin k k yk cos ik sin k yk cos ik sin k ik yk cos ik cos k k Yk xk sin k xk cos k k yk cos ik cos k yk sin ik sin k ik yk cos ik sin k k Z y sin i y sin i i k k k k k k
uk k uk rk A(1 e cos Ek ) rk i i IDOT t i k k k 0
xk rk cos uk rk sin uk uk yk rk sin uk rk cos uk uk
k e
1 e 2 sin Ek sin vk 1 e cos Ek cos v cos Ek e k uk Cus sin(2k ) Cuc cos(2k ) rk Crs sin(2k ) Crc cos(2k ) i C sin(2 ) C cos(2 ) is k ic k k uk k uk rk A(1 e cos Ek ) rk i i IDOT t i k k k 0
计算 GEO 卫星在自定义坐标系中 的速度
计算 GEO 卫星在 CGCS2000 坐标系 中的速度
0 0 RX ( ) 0 sin 0 cos sin RZ ( ) sin 0
cos sin cos 0 sin 0 0 0 0
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+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。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
#include<iostream>#include"math.h"#include"float.h"using namespace std;double kpro(double mk,double e){double f[2000],ek;int i=0;f[0]=mk;do{f[i+1]=mk+e*sin(f[i]);i++;}while(abs(f[i+1]-f[i])>=0.005);ek=f[i+1];return ek;}#define PI 3.14159265357589int main(void){doubletoc,a2,a1,a0,dotn,t1,m0,e,roota,toe,cic,crc,cis,crs,cuc,cus,omega0,omega,i0,omegadot,idot;double n,n0,dott,t,tk,mk,ek,vk,ak,u,r,i,uk,rk,ik,xk,yk,wk,we,x,y,z;int k,j;cout<<"输º?入¨?秒?数ºy据Y参?数ºy toc\n";cin>>toc;cout<<"输º?入¨?秒?数ºy据Y幂Y数ºy\n";cin>>k;toc=toc*pow(10,k);cout<<"输º?入¨?时º¡À钟¨®修T正y2 a2\n";cin>>a2;cout<<"输º?入¨?时º¡À钟¨®修T正y1 a1\n";cin>>a1;cout<<"输º?入¨?时º¡À钟¨®修T正y0 a0\n";cin>>a0;cout<<"输º?入¨?摄¦?动¡¥改?正y数ºy dn\n";cin>>dotn;cout<<"输º?入¨?观?测a时º¡À间? t'\n";cin>>t1;cout<<"输º?入¨?反¤¡ä常¡ê m0\n";cin>>m0;cout<<"输º?入¨?偏?心?率¨º e\n";cin>>e;cout<<"输º?入¨?长¡è半ã?径?平?方¤?根¨´ roota\n";cin>>roota;cout<<"输º?入¨?星?历¤¨²基¨´准Á?时º¡À间?系¦Ì数ºy toe\n";cin>>toe;cout<<"输º?入¨?星?历¤¨²基¨´准Á?时º¡À间?幂Y数ºy\n";cin>>j;toe=toe*pow(10,j);cout<<"输º?入¨?对?倾?角?余®¨¤弦¨©校¡ê正y值¦Ì cic\n";cin>>cic;cout<<"输º?入¨?对?轨¨¬道̨¤半ã?径?余®¨¤弦¨©校¡ê正y值¦Ì crc\n";cin>>crc;cout<<"输º?入¨?对?倾?角?正y弦¨©校¡ê正y值¦Ì cis\n";cin>>cis;cout<<"输º?入¨?对?轨¨¬道̨¤半ã?径?正y弦¨©校¡ê正y值¦Ì crs\n";cin>>crs;cout<<"输º?入¨?对?纬3度¨¨幅¤¨´角?余®¨¤弦¨©校¡ê正y值¦Ì cuc\n";cin>>cuc;cout<<"输º?入¨?对?纬3度¨¨幅¤¨´角?正y弦¨©校¡ê正y值¦Ì cus\n";cin>>cus;cout<<"输º?入¨?升¦y交?点Ì?经-度¨¨ omega0\n";cin>>omega0;cout<<"输º?入¨?近¨¹地Ì?点Ì?幅¤¨´角? omega\n";cin>>omega;cout<<"输º?入¨?倾?角? i0\n";cin>>i0;cout<<"输º?入¨?升¦y交?点Ì?经-度¨¨变À?化¡¥率¨º omegadot\n";cin>>omegadot;cout<<"输º?入¨?倾?角?率¨º Idot\n";cin>>idot;//计?算?平?均¨´角?速¨´度¨¨n0=sqrt(3.986005)*pow(10,7)/pow(roota,3);n=n0+dotn;cout<<"平?均¨´角?速¨´度¨¨是º?n：êo"<<n<<endl;//计?算?归¨¦化¡¥时º¡À间?dott=a0+a1*(t1-toc)+a2*(t1-toc)*a2*(t1-toc);t=t1-dott;tk=t-toe;if(tk>302400){tk=tk-604800;}else if(tk<-302400){tk=tk+604800;}cout<<"归¨¦化¡¥时º¡À间?是º?tk：êo"<<tk<<endl;//计?算?平?近¨¹点Ì?角?mk=m0+n*tk;cout<<"平?近¨¹点Ì?角?是º?mk：êo"<<mk<<endl;//计?算?偏?近¨¹点Ì?角?ek=kpro(mk,e);cout<<"偏?近¨¹点Ì?角?是º?ek "<<ek<<endl;//计?算?真?近¨¹点Ì?角?vk=atan((sqrt(1-e*e)*sin(ek))/(cos(ek)-e));if(cos(vk)>0 && sin(vk)>0 )vk=atan(sin(vk)/cos(vk));else if(sin(vk)>0 && cos(vk)<0)vk=acos(cos(vk));else if(sin(vk)<0 &&cos(vk)<0)vk=PI-asin(sin(vk));elsevk=2*PI-acos(cos(vk));cout<<"真?近¨¹点Ì?角?是º?vk：êo"<<vk<<endl;//计?算?距¨¤角?ak=vk+omega;cout<<"距¨¤角?是º?ak：êo"<<ak<<endl;//计?算?摄¦?动¡¥改?正y项?u=cuc*cos(2*ak)+cus*sin(2*ak);r=crc*cos(2*ak)+crs*sin(2*ak);i=cic*cos(2*ak)+cis*sin(2*ak);cout<<"摄¦?动¡¥改?正y项?:"<<"u:"<<u<<",r:"<<r<<",i:"<<i<<endl;//计?算?摄¦?动¡¥改?正y的Ì?升¦y交?距¨¤角?，ê?卫¨¤星?矢º?量¢?，ê?轨¨¬道̨¤倾?角?uk=ak+u;rk=roota*roota*(1-e*cos(ek))+r;ik=i0+tk*idot+i;cout<<"摄¦?动¡¥改?正y的Ì?升¦y交?距¨¤角?uk:"<<uk<<",卫¨¤星?矢º?量¢?rk:"<<rk<<",轨¨¬道̨¤倾?角?ik:"<<ik<<endl;//计?算?卫¨¤星?在¨²轨¨¬平?面?坐Á?标À¨º系¦Ì的Ì?坐Á?标À¨ºxk=rk*cos(uk);yk=rk*sin(uk);cout<<"卫¨¤星?在¨²轨¨¬平?面?坐Á?标À¨º系¦Ì的Ì?坐Á?标À¨º：êo"<<"xk:"<<xk<<",yk:"<<yk<<endl;//计?算?观?测a时º¡À刻¨¬升¦y交?点Ì?经-度¨¨we=7.292115147*pow(10,-5);wk=omega0+(omegadot-we)*tk-we*toe;cout<<"观?测a时º¡À刻¨¬升¦y交?点Ì?经-度¨¨：êo"<<wk<<endl;//计?算?卫¨¤星?在¨²低̨ª固¨¬系¦Ì下?坐Á?标À¨ºx=xk*cos(wk)-yk*cos(ik)*sin(wk);y=xk*sin(wk)+yk*cos(ik)*cos(wk);z=yk*sin(ik);cout<<"卫¨¤星?在¨²地Ì?固¨¬系¦Ì下?坐Á?标À¨º是º?：êo"<<"x:"<<x<<",y:"<<y<<",z:"<<z<<endl;return 0;}///////固定值版#include<iostream>#include"math.h"using namespace std;//迭̨¹代䨲法¤¡§计?算?偏?近¨¹点Ì?角?double kpro(double mk,double e){double f,t,m,ek;f=mk;do{m=f;t=mk+e*sin(m);f=t;}while((f-m)>10e-12);ek=f;return ek;}#define PI 3.141592655358979323846int main(void){//设¦¨¨置?参?数ºydoubletoc,a2,a1,a0,dotn,t1,m0,e,roota,toe,cic,crc,cis,crs,cuc,cus,omega0,omega,i0,omegadot,idot;double n,n0,dott,t,tk,mk,ek,vk,ak,u,r,i,uk,rk,ik,xk,yk,wk,we,x,y,z;//参?数ºy赋3值¦Ìtoc=1.801672*pow(10,-38);a2=0;a1=1.7053025658199999e-012;a0=0.00028600683435800001;dotn=4.1115998360700002e-009;t1=11404800.0;m0=1.2263973009600000;e=0.0053100715158500003;roota=5153.7127704599998;toe=0.0;cic=-9.6857547759999998e-008;crc=175.34375000000000;cis=-7.8231096267699997e-008;crs=-105.43750000000000;cuc=-5.5264681577700003e-006;cus=1.1192634701700000e-005;omega0=-2.9080127721900002;omega=-1.6819446292500000;i0=0.97432927738800001;omegadot=-7.7124641122299999e-009;ido t=1.8643633724999999e-010;//计?算?平?均¨´角?速¨´度¨¨n0=sqrt(3.986005)*pow(10,7)/pow(roota,3);n=n0+dotn;cout<<"平?均¨´角?速¨´度¨¨是º?n：êo"<<n<<endl;//计?算?归¨¦化¡¥时º¡À间?dott=a0+a1*(t1-toc)+a2*(t1-toc)*a2*(t1-toc);t=t1-dott;tk=t-toe;if(tk>302400){tk=tk-604800;}else if(tk<-302400){tk=tk+604800;}cout<<"归¨¦化¡¥时º¡À间?是º?tk：êo"<<tk<<endl;//计?算?平?近¨¹点Ì?角?mk=m0+n*tk;cout<<"平?近¨¹点Ì?角?是º?mk：êo"<<mk<<endl;//计?算?偏?近¨¹点Ì?角?ek=kpro(mk,e);cout<<"偏?近¨¹点Ì?角?是º?ek "<<ek<<endl;//计?算?真?近¨¹点Ì?角?vk=atan((sqrt(1-e*e)*sin(ek))/(cos(ek)-e));if(cos(vk)>0 && sin(vk)>0 )vk=atan(sin(vk)/cos(vk));else if(sin(vk)>0 && cos(vk)<0)vk=acos(cos(vk));else if(sin(vk)<0 &&cos(vk)<0)vk=PI-asin(sin(vk));elsevk=2*PI-acos(cos(vk));cout<<"真?近¨¹点Ì?角?是º?vk：êo"<<vk<<endl;//计?算?距¨¤角?ak=vk+omega;cout<<"距¨¤角?是º?ak：êo"<<ak<<endl;//计?算?摄¦?动¡¥改?正y项?u=cuc*cos(2*ak)+cus*sin(2*ak);r=crc*cos(2*ak)+crs*sin(2*ak);i=cic*cos(2*ak)+cis*sin(2*ak);cout<<"摄¦?动¡¥改?正y项?:"<<"u:"<<u<<",r:"<<r<<",i:"<<i<<endl;//计?算?摄¦?动¡¥改?正y的Ì?升¦y交?距¨¤角?，ê?卫¨¤星?矢º?量¢?，ê?轨¨¬道̨¤倾?角?uk=ak+u;rk=roota*roota*(1-e*cos(ek))+r;ik=i0+tk*idot+i;cout<<"摄¦?动¡¥改?正y的Ì?升¦y交?距¨¤角?uk:"<<uk<<",卫¨¤星?矢º?量¢?rk:"<<rk<<",轨¨¬道̨¤倾?角?ik:"<<ik<<endl;//计?算?卫¨¤星?在¨²轨¨¬平?面?坐Á?标À¨º系¦Ì的Ì?坐Á?标À¨ºxk=rk*cos(uk);yk=rk*sin(uk);cout<<"卫¨¤星?在¨²轨¨¬平?面?坐Á?标À¨º系¦Ì的Ì?坐Á?标À¨º：êo"<<"xk:"<<xk<<",yk:"<<yk<<endl;//计?算?观?测a时º¡À刻¨¬升¦y交?点Ì?经-度¨¨we=7.292115147*pow(10,-5);wk=omega0+(omegadot-we)*tk-we*toe;cout<<"观?测a时º¡À刻¨¬升¦y交?点Ì?经-度¨¨：êo"<<wk<<endl;//计?算?卫¨¤星?在¨²低̨ª固¨¬系¦Ì下?坐Á?标À¨ºx=xk*cos(wk)-yk*cos(ik)*sin(wk);y=xk*sin(wk)+yk*cos(ik)*cos(wk);z=yk*sin(ik);cout<<"卫¨¤星?在¨²地Ì?固¨¬系¦Ì下?坐Á?标À¨º是º?：êo"<<"x:"<<x<<",y:"<<y<<",z:"<<z<<endl;return 0;}。