如何计算视黄道Lat long和RA Dec in rust

uplii1fm  于 2023-11-19  发布在  其他
关注(0)|答案(2)|浏览(153)

我想在rust中计算黄道纬度和赤经，以前我使用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应用程序中使用相同的日期和时间：我得到的值与天空场相同。

*******************************************************************************
 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

2条答案

按热度按时间
ie3xauqp

ie3xauqp1#

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");

字符串

赞(0）分享  回复(0）举报  2023-11-19
2skhul33

2skhul332#

//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(())
}

字符串

赞(0）分享  回复(0）举报  2023-11-19
我来回答

相关问题

    查看更多

    热门标签

    更多
    JavaquerypythonNode开发语言requestUtil数据库Table后端算法LoggerMessageElementParser

    最新问答

    更多
    • 蜀ICP备2022004421号-2 NULL123 https://www.null123.com © All rights reserved
    • 本站内容来源互联网,如果侵犯您的权益请联系我们删除, 联系方式：448109455@qq.com