Asteroid opposition, Phase angle bisector longitude, etc?

[Errabundo66]

Hello again.

I've been looking for info about that, with no success. I'd like to calculate some data for an asteroid. I used previously the JPL Horizons ephemeris query in Python (I attach a capture image). But I prefer to calculate all, if possible, from the MPCORB.dat. In that way there is no need to be connected to internet for the calculations.

If you see the attached image, I look for:

• Magnitude
• Phase angle
• Phase angle bisector longitude and latitude

We are dealing with asteroids light curve inversion and need that data. Are those calculations made by Skyfield? Are they calculated by the ephemeris function?

And also I have an option for calculating next Opposition date for an asteroid (see second image) Can we use ALMANAC for an asteroid? Or is it only for main planets?

Thank you for your help.

[Errabundo66]

Hello. I've managed getting the phase angle for the asteroid, using "separation_from()".

ra: 01h 03m 39.32s
dec: +13deg 25' 44.3"
distance: 3.32236 au
phase angle: 16deg 14' 26.1"

I tried also using ICRF "phase angle" function, but I got some errors (in sume, I don't still know how to proceed with vectors in Skyfield)... I got the error:

    phase_angle_method2 = ICRF.phase_angle(asteroid, sun)
AttributeError: 'VectorSum' object has no attribute 't'

I'd like to know how to deal with vectors, because I need the Phase angle bisector Galactic longitude and latitude, and I don't really know how to calculate them.

That's my actual code:

from skyfield.api import Loader, Topos
from skyfield.api import load, wgs84
from skyfield.data import mpc
from skyfield.almanac import find_discrete
from skyfield.positionlib import ICRF
import pandas as pd
from skyfield.constants import GM_SUN_Pitjeva_2005_km3_s2 as GM_SUN

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

# Definir la ubicación del observador en la Tierra (por ejemplo, Madrid)
observatory = wgs84.latlon(40.4168, -3.7038)

# 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])

# 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['(1234) Elyna']
print (row)

# Generating a position.

ts = load.timescale()
t = ts.utc(2020, 5, 31)

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

asteroid = sun + mpc.mpcorb_orbit(row, ts, GM_SUN)

sun_position = asteroid.at(t).observe(sun)
earth_position = asteroid.at(t).observe(earth)

phase_angle = sun_position.separation_from(earth_position)

ra, dec, distance = earth.at(t).observe(asteroid).radec()
print(ra)
print(dec)
print (distance)

print(phase_angle)
print(phase_angle_method2)

[Errabundo66]

Another try:

I successfully got the Phase Angle.

But as I explained I need Phase Angle Bisector Galactic longitude and latitude.

I've tried getting the Bisector Vector = using both vectors Asteroid-Earth and Asteroid-sun, converting them to Numpy vectors, normalising them, adding one to another... and converting them again to a position.

But the result is not ok.

from skyfield.api import Loader, Topos
from skyfield.api import load, wgs84
from skyfield.data import mpc
from skyfield.almanac import find_discrete
from skyfield.positionlib import ICRF
import pandas as pd
from skyfield.constants import GM_SUN_Pitjeva_2005_km3_s2 as GM_SUN
from skyfield.vectorlib import VectorFunction, VectorSum
import numpy as np
from skyfield.framelib import galactic_frame
from skyfield.framelib import ecliptic_frame
from skyfield.positionlib import build_position

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

# Definir la ubicación del observador en la Tierra (por ejemplo, Madrid)
observatory = wgs84.latlon(40.4168, -3.7038)

# 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])

# 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['(1234) Elyna']
print (row)

# Generating a position.

ts = load.timescale()
t = ts.utc(2020, 5, 31)

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

asteroid = sun + mpc.mpcorb_orbit(row, ts, GM_SUN)

sun_position = asteroid.at(t).observe(sun)
earth_position = asteroid.at(t).observe(earth)

x1, y1, z1 = sun_position.frame_xyz(galactic_frame).au
x2, y2, z2 = earth_position.frame_xyz(galactic_frame).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

print(sum_vector)

position = build_position(sum_vector, t=t)

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

print(' {:.4f} latitude'.format(lat.degrees))
print(' {:.4f} longitude'.format(lon.degrees))
print(' {:.3f} au distant'.format(distance2.au))

phase_angle = sun_position.separation_from(earth_position)
print('phase angle:', phase_angle)

ra, dec, distance = earth.at(t).observe(asteroid).radec()
print('ra:',ra)
print('dec:',dec)
print ('distance:',distance)

I got:

[ 0.52355609 -1.10953336  1.55403235]
 13.8545 latitude  *** Wrong value. The correct one is: 6.8102 (see below)
 84.3732 longitude *** Wrong value. The correct one is: 11.7143 (see below)
 1.980 au distant 
phase angle: 16deg 14' 26.1" ****It's OK!!
ra: 01h 03m 39.32s
dec: +13deg 25' 44.3"
distance: 3.32236 au

The correct ones are:

Target body name: 1234 Elyna (1931 UF)            
 Date__(UT)__HR:MN      R.A._____(ICRF)  DEC______(ICRF)        phi   PAB-LON   PAB-LAT
 2020-May-31 00:00   m  01 03 49.488340  +13 26 53.51011    16.2468   11.7143    6.8102

What am I doing wrong? Thank you

[Errabundo66]

I closed it and moved to "Discussions" section! Sorry