5
$\begingroup$

I am calculating the apparent planetary position of planets using the cspice library (I am specifically using rust-spice, which is a wrapper of cspice in Rust).

when I calculate the same values in the Horizons app there is a discrepancy between the results.

The results in Horizons vs SPICE:

  1. RA : spice: 181.29145390230866, Horizons: 181.57546
  2. Dec: spice: 1.5006592306305064, Horizons: 1.37450
  3. Ecliptic Longitude: spice: 181.2848320891579, Horizons: 180.8987129
  4. Ecliptic Latitude: spice: 1.5006592306305062, Horizons: 1.8877102

As you can see, there is a significant difference in the values obtained, and I know I am missing something in the setup. I tried various values of aberration corrections, but no values returns me anything close to Horizons' results. How can I get the same results as Horizons'?

I am posting the results of Horizons as well as excerpts of my Rust code:

Kernels (kernels.furnsh):

\begindata
         KERNELS_TO_LOAD = (
            'JPL/de441.bsp'
            'JPL/naif0012.tls'
                           )
\begintext

My Rust code:

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

Results (for each value of 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

The Horizons results:

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

**********************************************************************************************************************
$\endgroup$
8
  • $\begingroup$ Differences look too big but are you calculating for the same time ? I can't see your definition of et in the rust code but 2023-10-05 12:25 UTC (in HORIZONS) and 2023-10-05 12:25 ET/TDB are not the same time and so won't give the same positions. Next thing to check after the time would be whether the code and HORIZONS are setup to give the results in the same reference frame $\endgroup$ Oct 9, 2023 at 15:37
  • $\begingroup$ Do you have a reference for 'cspice'? Is it part of this SPICE (Wikipedia page)? Is it the C Toolkit? Or something else? $\endgroup$ Oct 10, 2023 at 20:39
  • $\begingroup$ (SPICE is a heavily overloaded acronym. In this case, it may be for "spacecraft planet instrument C-matrix ("camera matrix") events".) $\endgroup$ Oct 10, 2023 at 20:42
  • $\begingroup$ Related: Discrepancy between HORIZONS and SPICE $\endgroup$ Oct 10, 2023 at 20:46
  • $\begingroup$ JPL Horizons: JPL Horizons On-Line Ephemeris System (not to be confused with New Horizons) $\endgroup$ Oct 10, 2023 at 20:49

1 Answer 1

7
$\begingroup$

The issue is with how you have configured HORIZONS to show the output. You have apparent RA, Dec selected as output columns which includes the precession and nutation, in addition to the light time and annual/"stellar" aberration being performed by spkpos. If you go into the Table Options for HORIZONS and select quantity 1 (Astrometric RA & DEC)as well then you should see: Output of JPL HORIZONS showing ICRF and apparent RA, Dec values

The positions under the R.A._______(ICRF)_______DEC column compares with the abcorr='LT' result from SPICE to a few millarcseconds.

$\endgroup$
1
  • $\begingroup$ I think (have not tested) that you might need the high precision Earth Orientation Parameters binary PCK file as detailed in Tutorial #23 $\endgroup$ Oct 11, 2023 at 4:39

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .