Skip to content

Instantly share code, notes, and snippets.

@StuartLittlefair
Last active July 18, 2023 17:00
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save StuartLittlefair/5aaf476c5d7b52d20aa9544cfaa936a1 to your computer and use it in GitHub Desktop.
Save StuartLittlefair/5aaf476c5d7b52d20aa9544cfaa936a1 to your computer and use it in GitHub Desktop.
How to do barycentric velocity correction with astropy v1.3
from astropy.time import Time
from astropy.coordinates import SkyCoord, solar_system, EarthLocation, ICRS
from astropy import units as u
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 = sc.icrs.represent_as(UnitSphericalRepresentation).represent_as(CartesianRepresentation)
return sc_cartesian.dot(velocity).to(u.km/u.s) # similarly requires PR5434
@janerigby
Copy link

janerigby commented Feb 15, 2017

Stuart, in line 46, sc. isn't defined. assuming it should be "skycoord', then I get this error:
"global name 'UnitSphericalRepresentation' is not defined".
Needed to add "from astropy import coordinates"

@telegraphic
Copy link

Thanks, this code ran once I renamed sc to skycoord and added an import line
from astropy.coordinates import UnitSphericalRepresentation, CartesianRepresentation.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment