RadioAstronomySoftwareGroup / pyuvsim

A ultra-high precision package for simulating radio interferometers in python on compute clusters.
https://pyuvsim.readthedocs.io/en/latest/
BSD 3-Clause "New" or "Revised" License
43 stars 7 forks source link

Does pyuvsim write telescope locations with sufficient precision? #400

Closed jpober closed 6 months ago

jpober commented 2 years ago

@StellanBechtold has been working to replicate the visibilities for a single source simulation using only the alt/az calculation in the SkyModel object and an analytic beam. He's found that when he loads in the UVData object that pyuvsim created, and uses get_telescope_location_lat_lon_alt_degrees to get the coordinates to initialize the SkyModel, he gets slightly different coordinates than what was used in the yaml file initializing pyuvsim. For example, if the yaml file has:

(-30.72152777777791, 21.428305555555557,1073.0000000093132)

The resultant UVData object instead has:

(-30.721527777778046, 21.428305555555557, 1073.0000000176951)

The effect when comparing the calculation with the simulation is clear. When the SkyModel is initialized using the values from get_telescope_location_lat_lon_alt_degrees the difference between his calculation and the simulated visibilities looks like this:

Screen Shot 2022-04-26 at 2 22 47 PM

If, instead, he copies the values directly from the yaml file to initialize the SkyModel, he gets residuals that look like this:

Screen Shot 2022-04-26 at 2 15 48 PM

I suspect that pyuvdata has been tested to correctly roundtrip telescope locations (but maybe not?), so the likely culprit seems to be that pyuvsim is not storing the telescope location with sufficient precision before setting it on the resultant UVData object. Or maybe there's something else going on.

bhazelton commented 2 years ago

What file type is he writing out to? This may be very file-type specific.

jpober commented 2 years ago

@StellanBechtold I think it's a uvh5 file, right?

jpober commented 2 years ago

Back of the envelope, this is sub-micron shift in the telescope position. There's no need for pyuvdata to support that level of precision, but if it's an easy fix in pyuvsim to cast the telescope location as a double when reading it out of the yaml, that's probably worth doing.

StellanBechtold commented 2 years ago

Yes, I only used uvh5 files.

bhazelton commented 6 months ago

I looked into this and casting to a double (float64) doesn't solve the problem. The change is caused by converting from lat/lon/alt to ECEF coordinates (which are what is actually stored on the UVData object) and then converting back. That conversion code is implemented in cython in pyuvdata's utils.pyx and it all seems to be double precision, but maybe @mkolopanis can confirm.

mkolopanis commented 6 months ago

a note as of this statement, uvdata's extension module only uses 64-bit precision in the geodetic <-> geocentric conversions