Compute Gaussian random field from power spectrum

[jasonmcewen]

We should have a function that takes a power spectrum and then compute a Gaussian random field with that power spectrum. While this is relatively straightforward, it would make things easier for users.

We already have functionality to compute Gaussian randoms fields on the sphere in signal_generator but this doesn't support a power spectrum. A good place for this functionality is probably in the same file.

[ntessore]

I'd be happy to give it a shot!

[ntessore]

Looking into implementing this. In principle, it's straightforward, since you have the flm generation routine already. However, s2fft generally works with the full set of coefficients $-\ell \le m \le \ell$ for the complex fields, not E and B modes. This requires some agreement on the API.

First, presumably you want to adopt a similar convention of using the "full" set of Cls, not E/B-mode Cls? This means a pair of complex-valued Cls for each combination of spin-weighted fields: variance $C_\ell^{EE} + C_\ell^{BB} + i (C_\ell^{EB} - C_\ell^{BE})$ and pseudo-variance $C_\ell^{EE} - C_\ell^{BB} + i (C_\ell^{EB} + C_\ell^{BE})$. (This would merit a helper function to pack/unpack E/B-mode Cls.)

In this case, the most natural input for generating $N$ complex-valued fields would be a stack of complex-valued Cls of shape (N, N, 2, L). Sampling then involves rolling the (N, N) axes to the end, doing a Cholesky decomposition, and multiplying the resulting stack of matrices with the generated flm.

For real-valued fields, the 2 complex-valued Cls collapse into a pair of identical real-valued Cls (as every combination with "B" is zero). In that case, the input for sampling $N$ real-valued fields could be a stack of real-valued Cls of shape (N, N, L). This could be controlled with a reality= parameter, as you have elsewhere.

Does this make sense in terms of an s2fft-like interface? If yes, I could give it a go!