Comparison of proposed barycentric corrections in astropy to Jason Eastman's web applets for barycentric velocity correction

test_bcorr.py

from astropy import units as u
from astropy import coordinates as coord
from astropy.time import Time

from astropy.coordinates.solar_system import get_body_barycentric_posvel
from astropy.coordinates.matrix_utilities import matrix_product
from astropy.coordinates.representation import CartesianRepresentation ,UnitSphericalRepresentation
from astropy.coordinates.builtin_frames import GCRS
from astropy.coordinates import SkyCoord, solar_system, EarthLocation, ICRS
KPS = u.km/u.s

from barycorr import bvc

def helio_vector(t, loc):
    """
    Compute the heliocentric velocity correction at a given time and place.
    Paramters
    ---------
    t : astropy.time.Time
        Time of the observation. Can be a Time array.
    loc : astropy.coordinates.EarthLocation
        The observer location at which to compute the correction.
    """
    vsun = get_body_barycentric_posvel('sun', t)[1]
    vearth = get_body_barycentric_posvel('earth', t)[1]

    vsunearth = vearth #  - vsun

    _, gcrs_v = loc.get_gcrs_posvel(t)

    return vsunearth + gcrs_v

def helio_corr(t, loc, target):
    """
    Compute the correction required to convert a radial velocity at a given
    time and place to a heliocentric velocity.
    Paramters
    ---------
    t : astropy.time.Time
        Time of the observation. Can be a Time array.
    loc : astropy.coordinates.EarthLocation
        The observer location at which to compute the correction.
    target : astropy.coordinates.SkyCoord
        The on-sky location at which to compute the correction.
    Returns
    -------
    vcorr : astropy.units.Quantity with velocity units
        The heliocentric correction with a positive sign.  I.e., *add* this
        to an observed radial velocity to get the heliocentric velocity.
    """
    vsuntarg_cartrepr = helio_vector(t, loc)
    gcrs_p, _ = loc.get_gcrs_posvel(t)

    gtarg = target.transform_to(GCRS(obstime=t, obsgeoloc=gcrs_p))
    targcart = gtarg.represent_as(UnitSphericalRepresentation).to_cartesian()
    res = targcart.dot(vsuntarg_cartrepr).to(KPS)
    return res
    
def velcorr(time, skycoord, location=None):
  """Barycentric velocity correction.
  
  Uses the ephemeris set with  ``astropy.coordinates.solar_system_ephemeris.set`` for corrections. 
  For more information see `~astropy.coordinates.solar_system_ephemeris`.
  
  Parameters
  ----------
  time : `~astropy.time.Time`
    The time of observation.
  skycoord: `~astropy.coordinates.SkyCoord`
    The sky location to calculate the correction for.
  location: `~astropy.coordinates.EarthLocation`, optional
    The location of the observatory to calculate the correction for.
    If no location is given, the ``location`` attribute of the Time
    object is used
    
  Returns
  -------
  vel_corr : `~astropy.units.Quantity`
    The velocity correction to convert to Barycentric velocities. Should be added to the original
    velocity.
  """
  
  if location is None:
    if time.location is None:
        raise ValueError('An EarthLocation needs to be set or passed '
                         'in to calculate bary- or heliocentric '
                         'corrections')
    location = time.location
    
  # ensure sky location is ICRS compatible
  if not skycoord.is_transformable_to(ICRS()):
    raise ValueError("Given skycoord is not transformable to the ICRS")
  
  ep, ev = solar_system.get_body_barycentric_posvel('earth', t) # ICRS position and velocity of Earth's geocenter
  op, ov = location.get_gcrs_posvel(t) # GCRS position and velocity of observatory
  # ICRS and GCRS are axes-aligned. Can add the velocities
  velocity = ev + ov # relies on PR5434 being merged
  
  # get unit ICRS vector in direction of SkyCoord
  sc_cartesian = skycoord.icrs.represent_as(UnitSphericalRepresentation).represent_as(CartesianRepresentation)
  return sc_cartesian.dot(velocity).to(u.km/u.s) # similarly requires PR5434
    
t = Time(2457784.960783724, format='jd')
lapalma = coord.EarthLocation.of_site('lapalma')
m31 = coord.SkyCoord.from_name('M31')
epoch = Time("J2000").jd

res = bvc(t.jd, m31.ra.deg, m31.dec.deg, lat=lapalma.latitude.deg, lon=lapalma.longitude.deg,
          elevation=lapalma.height.value, epoch=epoch)
print((res/1000)*KPS)
          
print(helio_corr(t, lapalma, m31))
print(velcorr(t, m31, lapalma))