Skip to content

Backport light deflection #10885

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
eteq merged 2 commits into from
Oct 19, 2020
Merged

Backport light deflection #10885

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
258e603
Revert "Merge pull request #10666 from StuartLittlefair/replace_ldsun"
mhvk Oct 16, 2020
7745d48
Better light deflection calcs in solar system
StuartLittlefair Aug 18, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
30 changes: 26 additions & 4 deletions astropy/coordinates/builtin_frames/icrs_cirs_transforms.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,11 @@ def icrs_to_cirs(icrs_coo, cirs_frame):
def cirs_to_icrs(cirs_coo, icrs_frame):
# set up the astrometry context for ICRS<->cirs and then convert to
# astrometric coordinate direction
astrom = erfa_astrom.get().apci(cirs_coo)
jd1, jd2 = get_jd12(cirs_coo.obstime, 'tt')
x, y, s = get_cip(jd1, jd2)
earth_pv, earth_heliocentric = prepare_earth_position_vel(cirs_coo.obstime)
# erfa.apci requests TDB but TT can be used instead of TDB without any significant impact on accuracy
astrom = erfa.apci(jd1, jd2, earth_pv, earth_heliocentric, x, y, s)
srepr = cirs_coo.represent_as(SphericalRepresentation)
i_ra, i_dec = aticq(srepr.without_differentials(), astrom)

Expand Down Expand Up @@ -109,8 +113,20 @@ def cirs_to_cirs(from_coo, to_frame):

@frame_transform_graph.transform(FunctionTransformWithFiniteDifference, ICRS, GCRS)
def icrs_to_gcrs(icrs_coo, gcrs_frame):
# first set up the astrometry context for ICRS<->GCRS.
astrom = erfa_astrom.get().apcs(gcrs_frame)
# first set up the astrometry context for ICRS<->GCRS. There are a few steps...
# get the position and velocity arrays for the observatory. Need to
# have xyz in last dimension, and pos/vel in one-but-last.
# (Note could use np.stack once our minimum numpy version is >=1.10.)
obs_pv = erfa.pav2pv(
gcrs_frame.obsgeoloc.get_xyz(xyz_axis=-1).to_value(u.m),
gcrs_frame.obsgeovel.get_xyz(xyz_axis=-1).to_value(u.m/u.s))

# find the position and velocity of earth
jd1, jd2 = get_jd12(gcrs_frame.obstime, 'tdb')
earth_pv, earth_heliocentric = prepare_earth_position_vel(gcrs_frame.obstime)

# get astrometry context object, astrom.
astrom = erfa.apcs(jd1, jd2, obs_pv, earth_pv, earth_heliocentric)

if icrs_coo.data.get_name() == 'unitspherical' or icrs_coo.data.to_cartesian().x.unit == u.one:
# if no distance, just do the infinite-distance/no parallax calculation
Expand Down Expand Up @@ -202,7 +218,13 @@ def gcrs_to_hcrs(gcrs_coo, hcrs_frame):

# set up the astrometry context for ICRS<->GCRS and then convert to ICRS
# coordinate direction
astrom = erfa_astrom.get().apcs(gcrs_coo)
obs_pv = erfa.pav2pv(
gcrs_coo.obsgeoloc.get_xyz(xyz_axis=-1).to_value(u.m),
gcrs_coo.obsgeovel.get_xyz(xyz_axis=-1).to_value(u.m/u.s))

jd1, jd2 = get_jd12(hcrs_frame.obstime, 'tdb')
earth_pv, earth_heliocentric = prepare_earth_position_vel(gcrs_coo.obstime)
astrom = erfa.apcs(jd1, jd2, obs_pv, earth_pv, earth_heliocentric)
srepr = gcrs_coo.represent_as(SphericalRepresentation)
i_ra, i_dec = aticq(srepr.without_differentials(), astrom)

Expand Down
8 changes: 4 additions & 4 deletions astropy/coordinates/builtin_frames/utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -223,8 +223,8 @@ def aticq(srepr, astrom):
eh = astrom['em'][..., np.newaxis] * astrom['eh']
# unit vector from Sun to object
q = eh + srepr_distance[..., np.newaxis].to_value(u.au) * before
sundist = np.linalg.norm(q, axis=-1, keepdims=True)
q /= sundist
sundist, q = erfa.pn(q)
sundist = sundist[..., np.newaxis]
# calculation above is extremely unstable very close to the sun
# in these situations, default back to ldsun-style behaviour,
# since this is reversible and drops to zero within stellar limb
Expand Down Expand Up @@ -294,8 +294,8 @@ def atciqz(srepr, astrom):
eh = astrom['em'][..., np.newaxis] * astrom['eh']
# unit vector from Sun to object
q = eh + srepr_distance[..., np.newaxis].to_value(u.au) * pco
sundist = np.linalg.norm(q, axis=-1, keepdims=True)
q /= sundist
sundist, q = erfa.pn(q)
sundist = sundist[..., np.newaxis]
# calculation above is extremely unstable very close to the sun
# in these situations, default back to ldsun-style behaviour,
# since this is reversible and drops to zero within stellar limb
Expand Down