Inconsistent distance calculations in dm_to_dist

[shutsch]

Hi,

It appears to me that there are some inconsistencies with the dm_to_dist function in pygedm.
For some particular lines of sights and radial positions, the predicted distances decrease with increasing DM, which is clearly unphysical, as it would correspond to negative electron densities between the points.

An example script:

import pygedm

gl = 243.49
gb = 45.782

dm_small = 11.21
dm_large = 11.22

print(pygedm.dm_to_dist(gl, gb, dm_small)[0])
print(pygedm.dm_to_dist(gl, gb, dm_large)[0])

This gives:

817.3250122070312 pc
803.8787231445312 pc 

The distance starts to increase again with rising DM from there on.
I have found a couple of dozens of such jumps, I can attach a list if needed.

The electron densities at these points are physical and (roughly) consistent with the increase in DM.

[telegraphic]

Hi @shutsch, that is indeed peculiar!

This behavior matches the compiled version of YMW16, so it is a quirk/bug of YMW16 itself:

./ymw16 Gal 243.49 45.782 11.21 1
Warning: YMW16_DIR set to local directory
Gal: gl= 243.490 gb=  45.782 DM=   11.21 DM_Gal:   11.21 Dist:    817.3 log(tau_sc): -7.983

./ymw16 Gal 243.49 45.782 11.22 1
Warning: YMW16_DIR set to local directory
Gal: gl= 243.490 gb=  45.782 DM=   11.22 DM_Gal:   11.22 Dist:    803.9 log(tau_sc): -7.982

Having a bit of a diagnostic look, it seems this is due to different step sizes being chosen for the integration along the line of sight within YMW16:

dm=11.21: dmstep=0.047034, ncount=191, dstep=4.311538
dm=11.22:  dmstep=0.047416, ncount=188, dstep=4.315385

Not sure if this is a bug or just a quirk of the numerical integration. But my feeling is that this is probably just a numerical error, and other uncertainties due to the model itself are larger.

YMW16 does not output uncertainty estimates, but in the pygedm paper I calculated that 87% of distance estimates lie
between 0.35-2.76x the actual distance (i.e. measured distances to pulsars through parallax or other techniques). Based on this, you could argue that the uncertainty bounds for a DM of 11.22 pc/cm3 would be 281--2216 pc. So the numerical error (if that is what it is) is much lower than the model error!

PS: I considered adding uncertainties to pygedm, but in practice errors are strongly dependent upon line of sight in a complex way. A target that is close to a pulsar used when fitting the model will have a more accurate estimate more accurate than a target that isn't near a pulsar used to fit the model.

[shutsch]

Hi @telegraphic,

Thanks for looking into this.
I agree that a numerical error in the integration is a plausible explanation.
I ran into the issue, as I use pygedm in a numerical optimization scheme which relies on the DM-distance relation being monotonously increasing.
This crashed when encountering these jumps, but I was able to alleviate this by smoothing the values given the dm_to_dist function.

I was a bit worried about systematic effects this would introduce, but the paper you reference makes it clear that this isn't really an issue.
If this is not worth fixing, I think at least having a warning somewhere in the documentation would be helpful!

[telegraphic]

Updating docs is a good plan 👍