Better way to calculate apparent angular separatio

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"

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