Ephemeris difference with JPL Horizons

[Errabundo66]

Hello.

I'm not sure if it's my error, or the result is OK.
With the code I attach at the botton, my results for the ra/dec position, for example, for asteroid (1234) Elyna on 2025/06/05 00:00 are:

La antigüedad del archivo MPCORB.DAT es: Wed May 24 20:36:18 2023
La antigüedad en días es: 16.89586230268633
1283767 minor planets loaded
designation_packed                            01234
magnitude_H                                   10.83
magnitude_G                                    0.15
epoch_packed                                  K232P
mean_anomaly_degrees                      162.15207
argument_of_perihelion_degrees             90.41267
longitude_of_ascending_node_degrees       304.33904
inclination_degrees                         8.52628
eccentricity                               0.086996
mean_daily_motion_degrees                  0.188264
semimajor_axis_au                          3.015026
uncertainty                                       0
reference                                 E2023-K85
observations                                 4248.0
oppositions                                      40
observation_period                        1931-2023
rms_residual_arcseconds                        0.51
coarse_perturbers                               M-v
precise_perturbers                              3Ek
computer_name                              MPCLINUX
hex_flags                                      0000
designation                            (1234) Elyna
last_observation_date                      20230521
Name: (1234) Elyna, dtype: object
Skyfield (1234) Elyna 
        Absolute Magnitude (H from data file): 10.83 
        Earth Body Distance AU: 3.032466379564208 
        Apparent Magnitude (calculated): 16.48225268603929 

ra: 00h 11m 44.86s ******* THIS ONE
dec: +07deg 02' 33.2"******* THIS ONE
distancia desde Observatory: 3.03247 au
phase angle: 19deg 32' 06.1"
PAB-LON 5.5668
PAB-LAT 356.0476
['2019-07-23T05:18:34Z', '2020-03-05T01:08:28Z', '2020-10-31T08:28:14Z', '2021-06-20T18:34:31Z', '2022-02-03T01:42:22Z', '2022-09-16T05:59:28Z', '2023-04-18T14:35:03Z', '2023-11-26T15:53:01Z', '2024-07-02T11:01:24Z']
[1 0 1 0 1 0 1 0 1]
2019-07-23T05:18:34Z
['2019-07-23T05:18:34Z', '2020-10-31T08:28:14Z', '2022-02-03T01:42:22Z', '2023-04-18T14:35:03Z', '2024-07-02T11:01:24Z']

But with JPL Horizons I got:

*******************************************************************************
JPL/HORIZONS                1234 Elyna (1931 UF)           2023-Jun-10 08:57:23
Rec #:    1234 (+COV) Soln.date: 2023-Apr-24_20:42:20   # obs: 4189 (1931-2023)
 
IAU76/J2000 helio. ecliptic osc. elements (au, days, deg., period=Julian yrs):
 
  EPOCH=  2457528.5 ! 2016-May-20.00 (TDB)         Residual RMS= .24411
   EC= .08789671871351162  QR= 2.748091614469111   TP= 2457209.7198939477
   OM= 304.7147484421288   W=  86.86544072817924   IN= 8.518395122054796
   A= 3.012917145296339    MA= 60.07786649368683   ADIST= 3.277742676123568
   PER= 5.22985            N= .188461781           ANGMOM= .029743418
   DAN= 2.97534            DDN= 3.00408            L= 31.545295
   B= 8.505556             MOID= 1.77134001        TP= 2015-Jul-06.2198939477
 
Asteroid physical parameters (km, seconds, rotational period in hours):
   GM= n.a.                RAD= 11.438             ROTPER= 5.4221
   H= 10.84                G= .150                 B-V= n.a.
                           ALBEDO= .162            STYP= K
 
