1

我正在尝试计算近地点和远地点（或一般给定第二个物体，如太阳和行星等）

from skyfield import api, almanac
from scipy.signal import argrelextrema
import numpy as np

e = api.load('de430t.bsp')

def apsis(year = 2019, body='moon'):
    apogees = dict()
    perigees = dict()
    planets = e
    earth, moon = planets['earth'], planets[body]

    t = ts.utc(year, 1, range(1,367)) 
    dt = t.utc_datetime()

    astrometric = earth.at(t).observe(moon)
    _, _, distance = astrometric.radec()

    #find perigees, at day precision
    localmaxes = argrelextrema(distance.km, np.less)[0]
    for i in localmaxes:
       # get minute precision
       t2 = ts.utc(dt[i].year, dt[i].month, dt[i].day-1, 0, range(2881))
       dt2 = t2.utc_datetime()  # _and_leap_second()
       astrometric2 = earth.at(t2).observe(moon)
       _, _, distance2 = astrometric2.radec()
       m = min(distance2.km)
       daindex = list(distance2.km).index(m)
       perigees[dt2[daindex]] = m

    #find apogees, at day precision
    localmaxes = argrelextrema(distance.km, np.greater)[0]
    for i in localmaxes:
        # get minute precision
        t2 = ts.utc(dt[i].year, dt[i].month, dt[i].day-1, 0, range(2881))
        dt2 = t2.utc_datetime()  
        astrometric2 = earth.at(t2).observe(moon)
        _, _, distance2 = astrometric2.radec()
        m = max(distance2.km)
        daindex = list(distance2.km).index(m)
        apogees[dt2[daindex]] = m

    return apogees, perigee

当我为 2019 年运行此程序时，下一个远地点计算在 2019-09-13 13:16。这与John Walker (13:33)、Fred Espenak (13:32)、Time and Date dot com (13:32)等表格相差几分钟。

由于舍入与截断秒数等原因，我希望在其他来源中看到一分钟的差异，但超过 15 分钟的差异似乎不寻常。我已经用 de431t 和 de421 ephemeris 进行了尝试，结果相似。

这里有什么区别？我正在计算每个身体中心的距离，对吧？我在搞什么鬼？

4

1 回答 1

1

经过更多研究并将天空场输出与 JPL 地平线的输出进行比较后，天空场的计算似乎是正确的，至少与 JPL 星历表相比（这并不奇怪）

我将上面的代码片段切换为使用与 HORIZONS 相同的（大量）de432t SPICE 内核。这与 HORIZONS 输出（见下文，标记的各种来源报告的远地点）一致，月球开始远离（deldot 或观察者（地心地球）和目标天体（地心月球）之间的距离率变为负数

Ephemeris / WWW_USER Fri Sep 13 17:05:39 2019 Pasadena, USA      / Horizons    
*******************************************************************************
Target body name: Moon (301)                      {source: DE431mx}
Center body name: Earth (399)                     {source: DE431mx}
Center-site name: GEOCENTRIC
*******************************************************************************
Start time      : A.D. 2019-Sep-13 13:10:00.0000 UT      
Stop  time      : A.D. 2019-Sep-13 13:35:00.0000 UT      
Step-size       : 1 minutes
*******************************************************************************
Target pole/equ : IAU_MOON                        {East-longitude positive}
Target radii    : 1737.4 x 1737.4 x 1737.4 km     {Equator, meridian, pole}    
Center geodetic : 0.00000000,0.00000000,0.0000000 {E-lon(deg),Lat(deg),Alt(km)}
Center cylindric: 0.00000000,0.00000000,0.0000000 {E-lon(deg),Dxy(km),Dz(km)}
Center pole/equ : High-precision EOP model        {East-longitude positive}
Center radii    : 6378.1 x 6378.1 x 6356.8 km     {Equator, meridian, pole}    
Target primary  : Earth
Vis. interferer : MOON (R_eq= 1737.400) km        {source: DE431mx}
Rel. light bend : Sun, EARTH                      {source: DE431mx}
Rel. lght bnd GM: 1.3271E+11, 3.9860E+05 km^3/s^2                              
Atmos refraction: NO (AIRLESS)
RA format       : HMS
Time format     : CAL 
EOP file        : eop.190912.p191204                                           
EOP coverage    : DATA-BASED 1962-JAN-20 TO 2019-SEP-12. PREDICTS-> 2019-DEC-03
Units conversion: 1 au= 149597870.700 km, c= 299792.458 km/s, 1 day= 86400.0 s 
Table cut-offs 1: Elevation (-90.0deg=NO ),Airmass (>38.000=NO), Daylight (NO )
Table cut-offs 2: Solar elongation (  0.0,180.0=NO ),Local Hour Angle( 0.0=NO )
Table cut-offs 3: RA/DEC angular rate (     0.0=NO )                           
*******************************************************************************
 Date__(UT)__HR:MN                delta      deldot
