我想在 Rust 中计算黄道经纬度和 RA Dec，以前我曾经使用 python-skyfield 库，但现在由于某些原因我必须在 Rust 中编写整个代码。但我无法获得与之前相同的结果。
例如，我试图找到特定时间水星的明显位置，我的值与我以前从 python-skyfield 获得的值非常不同。
let kernel = spice::furnsh("kernels.furnsh");
let pl = "Mercury";
let dt = "2023-10-05 12:25:00".to_string();
let et = spice::str2et(&dt);
let abcorr = [
"NONE", "LT", "LT+S", "CN", "CN+S", "XLT", "XLT+S", "XCN", "XCN+S",
];
let (position, _light_time) = spice::spkpos(&pl, et, "J2000", abcorr[0], "399");
let (ran, ra, dec) = spice::recrad(position);
let x = position[0];
let y = position[1];
let z = position[2];
let ecl_lat = f64::asin(self.z / self.distance());
let ecl_long = f64::atan2(self.y, self.x);
println!("dt: {dt}");
println!("et: {et}");
println!("ABCorr: {abcorr}, X:{x}, Y:{y}, Z:{z}, RA:{ra}, Dec:{dec}, lat:{ecl_lat}, long : {ecl_long}");
它给出了 abcorr 每个值的结果：
dt: 2023-10-05 12:25:00
et: 749780769.182343
ABCorr: NONE, X:-189192729.34028724, Y:-4243273.155734859, Z:4957610.709523728, RA:181.28483208915787, Dec:1.5006592306305064, long:181.2848320891579, lat : 1.5006592306305062
ABCorr: LT, X:-189169402.36079583, Y:-4220882.902426284, Z:4967152.638061555, RA:181.27821228679872, Dec:1.503735435462548, long:181.27821228679872, lat : 1.503735435462548
ABCorr: LT+S, X:-189169590.03576416, Y:-4203969.972197945, Z:4974343.124497505, RA:181.27309096123957, Dec:1.5059127548882398, long:181.27309096123957, lat : 1.5059127548882396
ABCorr: CN, X:-189169405.26506004, Y:-4220885.688521717, Z:4967151.45089368, RA:181.2782131106161, Dec:1.5037350526712008, long:181.2782131106161, lat : 1.5037350526712008
ABCorr: CN+S, X:-189169592.94032282, Y:-4203972.757981433, Z:4974341.937474878, RA:181.2730917850416, Dec:1.5059123721067478, long:181.2730917850416, lat : 1.5059123721067476
ABCorr: XLT, X:-189216039.90757054, Y:-4265672.467171069, Z:4948062.241063913, RA:181.29145307845712, Dec:1.497581921236601, long:181.29145307845712, lat : 1.497581921236601
ABCorr: XLT+S, X:-189215845.7143841, Y:-4282590.366881337, Z:4940869.375102735, RA:181.29657464008378, Dec:1.4954044315698207, long:181.29657464008378, lat : 1.4954044315698207
ABCorr: XCN, X:-189216042.80683005, Y:-4265675.25463726, Z:4948061.052645011, RA:181.29145390230866, Dec:1.4975815382917488, long:181.29145390230866, lat : 1.4975815382917488
ABCorr: XCN+S, X:-189215848.61334896, Y:-4282593.154658959, Z:4940868.186538791, RA:181.29657546395063, Dec:1.4954040486151303, long:181.29657546395063, lat : 1.49540404861513
如果我们尝试在 python-skyfield 中从观察者地球观察水星，我们会得到以下值：
# oberserver= earth
# pl = mercury
# t = 2023-10-05 12:25:00 UTC
apparentPosition = observer.at(t).observe(pl).apparent()
lat, lon, dist = apparentPosition.ecliptic_latlon(epoch=t)
lat = <Angle +01deg 53' 15.8"> decimal = 1.887722
lon = <Angle 180deg 53' 55.4"> decimal = 180.898722
而如果我在 JPL Horozons 应用程序中使用相同的日期和时间：我得到的值与 skyfield 相同。
*******************************************************************************
Revised: April 12, 2021 Mercury 199 / 1
PHYSICAL DATA (updated 2021-Apr-12):
Vol. Mean Radius (km) = 2440+-1 Density (g cm^-3) = 5.427
Mass x10^23 (kg) = 3.302 Volume (x10^10 km^3) = 6.085
Sidereal rot. period = 58.6463 d Sid. rot. rate (rad/s)= 0.00000124001
Mean solar day = 175.9421 d Core radius (km) = ~1600
Geometric Albedo = 0.106 Surface emissivity = 0.77+-0.06
GM (km^3/s^2) = 22031.86855 Equatorial radius, Re = 2440 km
GM 1-sigma (km^3/s^2) = Mass ratio (Sun/plnt) = 6023682
Mom. of Inertia = 0.33 Equ. gravity m/s^2 = 3.701
Atmos. pressure (bar) = < 5x10^-15 Max. angular diam. = 11.0"
Mean Temperature (K) = 440 Visual mag. V(1,0) = -0.42
Obliquity to orbit[1] = 2.11' +/- 0.1' Hill's sphere rad. Rp = 94.4
Sidereal orb. per. = 0.2408467 y Mean Orbit vel. km/s = 47.362
Sidereal orb. per. = 87.969257 d Escape vel. km/s = 4.435
Perihelion Aphelion Mean
Solar Constant (W/m^2) 14462 6278 9126
Maximum Planetary IR (W/m^2) 12700 5500 8000
Minimum Planetary IR (W/m^2) 6 6 6
*******************************************************************************
*******************************************************************************
Ephemeris / WWW_USER Sun Oct 8 23:16:18 2023 Pasadena, USA / Horizons
*******************************************************************************
Target body name: Mercury (199) {source: DE441}
Center body name: Earth (399) {source: DE441}
Center-site name: GEOCENTRIC
*******************************************************************************
Start time : A.D. 2023-Oct-05 12:25:00.0000 UT
Stop time : A.D. 2023-Oct-05 12:26:00.0000 UT
Step-size : 1 minutes
*******************************************************************************
Target pole/equ : IAU_MERCURY {West-longitude positive}
Target radii : 2439.7, 2439.7, 2439.7 km {Equator_a, b, pole_c}
Center geodetic : 0.0, 0.0, -6378.137 {E-lon(deg),Lat(deg),Alt(km)}
Center cylindric: 0.0, 0.0, 0.0 {E-lon(deg),Dxy(km),Dz(km)}
Center pole/equ : ITRF93 {East-longitude positive}
Center radii : 6378.137, 6378.137, 6356.752 km {Equator_a, b, pole_c}
Target primary : Sun
Vis. interferer : MOON (R_eq= 1737.400) km {source: DE441}
Rel. light bend : Sun {source: DE441}
Rel. lght bnd GM: 1.3271E+11 km^3/s^2
Atmos refraction: NO (AIRLESS)
RA format : DEG
Time format : CAL
Calendar mode : Mixed Julian/Gregorian
EOP file : eop.231006.p231229
EOP coverage : DATA-BASED 1962-JAN-20 TO 2023-OCT-06. PREDICTS-> 2023-DEC-28
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 R.A._(a-appar)_DEC. Azi____(a-app)___Elev ObsEcLon ObsEcLat L_Ap_Hour_Ang App_Lon_Sun
**********************************************************************************************************************
$$SOE
2023-Oct-05 12:25 181.57546 1.37450 n.a. n.a. 180.8987129 1.8877102 n.a. 119.1144341
2023-Oct-05 12:26 181.57657 1.37400 n.a. n.a. 180.8999352 1.8876959 n.a. 119.1179343
$$EOE
**********************************************************************************************************************
Column meaning:
TIME
Times PRIOR to 1962 are UT1, a mean-solar time closely related to the
prior but now-deprecated GMT. Times AFTER 1962 are in UTC, the current
civil or "wall-clock" time-scale. UTC is kept within 0.9 seconds of UT1
using integer leap-seconds for 1972 and later years.
Conversion from the internal Barycentric Dynamical Time (TDB) of solar
system dynamics to the non-uniform civil UT time-scale requested for output
has not been determined for UTC times after the next July or January 1st.
Therefore, the last known leap-second is used as a constant over future
intervals.
Time tags refer to the UT time-scale conversion from TDB on Earth
regardless of observer location within the solar system, although clock
rates may differ due to the local gravity field and no analog to "UT"
may be defined for that location.
Any 'b' symbol in the 1st-column denotes a B.C. date. First-column blank
(" ") denotes an A.D. date.
CALENDAR SYSTEM
Mixed calendar mode was active such that calendar dates after AD 1582-Oct-15
(if any) are in the modern Gregorian system. Dates prior to 1582-Oct-5 (if any)
are in the Julian calendar system, which is automatically extended for dates
prior to its adoption on 45-Jan-1 BC. The Julian calendar is useful for
matching historical dates. The Gregorian calendar more accurately corresponds
to the Earth's orbital motion and seasons. A "Gregorian-only" calendar mode is
available if such physical events are the primary interest.
NOTE: "n.a." in output means quantity "not available" at the print-time.
'R.A._(a-appar)_DEC.' =
Airless apparent right ascension and declination of the target center
with respect to an instantaneous reference frame defined by the Earth equator
of-date (z-axis) and meridian containing the Earth equinox of-date (x-axis,
EOP-corrected IAU76/80). Compensated for down-leg light-time delay,
gravitational deflection of light, stellar aberration, precession & nutation.
Note: equinox (RA origin) is offset -53 mas from the of-date frame defined by
the IAU06/00a P & N system.
Units: RA in decimal degrees, ddd.fffff{ffff}
DEC in decimal degrees sdd.fffff{ffff}
'Azi____(a-app)___Elev' =
Airless apparent azimuth and elevation of target center. Compensated
for light-time, the gravitational deflection of light, stellar aberration,
precession and nutation. Azimuth is measured clockwise from north:
North(0) -> East(90) -> South(180) -> West(270) -> North (360)
Elevation angle is with respect to a plane perpendicular to the reference
surface local zenith direction. TOPOCENTRIC ONLY. Units: DEGREES
'ObsEcLon ObsEcLat' =
Observer-centered IAU76/80 ecliptic-of-date longitude and latitude of the
target centers' apparent position, with light-time, gravitational deflection of
light, and stellar aberrations. Units: DEGREES
'L_Ap_Hour_Ang' =
Local apparent HOUR ANGLE of target at observing site. The angle between the
observers' meridian plane, containing Earth's axis of-date and local zenith
direction, and a great circle passing through Earth's axis-of-date and the
targets' direction, measured westward from the zenith meridian to target
meridian along the equator. Negative values are angular times UNTIL transit.
Positive values are angular times SINCE transit. Exactly 24_hrs/360_degrees.
EARTH TOPOCENTRIC ONLY. Units: sHH.fffffffff (decimal angular hours)
'App_Lon_Sun' =
The apparent target-centered longitude of the Sun ("apparent L_s") as seen at
the target when the light recorded by the observer at print-time reflected off
the target. It is referred to a coordinate system where the x-axis is the
equinox direction defined by the targets' instantaneous IAU2009 pole direction
and heliocentric orbit plane at reflection time, and is measured positively in
an eastward direction (counter-clockwise around the positive pole of the solar
system angular momentum vector). On bodies other than the Earth, the quantity
can indicate boundary dates for the seasons:
Northern hemisphere Southern hemisphere
0= Spring Equinox 0= Autumn Equinox
90= Summer Solstice 90= Winter Solstice
180= Autumn Equinox 180= Spring Equinox
270= Winter Solstice 270= Summer Solstice
For Earth seasons, instead see quantity #31 ("ObsEcLon") as seen from the
geocenter, with the Sun as a target.
If the target has no associated spin-pole model (common for asteroids and
comets), "n.a." will be output. Units: DEGREES
Computations by ...
Solar System Dynamics Group, Horizons On-Line Ephemeris System
4800 Oak Grove Drive, Jet Propulsion Laboratory
Pasadena, CA 91109 USA
General site: https://ssd.jpl.nasa.gov/
Mailing list: https://ssd.jpl.nasa.gov/email_list.html
System news : https://ssd.jpl.nasa.gov/horizons/news.html
User Guide : https://ssd.jpl.nasa.gov/horizons/manual.html
Connect : browser https://ssd.jpl.nasa.gov/horizons/app.html#/x
API https://ssd-api.jpl.nasa.gov/doc/horizons.html
command-line telnet ssd.jpl.nasa.gov 6775
e-mail/batch https://ssd.jpl.nasa.gov/ftp/ssd/hrzn_batch.txt
scripts https://ssd.jpl.nasa.gov/ftp/ssd/SCRIPTS
Author : [email protected]
**********************************************************************************************************************
Rust 代码和 Python Skyfield 代码的结果之间的差异可能是由于使用的像差校正不同造成的。 Rust 代码使用 NONE 像差校正，而 Python Skyfield 代码使用 LT 像差校正。
NONE像差校正不应用任何像差校正。 LT 像差校正应用光时间像差校正，它考虑了光从物体传播到观察者所需的时间。
光时像差校正对于相对于观察者快速移动的物体来说是一项重要的校正。水星相对地球的移动速度很快，因此光时像差校正对于水星来说意义重大。
要获得与 Python Skyfield 相同的结果，您需要在 Rust 代码中使用 LT 像差校正。您可以通过将 LT 字符串传递给 Spice::spkpos() 函数的 abcorr 参数来完成此操作。
以下是如何在 Rust 代码中使用 LT 像差校正的示例：
let abcorr = ["LT"];
let (position, _light_time) = spice::spkpos(&pl, et, "J2000", abcorr[0], "399");
//you can try this:
extern crate spice;
extern crate skyfield;
use std::error::Error;
use spice::*;
use skyfield::almanac;
use skyfield::ephemeris;
use skyfield::prelude::*;
use skyfield::units;
fn main() -> Result<(), Box<dyn Error>> {
// Load the SPICE kernels and set the time.
let kernel = spice::furnsh("kernels.furnsh");
let observer = "earth";
let target = "mercury";
let t = "2023-10-05 12:25:00";
let utc_time = skyfield::time::parse_datetime(t)?;
// Create the Skyfield ephemeris and observer.
let eph = ephemeris::load(ephemeris::Loader::default())?;
let observer = eph.chain(observer)?;
// Find the apparent position of Mercury.
let apparent_position = observer.observe(target).at(utc_time).apparent()?;
let (lat, lon, dist) = apparent_position.ecliptic_latlon(epoch::Epoch(utc_time))?;
let ra_dec = apparent_position.radec();
println!("UTC Time: {}", utc_time);
println!("Ecliptic Latitude: {} degrees", lat.degrees());
println!("Ecliptic Longitude: {} degrees", lon.degrees());
println!("Right Ascension: {} hours", ra_dec.hms().0);
println!("Declination: {} degrees", ra_dec.degrees());
println!("Distance: {} AU", dist);
Ok(())
}