ASTEROID comments: 
1: soln ref.= JPL#56, OCC=0
2: source=ORB
*******************************************************************************
 
 
*******************************************************************************
Ephemeris / WWW_USER Sat Jun 10 08:57:23 2023 Pasadena, USA      / Horizons    
*******************************************************************************
Target body name: 1234 Elyna (1931 UF)            {source: JPL#56}
Center body name: Earth (399)                     {source: DE441}
Center-site name: Observatorio Polop
*******************************************************************************
Start time      : A.D. 2025-Jun-05 00:00:00.0000 UT      
Stop  time      : A.D. 2025-Jun-05 12:00:00.0000 UT      
Step-size       : 60 minutes
*******************************************************************************
Target pole/equ : undefined
Target radii    : 11.438 km                                                    
Center geodetic : 359.8559, 38.6172241, .227616   {E-lon(deg),Lat(deg),Alt(km)}
Center cylindric: 359.8559, 4990.13617, 3959.347  {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                                          
Small-body perts: Yes                             {source: SB441-N16}
Atmos refraction: NO (AIRLESS)
RA format       : HMS
Time format     : CAL 
Calendar mode   : Mixed Julian/Gregorian
EOP file        : eop.230609.p230901                                           
EOP coverage    : DATA-BASED 1962-JAN-20 TO 2023-JUN-09. PREDICTS-> 2023-AUG-31
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 )                           
Table format    : Comma Separated Values (spreadsheet)
*******************************************************************************
Initial IAU76/J2000 heliocentric ecliptic osculating elements (au, days, deg.):
  EPOCH=  2457528.5 ! 2016-May-20.00 (TDB)         Residual RMS= .24411        
   EC= .08789671871351162  QR= 2.748091614469111   TP= 2457209.7198939477      
   OM= 304.7147484421288   W=  86.86544072817924   IN= 8.518395122054796       
  Equivalent ICRF heliocentric cartesian coordinates (au, au/d):
   X=-5.565583213271638E-01  Y= 2.536735593917346E+00  Z= 1.289177841465340E+00
  VX=-1.014677025136792E-02 VY=-4.497062132468969E-04 VZ=-1.656187129157970E-03
Asteroid physical parameters (km, seconds, rotational period in hours):        
   GM= n.a.                RAD= 11.438             ROTPER= 5.4221              
   H= 10.84                G= .150                 B-V= n.a.                   
                           ALBEDO= .162            STYP= K                     
*******************************************************************************
 Date__(UT)__HR:MN, , , R.A._____(ICRF), DEC______(ICRF),
*********************************************************
$$SOE
 **2025-Jun-05 00:00, ,m, 00 11 36.871867, +07 01 46.00405,**
 2025-Jun-05 01:00, ,m, 00 11 39.562063, +07 02 10.52364,

With Skyfiled:

ra: 00h 11m 44.86s
dec: +07deg 02' 33.2"

With JPL Horiz.

ra: 00 11 36.871867
dec: +07 01 46.00405,**

RA is almost OK, but in dec I get about 1 minute difference.

I'm not sure if the problem is the gravitational planets perturbations or not.
Does JPL gives the results for a planets perturbed orbit ? I thought it was only osculating elements calculations...

My actual code:

from skyfield.api import load, wgs84
from skyfield.data import mpc
from skyfield.almanac import find_discrete, oppositions_conjunctions
import pandas as pd
from skyfield.constants import GM_SUN_Pitjeva_2005_km3_s2 as GM_SUN
import math
import numpy as np
from skyfield.framelib import ecliptic_frame
from skyfield.positionlib import build_position
import os, time
from datetime import datetime


def getApparentMagnitude_HG( H_absoluteMagnitude, G_slope, bodyEarthDistanceAU, bodySunDistanceAU, earthSunDistanceAU ):
    beta = math.acos( \
                        ( bodySunDistanceAU * bodySunDistanceAU + \
                          bodyEarthDistanceAU * bodyEarthDistanceAU - \
                          earthSunDistanceAU * earthSunDistanceAU ) / \
                        ( 2 * bodySunDistanceAU * bodyEarthDistanceAU ) \
                    )

    psi_t = math.exp( math.log( math.tan( beta / 2.0 ) ) * 0.63 )
    Psi_1 = math.exp( -3.33 * psi_t )
    psi_t = math.exp( math.log( math.tan( beta / 2.0 ) ) * 1.22 )
    Psi_2 = math.exp( -1.87 * psi_t )

    # Have found a combination of G_slope, Psi_1 and Psi_2 can lead to a negative value in the log calculation.
    try:
        apparentMagnitude = H_absoluteMagnitude + \
                            5.0 * math.log10( bodySunDistanceAU * bodyEarthDistanceAU ) - \
                            2.5 * math.log10( ( 1 - G_slope ) * Psi_1 + G_slope * Psi_2 )

    except:
        apparentMagnitude = None

    return apparentMagnitude


filename = "MPCORB.DAT"
creation_time = os.path.getctime(filename)

# Obtener la fecha actual
current_time = time.time()

# Calcular la antigüedad en días
age_days = (current_time - creation_time) / (24 * 3600)

# Convertir el tiempo en formato legible para humanos
formatted_creation_time = time.ctime(creation_time)

print("La antigüedad del archivo", filename, "es:", formatted_creation_time)
print("La antigüedad en días es:", age_days)

# Cargar los datos de los planetas y la Tierra
#load = Loader('~/skyfield-data')
eph = load('de440s.bsp')

# Definir la ubicación del observador en la Tierra (por ejemplo, Benidorm)
observatory = wgs84.latlon(38, -0.2)

ts = load.timescale()
t = ts.utc(2025, 6, 5)

sun, earth = eph['sun'], eph['earth']

name_body = str('(1234) Elyna')

# #Cargar el archivo MPCORB.dat
mpcorb = load.open('MPCORB.DAT')

for i in range(43):
    repr(mpcorb.readline())

with mpcorb as f:
    minor_planets = mpc.load_mpcorb_dataframe(f)

print(minor_planets.iloc[0])

Save Dataframe in a binary archive for later use, quicker that reading MPCORB file -- Guardar el DataFrame en un archivo binario
minor_planets.to_pickle('data.pkl')

# Cargar el DataFrame desde el archivo binario
#minor_planets = pd.read_pickle('data.pkl')

print(minor_planets.shape[0], 'minor planets loaded')

# Index by designation for fast lookup.
minor_planets = minor_planets.set_index('designation', drop=False)

# Sample lookups.
row = minor_planets.loc[name_body]
print (row)

#Efemérides del asteroide
asteroid = sun + mpc.mpcorb_orbit(row, ts, GM_SUN)

#Efemerides en tiempo t
earthSunDistance = ( earth + observatory ).at( t ).observe( sun ).distance()
sunBodyDistance = ( sun ).at( t ).observe( asteroid ).distance()

ra_ast, dec_ast, earthBodyDistance = ( earth + observatory ).at( t ).observe( asteroid ).radec()

#Calculate apparent magnitude
apparentMagnitude = getApparentMagnitude_HG( row[ "magnitude_H" ], row[ "magnitude_G" ], earthBodyDistance.au, sunBodyDistance.au, earthSunDistance.au )

print( "Skyfield", name_body,
        "\n\tAbsolute Magnitude (H from data file):", row[ "magnitude_H" ],
        "\n\tEarth Body Distance AU:", earthBodyDistance.au,
        "\n\tApparent Magnitude (calculated):", apparentMagnitude,
        "\n" )

print('ra:',ra_ast)
print('dec:',dec_ast)
print ('distancia desde Observatory:',earthBodyDistance)

#Calculate PABL etc
#CReate asteroids vectors to earth and sun

ast_to_sun_vector = asteroid.at(t).observe(sun)
ast_to_earth_vector = asteroid.at(t).observe(earth)

#Calculate Phase Angle (Sun-Target-Earth angle)
phase_angle = ast_to_sun_vector.separation_from(ast_to_earth_vector)

print('phase angle:', phase_angle)


#Normalize those vectors for adding and creating bisector vector (sum_vector)
x1, y1, z1 = ast_to_sun_vector.xyz.au
x2, y2, z2 = ast_to_earth_vector.xyz.au

vector1 = np.array([x1, y1, z1])
vector2 = np.array([x2, y2, z2])

normalized_vector1 = vector1 / np.linalg.norm(vector1)
normalized_vector2 = vector2 / np.linalg.norm(vector2)

sum_vector = normalized_vector1 + normalized_vector2

#Create a position with (-sum_vector) centered at sun (minus, because we wnat the long in the direction to sun, not from sun)

position = build_position(center= sun.at(t), position_au= -sum_vector, t=t)

lat, lon, distance2 = position.frame_latlon(ecliptic_frame)


print('PAB-LON {:.4f}'.format(lat.degrees))
print('PAB-LAT {:.4f}'.format(lon.degrees))



#Find next opposition (for example between 2019 and 2025)-- Prueba para encontrar la proxima oposición

t0 = ts.utc(2019, 1, 1)
t1 = ts.utc(2025, 1, 1)
f_op = oppositions_conjunctions(eph, asteroid)
t, y = find_discrete(t0, t1, f_op)

print(t.utc_iso())
print(y)

# Encontrar el índice de la primera oposición con y = 1
first_opposition_index = np.where(y == 1)[0][0]

# Imprimir la fecha y hora de la primera oposición
print(t[first_opposition_index].utc_iso())

# Imprimir todas las fechas y horas de las oposiciones con y = 1
oppositions_indices = np.where(y == 1)[0]
oppositions_times = t[oppositions_indices]
print(oppositions_times.utc_iso())