# astrometry.py
# Defines Astrometry timing model class
import numpy
import astropy.coordinates as coords
import astropy.units as u
import astropy.constants as const
from astropy.coordinates.angles import Angle
from astropy.coordinates.angles import rotation_matrix
from astropy import log
from . import parameter as p
from .timing_model import TimingModel, MissingParameter, Cache
from ..utils import time_from_mjd_string, time_to_longdouble, str2longdouble
from pint.pulsar_ecliptic import PulsarEcliptic, OBL
from pint import ls
from pint import utils
import time
mas_yr = (u.mas / u.yr)
try:
from astropy.erfa import DAYSEC as SECS_PER_DAY
except ImportError:
from astropy._erfa import DAYSEC as SECS_PER_DAY
[docs]class Astrometry(TimingModel):
def __init__(self):
super(Astrometry, self).__init__()
self.add_param(p.MJDParameter(name="POSEPOCH",
description="Reference epoch for position"))
self.add_param(p.floatParameter(name="PX",
units="mas", value=0.0,
description="Parallax"))
self.delay_funcs['L1'] += [self.solar_system_geometric_delay,]
[docs] def setup(self):
super(Astrometry, self).setup()
self.register_deriv_funcs(self.d_delay_astrometry_d_PX, 'delay', 'PX')
@Cache.cache_result
[docs] def ssb_to_psb_xyz(self, epoch=None):
"""Returns unit vector(s) from SSB to pulsar system barycenter.
If epochs (MJD) are given, proper motion is included in the calculation.
"""
# TODO: would it be better for this to return a 6-vector (pos, vel)?
return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose()
@Cache.cache_result
[docs] def barycentric_radio_freq(self, toas):
"""Return radio frequencies (MHz) of the toas corrected for Earth motion"""
L_hat = self.ssb_to_psb_xyz(epoch=toas['tdbld'].astype(numpy.float64))
v_dot_L_array = numpy.sum(toas['ssb_obs_vel']*L_hat, axis=1)
return toas['freq'] * (1.0 - v_dot_L_array / const.c)
[docs] def solar_system_geometric_delay(self, toas):
"""Returns geometric delay (in sec) due to position of site in
solar system. This includes Roemer delay and parallax.
NOTE: currently assumes XYZ location of TOA relative to SSB is
available as 3-vector toa.xyz, in units of light-seconds.
"""
L_hat = self.ssb_to_psb_xyz(epoch=toas['tdbld'].astype(numpy.float64))
re_dot_L = numpy.sum(toas['ssb_obs_pos']*L_hat, axis=1)
delay = -re_dot_L.to(ls).value
if self.PX.value != 0.0 \
and numpy.count_nonzero(toas['ssb_obs_pos']) > 0:
L = ((1.0 / self.PX.value) * u.kpc)
# TODO: numpy.sum currently loses units in some cases...
re_sqr = numpy.sum(toas['ssb_obs_pos']**2, axis=1) * toas['ssb_obs_pos'].unit**2
delay += (0.5 * (re_sqr / L) * (1.0 - re_dot_L**2 / re_sqr)).to(ls).value
return delay
#@Cache.use_cache
[docs] def get_d_delay_quantities(self, toas):
"""Calculate values needed for many d_delay_d_param functions """
# TODO: Move all these calculations in a separate class for elegance
rd = dict()
# TODO: Should delay not have units of u.second?
delay = self.delay(toas)
# TODO: toas['tdbld'].quantity should have units of u.day
# NOTE: Do we need to include the delay here?
rd['epoch'] = toas['tdbld'].quantity * u.day #- delay * u.second
# Distance from SSB to observatory, and from SSB to psr
ssb_obs = toas['ssb_obs_pos'].quantity
ssb_psr = self.ssb_to_psb_xyz(epoch=numpy.array(rd['epoch']))
# Cartesian coordinates, and derived quantities
rd['ssb_obs_r'] = numpy.sqrt(numpy.sum(ssb_obs**2, axis=1))
rd['ssb_obs_z'] = ssb_obs[:,2]
rd['ssb_obs_xy'] = numpy.sqrt(ssb_obs[:,0]**2 + ssb_obs[:,1]**2)
rd['ssb_obs_x'] = ssb_obs[:,0]
rd['ssb_obs_y'] = ssb_obs[:,1]
rd['in_psr_obs'] = numpy.sum(ssb_obs * ssb_psr, axis=1)
# Earth right ascension and declination
rd['earth_dec'] = numpy.arctan2(rd['ssb_obs_z'], rd['ssb_obs_xy'])
rd['earth_ra'] = numpy.arctan2(rd['ssb_obs_y'], rd['ssb_obs_x'])
return rd
@Cache.use_cache
[docs] def d_delay_astrometry_d_PX(self, toas):
"""Calculate the derivative wrt PX
Roughly following Smart, 1977, chapter 9.
px_r: Extra distance to Earth, wrt SSB, from pulsar
r_e: Position of earth (vector) wrt SSB
u_p: Unit vector from SSB pointing to pulsar
t_d: Parallax delay
c: Speed of light
delta: Parallax
The parallax delay is due to a distance orthogonal to the line of sight
to the pulsar from the SSB:
px_r = sqrt( r_e**2 - (r_e.u_p)**2 ),
with delay
t_d = 0.5 * px_r * delta'/ c, and delta = delta' * px_r / (1 AU)
"""
rd = self.get_d_delay_quantities(toas)
px_r = numpy.sqrt(rd['ssb_obs_r']**2-rd['in_psr_obs']**2)
dd_dpx = 0.5*(px_r**2 / (u.AU*const.c)) * (u.mas / u.radian)
# We want to return sec / mas
return dd_dpx.decompose(u.si.bases) / u.mas
#@Cache.cache_result
[docs] def d_delay_astrometry_d_POSEPOCH(self, toas):
"""Calculate the derivative wrt POSEPOCH
"""
pass
[docs]class AstrometryEquatorial(Astrometry):
def __init__(self):
super(AstrometryEquatorial, self).__init__()
self.add_param(p.AngleParameter(name="RAJ",
units="H:M:S",
description="Right ascension (J2000)",
aliases=["RA"]))
self.add_param(p.AngleParameter(name="DECJ",
units="D:M:S",
description="Declination (J2000)",
aliases=["DEC"]))
self.add_param(p.floatParameter(name="PMRA",
units="mas/year", value=0.0,
description="Proper motion in RA"))
self.add_param(p.floatParameter(name="PMDEC",
units="mas/year", value=0.0,
description="Proper motion in DEC"))
self.set_special_params(['RAJ', 'DECJ', 'PMRA', 'PMDEC'])
[docs] def setup(self):
super(AstrometryEquatorial, self).setup()
# RA/DEC are required
for p in ("RAJ", "DECJ"):
if getattr(self, p).value is None:
raise MissingParameter("Astrometry", p)
# If PM is included, check for POSEPOCH
if self.PMRA.value != 0.0 or self.PMDEC.value != 0.0:
if self.POSEPOCH.quantity is None:
if self.PEPOCH.quantity is None:
raise MissingParameter("AstrometryEquatorial", "POSEPOCH",
"POSEPOCH or PEPOCH are required if PM is set.")
else:
self.POSEPOCH.quantity = self.PEPOCH.quantity
self.register_deriv_funcs(self.d_delay_astrometry_d_RAJ, 'delay', 'RAJ')
self.register_deriv_funcs(self.d_delay_astrometry_d_DECJ, 'delay', 'DEC')
self.register_deriv_funcs(self.d_delay_astrometry_d_PMRA, 'delay', 'PMRA')
self.register_deriv_funcs(self.d_delay_astrometry_d_PMDEC, 'delay', 'PMDEC')
#@Cache.cache_result
[docs] def coords_as_ICRS(self, epoch=None):
"""Returns pulsar sky coordinates as an astropy ICRS object instance.
If epoch (MJD) is specified, proper motion is included to return
the position at the given epoch.
If the ecliptic coordinates are provided,
"""
if epoch is None or (self.PMRA.value == 0.0 and self.PMDEC.value == 0.0):
return coords.ICRS(ra=self.RAJ.quantity, dec=self.DECJ.quantity)
else:
dt = (epoch - self.POSEPOCH.quantity.mjd) * u.d
dRA = dt * self.PMRA.quantity / numpy.cos(self.DECJ.quantity.radian)
dDEC = dt * self.PMDEC.quantity
return coords.ICRS(ra=self.RAJ.quantity+dRA, dec=self.DECJ.quantity+dDEC)
#@Cache.use_cache
[docs] def d_delay_astrometry_d_RAJ(self, toas):
"""Calculate the derivative wrt RAJ
For the RAJ and DEC derivatives, use the following approximate model for
the pulse delay. (Inner-product between two Cartesian vectors)
de = Earth declination (wrt SSB)
ae = Earth right ascension
dp = pulsar declination
aa = pulsar right ascension
r = distance from SSB to Earh
c = speed of light
delay = r*[cos(de)*cos(dp)*cos(ae-aa)+sin(de)*sin(dp)]/c
"""
rd = self.get_d_delay_quantities(toas)
psr_ra = self.RAJ.quantity
psr_dec = self.DECJ.quantity
geom = numpy.cos(rd['earth_dec'])*numpy.cos(psr_dec)*\
numpy.sin(psr_ra-rd['earth_ra'])
dd_draj = rd['ssb_obs_r'] * geom / (const.c * u.radian)
return dd_draj.decompose(u.si.bases)
#@Cache.use_cache
[docs] def d_delay_astrometry_d_DECJ(self, toas):
"""Calculate the derivative wrt DECJ
Definitions as in d_delay_d_RAJ
"""
rd = self.get_d_delay_quantities(toas)
psr_ra = self.RAJ.quantity
psr_dec = self.DECJ.quantity
geom = numpy.cos(rd['earth_dec'])*numpy.sin(psr_dec)*\
numpy.cos(psr_ra-rd['earth_ra']) - numpy.sin(rd['earth_dec'])*\
numpy.cos(psr_dec)
dd_ddecj = rd['ssb_obs_r'] * geom / (const.c * u.radian)
return dd_ddecj.decompose(u.si.bases)
#@Cache.use_cache
[docs] def d_delay_astrometry_d_PMRA(self, toas):
"""Calculate the derivative wrt PMRA
Definitions as in d_delay_d_RAJ. Now we have a derivative in mas/yr for
the pulsar RA
"""
rd = self.get_d_delay_quantities(toas)
psr_ra = self.RAJ.quantity
te = rd['epoch'] - time_to_longdouble(self.POSEPOCH.quantity) * u.day
geom = numpy.cos(rd['earth_dec'])*numpy.sin(psr_ra-rd['earth_ra'])
deriv = rd['ssb_obs_r'] * geom * te / (const.c * u.radian)
dd_dpmra = deriv * u.mas / u.year
# We want to return sec / (mas / yr)
return dd_dpmra.decompose(u.si.bases) / (u.mas / u.year)
#@Cache.use_cache
[docs] def d_delay_astrometry_d_PMDEC(self, toas):
"""Calculate the derivative wrt PMDEC
Definitions as in d_delay_d_RAJ. Now we have a derivative in mas/yr for
the pulsar DEC
"""
rd = self.get_d_delay_quantities(toas)
psr_ra = self.RAJ.quantity
psr_dec = self.DECJ.quantity
te = rd['epoch'] - time_to_longdouble(self.POSEPOCH.quantity) * u.day
geom = numpy.cos(rd['earth_dec'])*numpy.sin(psr_dec)*\
numpy.cos(psr_ra-rd['earth_ra']) - numpy.cos(psr_dec)*\
numpy.sin(rd['earth_dec'])
deriv = rd['ssb_obs_r'] * geom * te / (const.c * u.radian)
dd_dpmdec = deriv * u.mas / u.year
# We want to return sec / (mas / yr)
return dd_dpmdec.decompose(u.si.bases) / (u.mas / u.year)
[docs]class AstrometryEcliptic(Astrometry):
def __init__(self):
super(AstrometryEcliptic, self).__init__()
self.add_param(p.AngleParameter(name="ELONG",
units="deg",
description="Ecliptic longitude",
aliases=["LAMBDA"]))
self.add_param(p.AngleParameter(name="ELAT",
units="deg",
description="Ecliptic latitude",
aliases=["BETA"]))
self.add_param(p.floatParameter(name="PMELONG",
units="mas/year", value=0.0,
description="Proper motion in ecliptic longitude",
aliases=["PMLAMBDA"]))
self.add_param(p.floatParameter(name="PMELAT",
units="mas/year", value=0.0,
description="Proper motion in ecliptic latitude",
aliases=["PMBETA"]))
self.add_param(p.strParameter(name="ECL",
description="Obliquity angle value secetion"))
self.set_special_params(['ELONG', 'ELAT', 'PMELONG','PMELAT'])
[docs] def setup(self):
super(AstrometryEcliptic, self).setup()
# RA/DEC are required
for p in ("ELONG", "ELAT"):
if getattr(self, p).value is None:
raise MissingParameter("AstrometryEcliptic", p)
# If PM is included, check for POSEPOCH
if self.PMELONG.value != 0.0 or self.PMELAT.value != 0.0:
if self.POSEPOCH.quantity is None:
if self.PEPOCH.quantity is None:
raise MissingParameter("Astrometry", "POSEPOCH",
"POSEPOCH or PEPOCH are required if PM is set.")
else:
self.POSEPOCH.quantity = self.PEPOCH.quantity
self.register_deriv_funcs(self.d_delay_astrometry_d_ELAT, 'delay', 'ELAT')
self.register_deriv_funcs(self.d_delay_astrometry_d_ELONG, 'delay', 'ELONG')
self.register_deriv_funcs(self.d_delay_astrometry_d_PMELAT, 'delay', 'PMELAT')
self.register_deriv_funcs(self.d_delay_astrometry_d_PMELONG, 'delay', 'PMELONG')
@Cache.cache_result
[docs] def coords_as_ICRS(self, epoch=None):
"""Returns pulsar sky coordinates as an astropy ICRS object instance.
Pulsar coordinates will be transform from ecliptic coordinates to ICRS
If epoch (MJD) is specified, proper motion is included to return
the position at the given epoch.
If the ecliptic coordinates are provided,
"""
if epoch is None or (self.PMELONG.value == 0.0 and self.PMELAT.value == 0.0):
pos_ecl = PulsarEcliptic(lon=self.ELONG.quantity, lat=self.ELAT.quantity)
else:
dt = (epoch - self.POSEPOCH.quantity.mjd) * u.d
dELONG = dt * self.PMELONG.quantity / numpy.cos(self.ELAT.quantity.radian)
dELAT = dt * self.PMELAT.quantity
try:
PulsarEcliptic.obliquity = OBL[self.ECL.value]
except KeyError:
raise ValueError("No obliquity " + self.ECL.value + " provided. "
"Check your pint/datafile/ecliptic.dat file.")
pos_ecl = PulsarEcliptic(lon=self.ELONG.quantity+dELONG, lat=self.ELAT.quantity+dELAT)
return pos_ecl.transform_to(coords.ICRS)
[docs] def get_d_delay_quantities_ecliptical(self, toas):
"""Calculate values needed for many d_delay_d_param functions """
# TODO: Move all these calculations in a separate class for elegance
rd = dict()
# From the earth_ra dec to earth_elong and elat
try:
PulsarEcliptic.obliquity = OBL[self.ECL.value]
except KeyError:
raise ValueError("No obliquity " + self.ECL.value + " provided. "
"Check your pint/datafile/ecliptic.dat file.")
rd = self.get_d_delay_quantities(toas)
coords_icrs = coords.ICRS(ra=rd['earth_ra'], dec=rd['earth_dec'])
coords_elpt = coords_icrs.transform_to(PulsarEcliptic)
rd['earth_elong'] = coords_elpt.lon
rd['earth_elat'] = coords_elpt.lat
return rd
#@Cache.use_cache
[docs] def d_delay_astrometry_d_ELONG(self, toas):
"""Calculate the derivative wrt RAJ
For the RAJ and DEC derivatives, use the following approximate model for
the pulse delay. (Inner-product between two Cartesian vectors)
de = Earth declination (wrt SSB)
ae = Earth right ascension
dp = pulsar declination
aa = pulsar right ascension
r = distance from SSB to Earh
c = speed of light
delay = r*[cos(de)*cos(dp)*cos(ae-aa)+sin(de)*sin(dp)]/c
elate = Earth elat (wrt SSB)
elonge = Earth elong
elatp = pulsar elat
elongp = pulsar elong
r = distance from SSB to Earh
c = speed of light
delay = r*[cos(elate)*cos(elatp)*cos(elonge-elongp)+sin(elate)*sin(elatp)]/c
"""
rd = self.get_d_delay_quantities_ecliptical(toas)
psr_elong = self.ELONG.quantity
psr_elat = self.ELAT.quantity
geom = numpy.cos(rd['earth_elat'])*numpy.cos(psr_elat)*\
numpy.sin(psr_elong-rd['earth_elong'])
dd_delong = rd['ssb_obs_r'] * geom / (const.c * u.radian)
return dd_delong.decompose(u.si.bases)
#@Cache.use_cache
[docs] def d_delay_astrometry_d_ELAT(self, toas):
"""Calculate the derivative wrt DECJ
Definitions as in d_delay_d_RAJ
"""
rd = self.get_d_delay_quantities_ecliptical(toas)
psr_elong = self.ELONG.quantity
psr_elat = self.ELAT.quantity
geom = numpy.cos(rd['earth_elat'])*numpy.sin(psr_elat)*\
numpy.cos(psr_elong-rd['earth_elong']) - numpy.sin(rd['earth_elat'])*\
numpy.cos(psr_elat)
dd_delat = rd['ssb_obs_r'] * geom / (const.c * u.radian)
return dd_delat.decompose(u.si.bases)
#@Cache.use_cache
[docs] def d_delay_astrometry_d_PMELONG(self, toas):
"""Calculate the derivative wrt PMRA
Definitions as in d_delay_d_RAJ. Now we have a derivative in mas/yr for
the pulsar RA
"""
rd = self.get_d_delay_quantities_ecliptical(toas)
psr_elong = self.ELONG.quantity
psr_elat = self.ELAT.quantity
te = rd['epoch'] - time_to_longdouble(self.POSEPOCH.quantity) * u.day
geom = numpy.cos(rd['earth_elat'])*numpy.sin(psr_elong-rd['earth_elong'])
deriv = rd['ssb_obs_r'] * geom * te / (const.c * u.radian)
dd_dpmelong = deriv * u.mas / u.year
# We want to return sec / (mas / yr)
return dd_dpmelong.decompose(u.si.bases) / (u.mas / u.year)
#@Cache.use_cache
[docs] def d_delay_astrometry_d_PMELAT(self, toas):
"""Calculate the derivative wrt PMDEC
Definitions as in d_delay_d_RAJ. Now we have a derivative in mas/yr for
the pulsar DEC
"""
rd = self.get_d_delay_quantities_ecliptical(toas)
psr_elong = self.ELONG.quantity
psr_elat = self.ELAT.quantity
te = rd['epoch'] - time_to_longdouble(self.POSEPOCH.quantity) * u.day
geom = numpy.cos(rd['earth_elat'])*numpy.sin(psr_elat)*\
numpy.cos(psr_elong-rd['earth_elong']) - numpy.cos(psr_elat)*\
numpy.sin(rd['earth_elat'])
deriv = rd['ssb_obs_r'] * geom * te / (const.c * u.radian)
dd_dpmelat = deriv * u.mas / u.year
# We want to return sec / (mas / yr)
return dd_dpmelat.decompose(u.si.bases) / (u.mas / u.year)