天和空间站凌月凌日时间计算的尝试

2023-10-11 21:30
文章标签 计算 时间 尝试 凌日 凌月 空间站

本文主要是介绍天和空间站凌月凌日时间计算的尝试，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

最近对空间站过境、凌日月时间计算产生了兴趣，之前写了空间站过境的代码，接着该空间站凌月时间计算了。在网上搜索没有现成的，暂时也没找到相关算法的介绍。网上找到只有这个可以查询国际空间站凌日月时间的网站。找不到就自己尝试办法吧。凌月无非就是在地面观察点看到两者的位置重合了，两条思路：1.空间站和月球在天球的赤经和赤纬差值极小即可，2.两者从观察点看的夹角为0，根据这两条思路尝试写代码。首先可以把查询时间段内能过境本地的时间筛选出来，这个已经在空间站过境时间计算的代码里实现，并且可已简化计算，只要把过境最高点大于10°（太低了看不到）的升降到地平线时的时间算出来即可。然后分别计算各时间段内两者赤经赤纬差值最小值或夹角最小值。值越小说明两者距离越近。
下面是自己尝试写的代码，代码比较乱，没怎么整理，计算后在stellarium上模拟。（编程和天文都是纯业余爱好，轻喷）

"""
使用skyfield计算空间站凌日月时间
delcomp 2021-06-25
"""
# Topos 已弃用 建议使用wgs84
from skyfield.api import load, wgs84
from skyfield.searchlib import find_maxima, find_minima, _find_discrete, find_discrete
from skyfield.trigonometry import position_angle_of
from skyfield.constants import tau, DAY_S
from pytz import timezone
import numpy as np
import os
import time
from datetime import datetime, timedelta
from my_function import BD09_to_GCJ02, LatLng2DMS, GCJ02_to_GPS84, BD09_to_WGS84start_time = time.perf_counter()# 百度地图经纬度
LAT = xx.xx1
LNG = xxx.xxSTATION = 'TIANHE'
# STATION = 'ISS (ZARYA)'
DAYS = 300
ELEVATION = 0  # 观察点高度
ALT_DEF = 10  # 指定的最低高度
DT_FSTR = '%Y-%m-%d %H:%M:%S'  # 时间格式化串beijing = timezone('Asia/Shanghai')  # 中国时区
stations_url = 'http://celestrak.com/NORAD/elements/stations.txt'
stations_file = 'stations.txt'
ts = load.timescale()
# satellites = load.tle_file(stations_url, reload=True)
if not os.path.exists(stations_file):print('本地文件不存在，下载stations文件...')satellites = load.tle_file(stations_url, reload=True)
else:satellites = load.tle_file(stations_file)print('加载本地stations文件！')# 根据卫星名称构建卫星字典by_name = {sat.name: sat for sat in satellites}tianhe = by_name[STATION]# 文件旧了，更新t = ts.now()days = t - tianhe.epochprint('{:.3f} days away from epoch'.format(days))if abs(days) > 0.5:satellites = load.tle_file(stations_url, reload=True)print('stations文件已更新！')# 根据卫星编号构建卫星字典
# by_number = {sat.model.satnum: sat for sat in satellites}
# # 25544：ISS
# # 48274：TIANHE
# satellite = by_number[25544]by_name = {sat.name: sat for sat in satellites}
tianhe = by_name[STATION]print(tianhe)
print('BD09坐标: %.5f %.5f' % (LAT, LNG))
# 转为GCJ02坐标
lat_gcj, lng_gcj = BD09_to_GCJ02(LAT, LNG)
lat_d, lat_m, lat_s = LatLng2DMS(lat_gcj)
lon_d, lon_m, lon_s = LatLng2DMS(lng_gcj)
print('GCJ02坐标: %.5f %.5f' % (lat_gcj, lng_gcj),'(%d°%d′%.1f″  %d°%d′%.1f″)' % (lat_d, lat_m, lat_s, lon_d, lon_m, lon_s))
# 转为WGS84
lat_wgs, lng_wgs = BD09_to_WGS84(LAT, LNG)
lat_d, lat_m, lat_s = LatLng2DMS(lat_wgs)
lon_d, lon_m, lon_s = LatLng2DMS(lng_wgs)
print('WGS84坐标: %.5f %.5f' % (lat_wgs, lng_wgs),'(%d°%d′%.1f″  %d°%d′%.1f″)' % (lat_d, lat_m, lat_s, lon_d, lon_m, lon_s))mypos = wgs84.latlon(lat_wgs, lng_wgs, elevation_m=ELEVATION)
eph = load('de/de421.bsp')
sun, earth, moon = eph['sun'], eph['earth'], eph['moon']# 使用本地时区
now = beijing.localize(datetime.now())
# t0 = now.replace(year= 2022, month=1, day=14, hour=0, minute=0, second=0, microsecond=0)
t0 = now.replace(hour=0, minute=0, second=0, microsecond=0)
t1 = t0 + timedelta(days=DAYS)
print('查询时间间隔(北京)：%s 至 %s' % (t0, t1))
# find_events参数为UTC，将北京时间转为UTC
t0 = ts.from_datetime(t0)
t1 = ts.from_datetime(t1)
print('查询时间间隔(UTC)：%s 至 %s' % (t0.utc_strftime(DT_FSTR), t1.utc_strftime(DT_FSTR)))
difference= tianhe - myposdef cheat(t):"""Avoid computing expensive values that cancel out anyway."""t.gast = t.tt * 0.0t.M = t.MT = np.identity(3)def below_horizon_at(t):cheat(t)return difference.at(t).altaz()[0].degrees < ALT_DEFdef altitude_at(t):cheat(t)return difference.at(t).altaz()[0].degrees# 找到符合条件（中天高度至少10度）的过境时间
half_second = 0.5 / DAY_S
orbits_per_minute = tianhe.model.no_kozai / tau
orbits_per_day = 24 * 60 * orbits_per_minute
step_days = 0.05 / max(orbits_per_day, 1.0)
if step_days > 0.25:step_days = 0.25
# step_days = 1 / DAY_S
altitude_at.step_days = step_days
# 找到最高点的时间及高度
tmax, altitude = find_maxima(t0, t1, altitude_at, half_second, 12)# 筛选高度>10的tmax 对应_find_discrete(高度<10)，
tmax = tmax.tt[altitude >= ALT_DEF]doublets_0 = np.repeat(np.concatenate(((t0.tt,), tmax, (t1.tt,))), 2)
jdo_0 = (doublets_0[:-1] + doublets_0[1:]) / 2.0# 在地平线的时间，一次升起(event=0)，一次落下(event=1)
# 注意：可能碰到一种特殊情况，在指定的时间段内，有可能第一个时间为降落至10°，即event=1 例如：[1,0,1,0,1,0,1,0,1]
# 或者最后一个时间为升高至10°，即event=0，例如：[0,1,0,1,0,1,0,1,0]
# 当上述两种情况，不同时出现时，在后面计算pass_pairs时，无法配对过境时间
# 如果两种情况同时存在时，[1,0,1,0,1,0,1,0,1,0]，可以配对过境时间(升-降，0-1)，但配对将是错误的, [1-0, 1-0, 1-0, ... ]
# 正确配对应该是[0-1, 0-1, 0-1, ... ]
trs_0, event_0 = _find_discrete(t0.ts, jdo_0, below_horizon_at, half_second, 8)
# 判读是否出现上面注释提到的特殊情况，删除
if event_0[0] == 1:  # 第一种情况trs_0 = trs_0[1:]
if event_0[-1] == 0:  # 第二种情况trs_0 = trs_0[:-1]
# 将符合条件的过境时间配对
pass_pairs = trs_0.reshape(len(trs_0) // 2, 2)
# exit(-1)
# 计算在这些时间区间内是否凌月
pos = earth + mypos# def calc_transit(t):
#     m = pos.at(t).observe(moon).apparent()
#     sat = pos.at(t).observe(earth + tianhe).apparent()
#     return position_angle_of(m.altaz(), sat.altaz()).degrees# 角距
def calc_transit_moon(t):m = pos.at(t).observe(moon)sat = pos.at(t).observe(earth + tianhe)return m.separation_from(sat).degreesdef calc_transit_sun(t):m = pos.at(t).observe(sun)sat = pos.at(t).observe(earth + tianhe)return m.separation_from(sat).degrees# 赤经赤纬差
# def calc_transit_1(t):
#     moon_ra, moon_dec, _ = pos.at(t).observe(moon).apparent().radec()
#     station_ra, stations_dec, _ = pos.at(t).observe(earth + tianhe).apparent().radec()
#     return abs(moon_ra.hours - station_ra.hours) + abs(moon_dec.degrees - stations_dec.degrees)print('总过境次数', len(pass_pairs))
calc_transit_moon.step_days = 1 / DAY_S
for p in pass_pairs:tm0 = ts.tt_jd(p[0])tm1 = ts.tt_jd(p[1])times, angle = find_minima(tm0, tm1, calc_transit_moon)for t, ang in zip(times, angle):# 挑选夹角小于5的显示出来，角度太大的肯定不会凌月if ang < 1:t_str = str(t.astimezone(beijing))[:19]alz, az, dist = difference.at(t).altaz()unit = '°'print('{} \t角距:{:>5.3f}{} \t高度:{:>5.1f}° \t方位角:{:>6.1f}° \t距离:{:>7.1f}km'.format(t_str, ang, unit, alz.degrees, az.degrees, dist.km))print('elapsed time:', time.perf_counter() - start_time)

简单说明一下，代码里查询范围是当日至以后300天，这个明显是不可能的，空间站的轨道信息每天都会更新，所以关于空间站的计算，即使没有大的变轨操作，计算结果在1周后误差也会越来越大，这里只当是做个数学游戏。计算后在stellarium里看看是不是符合。另外在百度地图里查到的坐标，需要转换为WGS84坐标。在stellarium里设定自己坐标时，填这个转换后的坐标。
下面是运行结果：
在这里插入图片描述

挑一个数值小的，2021-12-11 22:09:14 角距离:0.400°
在这里插入图片描述

蓝色框里就是天和空间站，放大后
在这里插入图片描述

再看2021-08-18 22:49:18 角距离:0.348°
在这里插入图片描述

这次反而没有凌月，只是接近月亮
还有一个最小的：2021-12-16 15:22:44 角距离:0.072° 高度: 7.4°
在这里插入图片描述

角度这么小，照理说肯定凌月了，但模拟显示没有，不知道为什么。

另外，没有考虑月相和凌月时的高度等其它影响观测的因素，比如2021-10-04 03:34:51 角距离:0.089° 高度: 6.4° （7月10日更新天和轨道数据后计算结果，和昨天不完全一样），模拟显示：
在这里插入图片描述
这一天是其实是残月，并不适合观测凌月

之所以写出来，只是想和有共同爱好的网友交流一下，希望有高手指定一下。代码可以初步计算出凌日月时间，但没办法做到向上面提到的网站那样，可以计算出在观察点附件一定范围内的最佳观测点。

下面是https://transit-finder.com显示的最佳观察点地图，这个是国际空间站凌日的，中心红线是最佳点，红色范围都可以观测到。
在这里插入图片描述
目前正在尝试实现这个功能。希望有兴趣的一起交流。

这篇关于天和空间站凌月凌日时间计算的尝试的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！


http://www.chinasem.cn/article/190937。

相关文章

服务器集群同步时间手记

服务器集群同步时间手记

1.时间服务器配置（必须root用户） （1）检查ntp是否安装 [root@node1 桌面]# rpm -qa|grep ntpntp-4.2.6p5-10.el6.centos.x86_64fontpackages-filesystem-1.41-1.1.el6.noarchntpdate-4.2.6p5-10.el6.centos.x86_64 （2）修改ntp配置文件 [r
阅读更多...
poj 1113 凸包+简单几何计算

poj 1113 凸包+简单几何计算

题意： 给N个平面上的点，现在要在离点外L米处建城墙，使得城墙把所有点都包含进去且城墙的长度最短。 解析： 韬哥出的某次训练赛上A出的第一道计算几何，算是大水题吧。 用convexhull算法把凸包求出来，然后加加减减就A了。 计算见下图： 好久没玩画图了啊好开心。 代码： #include <iostream>#include <cstdio>#inclu
阅读更多...
uva 1342 欧拉定理（计算几何模板）

uva 1342 欧拉定理（计算几何模板）

题意： 给几个点，把这几个点用直线连起来，求这些直线把平面分成了几个。 解析： 欧拉定理： 顶点数 + 面数 - 边数= 2。 代码： #include <iostream>#include <cstdio>#include <cstdlib>#include <algorithm>#include <cstring>#include <cmath>#inc
阅读更多...
uva 11178 计算集合模板题

uva 11178 计算集合模板题

题意： 求三角形行三个角三等分点射线交出的内三角形坐标。 代码： #include <iostream>#include <cstdio>#include <cstdlib>#include <algorithm>#include <cstring>#include <cmath>#include <stack>#include <vector>#include <
阅读更多...
XTU 1237 计算几何

XTU 1237 计算几何

题面： Magic Triangle Problem Description: Huangriq is a respectful acmer in ACM team of XTU because he brought the best place in regional contest in history of XTU. Huangriq works in a big compa
阅读更多...
MiniGPT-3D, 首个高效的3D点云大语言模型，仅需一张RTX3090显卡，训练一天时间，已开源

MiniGPT-3D, 首个高效的3D点云大语言模型，仅需一张RTX3090显卡，训练一天时间，已开源

项目主页：https://tangyuan96.github.io/minigpt_3d_project_page/ 代码：https://github.com/TangYuan96/MiniGPT-3D 论文：https://arxiv.org/pdf/2405.01413 MiniGPT-3D在多个任务上取得了SoTA，被ACM MM2024接收，只拥有47.8M的可训练参数，在一张RTX
阅读更多...
音视频入门基础：WAV专题（10）——FFmpeg源码中计算WAV音频文件每个packet的pts、dts的实现

音视频入门基础：WAV专题（10）——FFmpeg源码中计算WAV音频文件每个packet的pts、dts的实现

一、引言 从文章《音视频入门基础：WAV专题（6）——通过FFprobe显示WAV音频文件每个数据包的信息》中我们可以知道，通过FFprobe命令可以打印WAV音频文件每个packet（也称为数据包或多媒体包）的信息，这些信息包含该packet的pts、dts： 打印出来的“pts”实际是AVPacket结构体中的成员变量pts，是以AVStream->time_base为单位的显
阅读更多...
批处理以当前时间为文件名创建文件

批处理以当前时间为文件名创建文件

批处理以当前时间为文件名创建文件 批处理创建空文件 有时候，需要创建以当前时间命名的文件，手动输入当然可以，但是有更省心的方法吗？ 假设我是 windows 操作系统，打开命令行。 输入以下命令试试： echo %date:~0,4%_%date:~5,2%_%date:~8,2%_%time:~0,2%_%time:~3,2%_%time:~6,2% 输出类似： 2019_06
阅读更多...
【MRI基础】TR 和 TE 时间概念

【MRI基础】TR 和 TE 时间概念

重复时间 (TR) 磁共振成像 (MRI) 中的 TR（重复时间，repetition time）是施加于同一切片的连续脉冲序列之间的时间间隔。具体而言，TR 是施加一个 RF（射频）脉冲与施加下一个 RF 脉冲之间的持续时间。TR 以毫秒 (ms) 为单位，主要控制后续脉冲之前的纵向弛豫程度（T1 弛豫），使其成为显著影响 MRI 中的图像对比度和信号特性的重要参数。 回声时间 (TE)
阅读更多...
计算数组的斜率，偏移，R2

计算数组的斜率，偏移，R2

模拟Excel中的R2的计算。         public bool fnCheckRear_R2(List<double[]> lRear, int iMinRear, int iMaxRear, ref double dR2)         {             bool bResult = true;             int n = 0;             dou
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2