Update NSIDE_MAX

[andrew-saydjari]

I am using this for a photometry project where the base resolution of the images is 0.26 arcsec (2k x 4k detectors). So, reasonable maps really need at least NSIDE = 2^14 in order to visualize CCD scale features. It is not obvious how to increase NSIDE_MAX from the user side? If it is easy to change, maybe adding to the docs to tell users how to do this would be sufficient. Otherwise, I would appreciate merging this PR to increase the max to at least 2^15. The code seems to have been written well enough that it generalizes just fine and nothing breaks when I changed this on my fork. However, let me know if there is some functionality I missed that I would need to edit to enable this change.

[ziotom78]

Hi @andrew-saydjari , the reason why Healpix.jl set the maximum value for NSIDE to 8192 is because this is the same constraint that has historically been set on other Healpix implementations.

This limit is due to the fact that when NSIDE>8192, the index of the last pixel in the map cannot be stored in a 32-bit integer, and thus 32-bit machines cannot use functions like ang2pix. (It is the reason why Healpy stopped supporting 32-bit machines).

Since Healpix.jl use the Int type to index pixels in a map, the general rule is that the following condition be always valid:

log2(nside2npix(Int128(2)^i)) < (sizeof(Int) * 8 - 1)

This is the reason why at the moment i must be ≤13.

Perhaps we might adapt the value of NSIDE_MAX to the size of an Int, so that on 32-bit architecture is set to 2^13, and on 64-bit architecture is set to 2^29. In this way, Healpix.jl would transparently support codes that use higher resolutions on 64-bit machines, while at the same time refusing to run these codes on older architectures:

NSIDE_MAX = 2^(floor(Int, 0.5 * ((sizeof(Int) * 8 - 1) - log2(12))))

@mreineck, what do you think?

[andrew-saydjari]

I pushed the requested change. Could a maintainer approve running the workflows?

[andrew-saydjari]

It looks like datatables needs to be modified to handle the change in NSIDE max, especially for ang2pixNest which references a lot of its index calculations to NSIDE_MAX.

[hsgg]

For readability, may I suggest using typemax()? I.e.,

NSIDE_MAX = 2^(floor(Int, 0.5 * log2(typemax(Int)) - log2(12)))

Also, I'm slightly confused why it wouldn't be

NSIDE_MAX = 2^floor(Int, 0.5 * log2(typemax(Int) / 12))

e.g.,

julia> 2^floor(Int, 0.5 * log2(typemax(Int32) / 12))
8192

[ziotom78]

I need to check how those data tables were produced, this might take some time because at the moment I am busy with end-of-semester duties.

BTW, @andrew-saydjari, why do you need a spherical projection like Healpix for representing those CCD images? Isn't a planar projection enough for your purpose? I'm curious!

[andrew-saydjari]

@ziotom78 I would help out, but I think it would be easier for someone with access to how they were produced. I am happy to wait until you have some more time.

Each individual CCD image is probably best viewed as a planar projection, but we have a very large survey area (>12% of the sky) on which we want to the survey coverage (and other derived stats). However, in order to do that at a resolution where gaps between CCDs and the like are visible, a high res Healpix map seems ideal in order to both wrangle projection distortions for the large survey footprint and maintain the details.

[ziotom78]

I am closing this because I was able to refactor the code of the library and add support for NSIDE>8192 in #75 . Let's track the progress there.