***************************************************
$$SOE
 2019-Sep-13 13:10     0.00271650099697   0.0000340
 2019-Sep-13 13:11     0.00271650100952   0.0000286
 2019-Sep-13 13:12     0.00271650101990   0.0000232
 2019-Sep-13 13:13     0.00271650102812   0.0000178
 2019-Sep-13 13:14     0.00271650103417   0.0000124
 2019-Sep-13 13:15     0.00271650103805   0.0000070
 2019-Sep-13 13:16     0.00271650103977   0.0000016 <----- Skyfield, HORIZONS
 2019-Sep-13 13:17     0.00271650103932  -0.0000038
 2019-Sep-13 13:18     0.00271650103670  -0.0000092
 2019-Sep-13 13:19     0.00271650103191  -0.0000146
 2019-Sep-13 13:20     0.00271650102496  -0.0000200
 2019-Sep-13 13:21     0.00271650101585  -0.0000254
 2019-Sep-13 13:22     0.00271650100456  -0.0000308
 2019-Sep-13 13:23     0.00271650099112  -0.0000362
 2019-Sep-13 13:24     0.00271650097550  -0.0000416
 2019-Sep-13 13:25     0.00271650095772  -0.0000470
 2019-Sep-13 13:26     0.00271650093778  -0.0000524
 2019-Sep-13 13:27     0.00271650091566  -0.0000578
 2019-Sep-13 13:28     0.00271650089139  -0.0000632
 2019-Sep-13 13:29     0.00271650086494  -0.0000686
 2019-Sep-13 13:30     0.00271650083633  -0.0000740
 2019-Sep-13 13:31     0.00271650080556  -0.0000794
 2019-Sep-13 13:32     0.00271650077262  -0.0000848  <------ Espenak, T&D.com
 2019-Sep-13 13:33     0.00271650073751  -0.0000902 
 2019-Sep-13 13:34     0.00271650070024  -0.0000956
 2019-Sep-13 13:35     0.00271650066081  -0.0001010

$$EOE

再看看 Espenak 的页面，他的计算是基于 Jean Meeus 的天文算法书（任何玩过这些东西的人都必须拥有的）。那本书中的月球星历来自Jean Chapront 的ELP2000/82。虽然这已安装到 DE430（以及其他）中，

果然，在今天 2019 年 9 月 13 日使用 ELP2000 模型找到最大月球距离时。你得到 2019-09-13 13:34。请参阅下面的代码。

Meeus 的公式基于 1982 年版的 Ephemeride Lunaire Parisienne，下面的源代码利用了 Chapront 的 2002 年更新，但几乎是其他来源提出的。

所以我认为我的答案是，它们是不同的答案，因为它们使用不同的模型。Skyfield 正在利用 JPL Development 星历表示为数值积分的模型，而 ELP 是一种更具分析性的方法。

最后我意识到这是一个挑剔的选择，我只是想更好地了解我正在使用的工具。但它引出了一个问题，哪种方法更准确？

据我所知，DE430 及其同位素已适合观测数据，即月球激光测距 (LLR) 测量。如果只是出于 LLR 的考虑，我想我会坚持使用 Skyfield 来计算月球距离。

from elp_mpp02 import mpp02 as mpp
import julian
import pytz
import datetime

def main():
    mpp.dataDir = 'ELPmpp02'
    mode = 1  # Historical mode
    jd = 2451545
    data = dict()
    maxdist = 0
    apogee = None
    for x in range(10,41):
        dt = datetime.datetime(2019, 9, 13, 13, x, tzinfo=pytz.timezone("UTC"))
        jd = julian.to_jd(dt, fmt='jd')
        lon, lat, dist = mpp.compute_lbr(jd, mode)
        if dist > maxdist:
            maxdist = dist
            apogee = dt
    print(f"{maxdist:.2} {apogee}")
于 2019-09-13T21:08:29.737 回答