PintPulsar needs to pull SSB vectors from Pint

[Hazboun6]

Currently PintPulsar populates _get_planetssb and _get_sunssb with arrays of zeros. These are the vectors needed from the timing package in order to use BayesEphem and the solar wind models, respectively. These are available in Pint, but the code in enterprise needs to be changed in order to be able to do ephemeris and solar wind modeling using Pint as the basal timing package.

https://github.com/nanograv/enterprise/blob/master/enterprise/pulsar.py#L341

[Hazboun6]

@vallis While #244 added the needed vectors, @AaronDJohnson's testing has not been successful.

@AaronDJohnson can you resubmit a PR, and link it here, when you get a chance to investigate the issue up close?

[Hazboun6]

Okay, we've diagnosed one major issue, which is that the PintPulsar object seems to lack a meaningful pos_t attribute, which is needed for the BayesEphem model.

# TODO: pos_t not currently implemented
self._pos_t = np.zeros((len(self._toas), 3))

We need to define the vector from Pint, but I am not quite sure how to do that. The libstempo/Tempo2 version looks like this:

self._pos_t = t2pulsar.psrPos.copy()
if "ELONG" and "ELAT" in np.concatenate((t2pulsar.pars(which="fit"), t2pulsar.pars(which="set"))):
    self._pos_t = utils.ecl2eq_vec(self._pos_t)

@paulray @aarchiba @scottransom is there an easy way to access this? I didn't see anything in the table that seemed pertinent. But I might not be recognizing an obvious name.

[scottransom]

@Hazboun6 What is the pos_t attribute? Is that the topocentric position of the pulsar?

[Hazboun6]

@scottransom Yes, sorry for not saying that explicitly!

[scottransom]

I'm pretty sure that is the ssb_obs_pos in the TOA table. For instance:

In [5]: t.table['ssb_obs_pos'][0]
Out[5]: array([ 4.08342950e+07, -1.34788757e+08, -5.84439627e+07])

[Hazboun6]

Thanks! Just out of curiosity, how would I translate 'ssb_obs_pos'? I was thinking that was the position vector from the SSB to the observatory.

[vallis]

I do think the tempo2 version is a unit vector… you can normalize it in PintPulsar. Michele

…

On Oct 5, 2021, at 1:12 PM, Jeff Hazboun ***@***.***> wrote: Thanks! Just out of curiosity, how would I translate 'ssb_obs_pos'? I was thinking that was the position vector from the SSB to the observatory. — You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <#238 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AADYJME2HZ4LBO7MU6ID7BDUFNL4VANCNFSM4T6INWFA>. Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.

[scottransom]

Oh, sorry! I misstated above. Yes, 'ssb_obs_pos' is the topocentric position of the observatory not the pulsar. So if we need the vector for the pulsar, that is gotten from the astrometry modules. See the ssb_to_psb* methods in astrometry.py.

[Hazboun6]

@vallis I'm confused about Pulsar.pos and Pulsar.pos_t. The latter is just the same value repeated for NTOAs, but the unit vector in .pos is close but does not match the values in pos_t (in a TempoPulsar). In the code it looks like pos is calculated from the RA and DEC. I'm wondering if I can just repeat pos or repeat the values given by the functions that Scott referenced above? They also do not match in PINT. In both cases they are close, but slightly off.

[scottransom]

