Sign error in solar angle calculation

[Matthew Kerr]

There seems to be a sign error in the calculation of the solar angle (angle between sun and pulsar as viewed by observer.) This results in an error of a few tenths of a degree (order the difference in SSB vs. heliocenter).

This came up in testing the PINT implementation comparing it with tempo2. Both using the back of a napkin and by comparison to a direct calculation with pyephem, I’ve convinced myself that the tempo2 calculation has a sign error. Below is a script and the resulting plot showing this, largely based on one written by Paul Ray and Scott Ransom.

The error appears both in the general2 plugin:

317c317
< rsa[j] = (-psr[0].obsn[varN].sun_ssb[j] - psr[0].obsn[varN].earth_ssb[j] + psr[0].obsn[varN].observatory_earth[j]);

> rsa[j] = psr[0].obsn[varN].sun_ssb[j] - psr[0].obsn[varN].earth_ssb[j] - psr[0].obsn[varN].observatory_earth[j];

and similarly in lines 70 and 72 of dm_delays.C for the solar wind delay calculation.

After applying this patch, the calculation in PINT and tempo2 is identical. Residual differences relative to pyephem are down at the arcminute level, and probably indicate a need for a refined test.

‌

[Michael Keith]

The definition of vectors in tempo2 seems somewhat suspect at the moment. Perhaps a full review is needed of what exactly these parameters are supposed to encode. I’m a bit nervous to change these things in case it comes up elsewhere, but I suspect that you are right.

In the code I have in dm_delays (which is more important than general2) it has
rsa[j] = -psr[p].obsn[i].sun_ssb[j] + psr[p].obsn[i].earth_ssb[j] + psr[p].obsn[i].observatory_earth[j];

Which doesn’t seem to match either of your lines… which is a bit confusing. Your “correct” line seems to be the the same but negative of the one in dm_delays.C, but… wouldn’t that put the solar wind at entirely the wrong phase (the observatory->sun vector would point totally backwards, so surely would be out of phase by 6 months).

I’ll try to investigate a bit.

Indeed, general2 to use the logic from dm_delays.C gives the following output from our script. I suspect that there must be a second minus sign somewhere that undoes this bug or otherwise the solar wind computations would be totally wrong in tempo2.

(Or… there is some deeper bug which explains both this and the other sign error we recently found)

[Michael Keith]

So… plotting the actual dm delay calculated vs MJD, the peak is well aligned with the minumum angle printed using your proposed new code. I think this comes from the fact that the “solar angle” used for the DM model is defined in a non-intuitive way. Internally for the dm delay it uses a value ‘ctheta’, i.e. cos theta, but theta is not exactly the angle between the sun and the pulsar, rather this is the angle rho from the tempo2 paper equation 30 where the DM due to the solar wind is rho/(r.sin(rho)).

This is a maximum at pi radians, where rho is large and sin-rho is small, hence why the peak is at the largest negative value of ctheta (and therefore theta (or rho) close to 180 degrees).

Nevertheless I think that it IS a bug that general2 had negated some values but not others when printing the ‘solar angle’, and so I think your fix makes sense for general2 plugin, but the dm_delays are fine.

[Michael Keith]

Patch pushed to 3f987a1 to use Matthew's definition of solar angle in general2 plugin.