Phase angle bisector Ecliptic Longitude and Latitude

[Errabundo66]

Sorry to everybody!
I didn't realised there is a "Discussions" sections.
I posted that in the "Issues" section. I'm going to close that one.

I successfully got the Phase Angle.

But as I explained I need Phase Angle Bisector Ecliptic 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 answer myself:

I managed to get what I wanted: Phase angle, Phase angel bisector (PAB) ecliptic Latitude and longitude.
I didn't understood quite well how to deal with vectors in Skyfield, but that code works!!

Perhaps there is another more direct way instead using Numpy for working with vectors, but it's not terrible.

I let you the final routine for that part of my needs:

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

import numpy as np
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(2025, 4, 2)

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

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

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

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

[brandon-rhodes]

Thank you for sharing this code, @Errabundo66! I had not heard of a ‘phase angle bisector’ before, but I see that it's used in papers involving asteroid brightness, and the math you do does indeed produce a vector that bisects two others.

I can't think offhand of a quicker way to get your result with Skyfield, but if I ever run across one I'll try to remember to come back here and make a comment. In the meantime, thanks for sharing this and hopefully it will be helpful to other folks if they search for PAB and Skyfield!