Better way to calculate apparent angular separation of two objects in Skyfield?

Question

UPDATE: Skyfield has just had a significant revision, including expanded documentation and a method for angular separation - see the accepted answer.

I'm calculating the apparent angular separation between two objects using Skyfield. I didn't find a method within the package, so I've "invented" a way by calculating the dot product between the two apparent position vectors.

Is this the best way to do it currently? Is it essentially correct, within Skyfield's scope?

def separation(seconds, lat, lon):

    lat, lon, seconds = float(lat), float(lon), float(seconds) # necessary it seems

    place = earth.topos(lat, lon)

    jd = JulianDate(utc=(2016, 3, 9, 0, 0, seconds))

    mpos = place.at(jd).observe(moon).apparent().position.km
    spos = place.at(jd).observe(sun).apparent().position.km

    mlen = np.sqrt((mpos**2).sum())
    slen = np.sqrt((spos**2).sum())

    sepa = ((3600.*180./np.pi) *
            np.arccos(np.dot(mpos, spos)/(mlen*slen)))

    return sepa


from skyfield.api import load, now, JulianDate
import numpy as np
from scipy.optimize import minimize

data = load('de421.bsp')

sun   = data['sun']
earth = data['earth']
moon  = data['moon']

sep = separation(12000, 32.5, 215.1)

print "sun-moon aparent separation: ", sep, " arcsec"

Designing

Here is a related Skyfield question. – Designing

I don't know Skyfield (but it looks interesting), but your code looks ok to me. Alternatively, use the cos rule of spherical trig on the RA and declination of the two apparent positions. There will be a tiny difference between the two methods due to rounding errors. Also see en.wikipedia.org/wiki/Great-circle_distance – Chemesh

Skyfield rocks! @PM2Ring thanks for the suggestion. I'm sticking with cartesian here because that's what I'm using in the application. I feel safer using cartesian coordinates and avoiding antipodals and zeros in denominators (unless the moon actually hits the earth some day) although I don't know if that's justified. – Designing

Answer 1

Skyfield now supports a method that returns the angular separation between two positions:

http://rhodesmill.org/skyfield/api-position.html#skyfield.positionlib.ICRF.separation_from

Rhiamon