@Hazboun6 Maybe one of them is the position wrt the SSB and the other (that varies) is the pulsar position wrt the observatory? (although with proper motion, I'd think that both should vary!) Seems like those variables need to be better documented! :-)

[Hazboun6]

You don't think these descriptions would fly in an astrometry paper? ;-)

@property
def pos(self):
    """Return unit vector to pulsar."""
    return self._pos

@property
def pos_t(self):
    """Return unit vector to pulsar as function of time."""
    return self._pos_t[self._isort, :]

[Hazboun6]

and the other (that varies)

@scottransom That is why I am really confused, because pos_t doesn't vary. It is the same number at all times. That's why I don't understand why it is different.

[vallis]

Pulsar.pos is evaluated from the estimated timing-model parameters at the fiducial epoch. pos_t is supposed to include the effects of proper motion. Both from SSB.

Perhaps the proper-motion shift in the vector is so small that is vanishes in double precision? Or perhaps tempo2 is lying about doing that? (You can check by imposing a huge proper motion in the parameters.) Either way, since the pos is dotted into a very small ephemeris delta, I think the time dependence would be negligible. We probably included it because why not...

[Hazboun6]

Okay, that helps quite a bit. In fact when I look at the first and last values in the array, fully displayed, as opposed to the shorten version you see when calling the array, I do see differences in the values.

So for consistency I can then do the same with the PintPulsar.
@scottransom Am I right in thinking that to get the proper motion effects I can do:

ssb_to_psb_xyz_ICRS(toas.get_mjds())

We convert everything into equatorial coordinates, so it is nice that this ICRS method is in here!
Also, I am assuming that I do not need get_mjds(high_precision=True) correct?

[scottransom]

Basically. Although that method is part of whatever "Astrometry*" model component that you have. So with a model m, toas object t, and using Equatorial coords, you would do this with:

In [17]: m.components['AstrometryEquatorial'].ssb_to_psb_xyz_ICRS(t.get_mjds())
Out[17]: 
<Quantity [[-0.04721326, -0.90671017, -0.41910329],
           [-0.04721326, -0.90671017, -0.41910329],
           [-0.04721326, -0.90671017, -0.41910329],
           ...,
           [-0.04721327, -0.90671017, -0.41910329],
           [-0.04721327, -0.90671017, -0.41910329],
           [-0.04721327, -0.90671017, -0.41910329]]>

However, a couple points. Those are topocentric times, and not barycentric times which are probably the correct thing to have. I suspect that will have a very small effect, though. We can easily give you the barycentric times if you want them instead. For example:

In [19]: m.components['AstrometryEquatorial'].ssb_to_psb_xyz_ICRS(m.get_barycentric_toas(t))
Out[19]: 
<Quantity [[-0.04721326, -0.90671017, -0.41910329],
           [-0.04721326, -0.90671017, -0.41910329],
           [-0.04721326, -0.90671017, -0.41910329],
           ...,
           [-0.04721327, -0.90671017, -0.41910329],
           [-0.04721327, -0.90671017, -0.41910329],
           [-0.04721327, -0.90671017, -0.41910329]]>

[Hazboun6]

Great! Thanks @scottransom that is all very useful. The pulsar I had loaded (not a 15 yr) has 'AstrometryEquatorial', but it looks like both classes have thessb_to_psb_xyz_ICRS function. To keep this generic, should I use which ever Astrometry class is part of the model?

[scottransom]

Yeah, that's exactly right. Would be nice if those methods were abstracted out to be part of a generic "Astrometry" interface. Maybe there is a way to do that. Just pinging @luojing1211 to check.

[luojing1211]

Internally, all calculation is done under the ICRS coordinate. If the input is ecliptic coordinate, it will first convert the position to the ICRS coordinate.

[scottransom]

@luojing1211 Right. I was just wondering if there wasn't a way that could use a generic astrometry class, rather than the specific Equatorial or Ecliptic component. In otherwords, something like:
m.components['Astrometry'].ssb_to_psb_xyz_ICRS(...)
rather than having to check which component the module has and then using:
m.components['AstrometryEquatorial'].ssb_to_psb_xyz_ICRS(...)

[luojing1211]

@scottransom I see what you mean. That is a new feature request which is going to be very useful. We can have an issue on PINT github.

[AaronDJohnson]

The BE signals are now exactly the same between PINT and TEMPO2. Here is a plot of 5000 points of a 43 pulsar run (w/out J1713 and B1937) with no correlations with BE on the 12.5 year data.

Subtracting the two, they appear to have a constant offset. This offset appears both with and without BayesEphem.