- VisualStudio2022插件的安装及使用-编程手把手系列文章
- pprof-在现网场景怎么用
- C#实现的下拉多选框,下拉多选树,多级节点
- 【学习笔记】基础数据结构:猫树
本代码使用numpy,pandas,gnss_lib_py,matplotlib四个函数库,请提前安装.
#两条命令根据使用环境进行选择
pip install gnss-lib-py pandas matplotlib #Python环境安装代码
conda install gnss-lib-py pandas matplotlib -c conda-forge #conda环境安装代码
GitHub主页:https://github.com/Stanford-NavLab/gnss_lib_py?tab=readme-ov-file 文档主页:https://gnss-lib-py.readthedocs.io/en/latest/index.html 本文主要使用该库的读取以及转化为DataFrame功能,其中参数的命名规则以及时间转换规则可以在文档中找到.
建议使用jupyter进行执行,下方代码为分单元块格式,如果没有jupyter环境可以直接粘贴到一个python文件进行运行.
import gnss_lib_py as glp
import datetime
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
# 导入23n文件
file_path = 'brdc2550.23n'
data = glp.RinexNav(file_path)
data_df = data.pandas_df()
# 寻找最小差值的参考时刻
def find(inweekmilli, refers, insv):
filter_sv = refers[refers['gnss_sv_id'] == insv]
defference = np.abs(inweekmilli - filter_sv['t_oe'])
return defference.idxmin()
times = np.array([None] * 24 * 4)
gpsmillis = np.array([None] * 24 * 4)
n = 0
for hour in range(0, 24):
minut = 0
while minut < 60:
times[n] = datetime.datetime(2023, 9, 12, hour, minut, 0, tzinfo=datetime.timezone.utc)
gpsmillis[n] = glp.datetime_to_gps_millis(times[n])
minut += 15
n += 1
gpsmillis = np.array(gpsmillis)
GM=3.986005E+14
sqrtGM = np.sqrt(GM)
sv_list = [f'G{str(i).zfill(2)}' for i in range(1, 33)]
outdata = pd.DataFrame(columns=['data', 'gnss_sv_id', 'X', 'Y', 'Z'], index=range(24 * 4 * 32))
orbit = pd.DataFrame(columns=['data', 'gnss_sv_id', 'x', 'y'], index=range(24 * 4 * 32))
m = 0
j = 0
print("正在计算,请稍候。")
for gpsmilli in gpsmillis:
for sv in sv_list:
week, milli_week = glp.gps_millis_to_tow(gpsmilli)
milli_week=milli_week-18
index = find(milli_week, data_df, sv)
print(f'sv:{sv},time:{times[j]},index:{index}')
print(milli_week,data_df.iloc[index]['t_oe'])
a = np.power(data_df.iloc[index]['sqrtA'], 2)
n0 = sqrtGM / np.power(a, 3 / 2)
n = n0 + data_df.iloc[index]['deltaN']
tk = milli_week - data_df.iloc[index]['t_oe']
M = data_df.iloc[index]['M_0'] + n * tk
e = data_df.iloc[index]['e']
# 打印中间结果M和e
print(f"M: {M}, e: {e}")
# 解开普勒方程
E = M
for _ in range(50): # 使用迭代方法求解E
E = M + e * np.sin(E)
# 打印中间结果E
print(f"E: {E}")
f = np.arctan((np.sqrt(1 - e**2) * np.sin(E)) / (np.cos(E) - e))
if E > np.pi*0.5:
f=f+np.pi
if E < -np.pi*0.5:
f=f-np.pi
if np.pi*0.5 > E > 0 > f:
f=f+np.pi
if -np.pi*0.5 < E < 0 < f:
f=f-np.pi
print(f"arctan({(np.sqrt(1 - e**2) * np.sin(E)) / (np.cos(E) - e)}),f:{f}")
u_pie = data_df.iloc[index]['omega'] + f
r_pie = a * (1 - e * np.cos(E))
C_uc = data_df.iloc[index]['C_uc']
C_us = data_df.iloc[index]['C_us']
C_rc = data_df.iloc[index]['C_rc']
C_rs = data_df.iloc[index]['C_rs']
C_ic = data_df.iloc[index]['C_ic']
C_is = data_df.iloc[index]['C_is']
delta_u = C_uc * np.cos(2 * u_pie) + C_us * np.sin(2 * u_pie)
delta_r = C_rc * np.cos(2 * u_pie) + C_rs * np.sin(2 * u_pie)
delta_i = C_ic * np.cos(2 * u_pie) + C_is * np.sin(2 * u_pie)
u = u_pie + delta_u
r = r_pie + delta_r
i = data_df.iloc[index]['i_0'] + delta_i + data_df.iloc[index]['IDOT'] * tk
print(f'u:{u}')
x = r * np.cos(u)
y = r * np.sin(u)
w_e = 7.292115E-5
L = data_df.iloc[index]['Omega_0'] + (data_df.iloc[index]['OmegaDot']- w_e )* milli_week - data_df.iloc[index]['OmegaDot']*data_df.iloc[index]['t_oe']
X = x * np.cos(L) - y * np.cos(i) * np.sin(L)
Y = x * np.sin(L) + y * np.cos(i) * np.cos(L)
Z = y * np.sin(i)
orbit.iloc[m,:] = [times[j],sv,x,y]
outdata.iloc[m, :] = [times[j], sv, X, Y, Z]
m += 1
j += 1
print("由于结果较长,请到Excel中查看,文件位于代码同级目录下outdata.csv。")
outdata.to_csv('outdata.csv')
print("导出成功。")
# 三维坐标可视化显示
out_sv = 'G20'
fig = plt.figure()
ax = plt.axes(projection='3d')
X = outdata[outdata['gnss_sv_id'] == out_sv]['X']
Y = outdata[outdata['gnss_sv_id'] == out_sv]['Y']
Z = outdata[outdata['gnss_sv_id'] == out_sv]['Z']
ax.plot(X, Y, Z, label=out_sv)
ax.legend()
plt.show()
最后此篇关于使用广播星历计算卫星坐标(Python)的文章就讲到这里了,如果你想了解更多关于使用广播星历计算卫星坐标(Python)的内容请搜索CFSDN的文章或继续浏览相关文章,希望大家以后支持我的博客! 。
SQL 和一般开发的新手,我有一个表(COUNTRIES),其中包含字段(INDEX、NAME、POPULATION、AREA) 通常我添加一个客户端(Delphi)计算字段(DENSITY)和 On
我想使用 calc(100%-100px),但在我的 demo 中不起作用由于高度只接受像素,因此如何将此百分比值转换为像素。 最佳答案 以下将为您提供高度: $(window).height();
我正在尝试在 MySQL 中添加列并动态填充其他列。 例如我有一张表“数字”并具有第 1 列、第 2 列、第 3 列,这些总数应填充在第 4 列中 最佳答案 除非我误解了你的问题,否则你不只是在寻找:
我想返回简单计算的结果,但我不确定如何执行此操作。我的表格如下: SELECT COUNT(fb.engineer_id) AS `total_feedback`, SUM(fb.ra
我一直在尝试做这个程序,但我被卡住了,我仍然是一个初学者,任何帮助将不胜感激。我需要程序来做 打印一个 10 X 10 的表格,其中表格中的每个条目都是行号和列号的总和 包含一个累加器,用于计算所有表
这个计算背后一定有一些逻辑。但我无法得到它。普通数学不会导致这种行为。谁能帮我解释一下原因 printf ("float %f\n", 2/7 * 100.0); 结果打印 1.000000 为什么会
我想计算从 0 到 (n)^{1/2} - 1 的数字的 AND每个数字从 0 到 (n)^{1/2} - 1 .我想在 O(n) 中执行此操作时间,不能使用 XOR、OR、AND 运算。 具体来说,
如何在 Excel 中将公式放入自定义数字格式?例如(出于说明目的随机示例), 假设我有以下数据: 输入 输出 在不编辑单元格中的实际数据的情况下,我想显示单元格中的值除以 2,并保留两位小数: 有没
每次我在 Flutter 应用程序中调用计算()时,我都会看到内存泄漏,据我所知,这基本上只是一种生成隔离的便捷方法。我的应用程序内存占用增加并且在 GC 之后永远不会减少。 我已将我的代码简化为仅调
我有数字特征观察 V1通过 V12用于目标变量 Wavelength .我想计算 Vx 之间的 RMSE列。数据格式如下。 每个变量“Vx”以 5 分钟的间隔进行测量。我想计算所有 Vx 变量的观测值
我正在寻找一种使用 C 语言计算文件中未知字符数的简单方法。谢谢你的帮助 最佳答案 POSIX 方式(可能是您想要的方式): off_t get_file_length( FILE *file ) {
我正在使用 Postgres,并且我正试图围绕如何在连续日期跨度中得出第一个开始日期的问题进行思考。例如 :- ID | Start Date | End Date =================
我有一个订单表格,我在其中使用 jQuery 计算插件来汇总总数。 此求和工作正常,但生成的“总和”存在问题。总之,我希望用逗号替换任何点。 代码的基础是; function ($this) {
我在使用 double 变量计算简单算术方程时遇到问题。 我有一个具有 double 属性 Value 的组件,我将此属性设置为 100。 然后我做一个简单的减法来检查这个值是否真的是 100: va
我在这里看到了一些关于 CRC 32 计算的其他问题。但没有一个让我满意,因此是这样。 openssl 库是否有任何用于计算 CRC32 的 api 支持?我已经在为 SHA1 使用 openssl,
当我在PHP日期计算中遇到问题时,我感到惊讶。 $add = '- 30 days'; echo date('Y-m-01', strtotime($add)); // result is 2017-
我正在使用 javascript 进行练习,我编写了这个脚本来计算 2 个变量的总和,然后在第三个方程中使用这个总和!关于如何完成这项工作的任何想法都将非常有用! First Number:
我有一个来自EAC的提示单和一个包含完整专辑的FLAC文件。 我正在尝试制作一些python脚本来播放文件,因为我需要能够设置在flac文件中开始的位置。 如何从CueSheet格式MM:SS:FF转
这个问题已经有答案了: Adding two numbers concatenates them instead of calculating the sum (24 个回答) 已关闭去年。 我有一个
4000 我需要上面字段 name="quantity" 和 id="price" 中的值,并使用 javascript 函数进行计算,并将其显示在字段 id= 中仅当我单击计算按钮时才显示“总
我是一名优秀的程序员,十分优秀!