search

Libsharp / libsharp

Library for fast spherical harmonic transforms, see http://arxiv.org/abs/1303.4945
GNU General Public License v2.0
24 stars 17 forks source link

Normalization of alms #23

Open zonca opened 5 years ago

zonca commented 5 years ago

I am trying to read alms I created with healpy from disk and transform them in the format expected by libsharp. However it seems there is a different normalization between the alm I create with healpy and with libsharp. Investigating this issue I even find some strange behavior, could you please take a look at this notebook and check if you find any mistake/algorithmic errors?

https://nbviewer.jupyter.org/gist/zonca/a481471543cd6f658dea5bd48d0f94e9

mreineck commented 5 years ago

If I understand correctly, you are requesting real-valued a_lm from libsharp:

libsharp_order = libsharp.packed_real_order(lmax, ms=local_m_indices) I don't understand how the analysis call can return complex-valued a_lm in such a case. (Disclaimer: I don't know the current Python interface of libsharp at all.) In any case: for real-valued a_lm all entries with m!=0 are scaled by a factor sqrt(2), so this is probably the discrepancy you are seeing.

mreineck commented 5 years ago

Sorry, I misread the notebook... the libsharp a_lm are not complex-valued after all.

For a description of the normalization, see e.g. https://commander.readthedocs.io/en/latest/library_guide/sphere.html

zonca commented 5 years ago

thanks, now normalization is clear. What about that signal at m=1 which only appears in libsharp and not in HEALPix?

mreineck commented 5 years ago

I can't try this myself at the moment, but I think you are interpreting the storage scheme of the real package a_lm incorrectly. They are storing only the values that can be different from zero, i.e. the imaginary value for (l=0, m) is left out, leading to a data set that has minimal size, but can be quite complicated to interpret.

mreineck commented 5 years ago

To be explicit, the storage order for a given m is Re(alm(0,m)), Re(alm(1,m)), Im(alm(1,m)), Re(alm(2,m))...

dagss commented 5 years ago

@zonca what Martin says is correct, but this discussion seems a bit on the head for me.

The real-packed storage is very convenient for some things, but very different from what Healpix-the-software is using.

If you want the Healpix storage scheme then libsharp certainly supports that directly. Perhaps support is needed for that in the Python wrapper though I cannot remember.

zonca commented 5 years ago

I didn't know libsharp supported the HEALPix storage scheme, I see there is a triangular and a rectangular order available in Python. Is one of them the same as HEALPix?

mreineck commented 5 years ago

Yes, triangular.