Add accurate computations for ell>29

[moble]

As explained in issue #7, the alternating sign in the sum for Wigner D elements causes huge numerical problems. It would be easy enough (now that I've changed how that sum is evaluated) to break this out into a separate C function that uses gmp to evaluate the sum. Hopefully, it would also be easy to link against gmp if gmpy2 is installed.

[raziel2646]

Hi,

Could you please explain how I can use the Wigner_D_element function for large ell?
For some reason I get an error when I evaluate the function for ell>=33:

line 113, in Wigner_D_element
    "(ell,mp,m)=({0},{1},{2})".format(twoell/2, twomp/2, twom/2)
UnboundLocalError: local variable 'twoell' referenced before assignment

[moble]

I really don't have time to work on the method I suggested above, but there is another way you could try, which I've already coded but haven't really brought into the code. I made a pull request #15 some time ago, but abandoned it because I didn't actually need the large ell values. If you'd like to try it, you'll need some git and python skills. The idea is that you should check out that pull request and use the code, so it's not trivial but I think this should do what you need.

1. Check out the package in git
2. cd into the directory that gets created
3. Run git reset --hard 0ab63aa
4. Run this code (with ell_max increased to whatever you want; this takes a while and could get ridiculous if you use a really large value):

import numpy as np
from sympy import pi
from sympy.physics.quantum.spin import Rotation
import mpmath
mpmath.mp.dps=128
ell_max=32

_binomial_coefficients = np.empty((((2*ell_max+1)*(2*ell_max+2))//2,), dtype=float)
_ladder_operator_coefficients = np.empty((((ell_max+2)*ell_max+1),), dtype=float)
_Delta = np.empty(((4*ell_max**3 + 12*ell_max**2 + 11*ell_max + 3)//3,), dtype=float)

i=0
for n in range(2*ell_max+1):
    for k in range(n+1):
        _binomial_coefficients[i]=float(mpmath.binomial(n,k))
        i+=1
print(i, _binomial_coefficients.shape)
np.save('binomial_coefficients',_binomial_coefficients)

i=0
for ell in range(ell_max+1):
    for m in range(-ell,ell+1):
        _ladder_operator_coefficients[i]=mpmath.sqrt(ell*(ell+1)-m*(m+1))
        i+=1
print(i, _ladder_operator_coefficients.shape)
np.save('ladder_operator_coefficients',_ladder_operator_coefficients)

def real(x):
    try:
        return x.real
    except AttributeError:
        return x
i=0
for ell in range(ell_max+1):
    for mp in range(-ell,ell+1):
        for m in range(-ell,ell+1):
            _Delta[i] = real(complex(Rotation.D(ell, mp, m, 0,pi/2,0).doit().evalf(n=32)))
            i += 1
print(i, _Delta.shape)
np.save('Delta',_Delta)

5. Edit the __init__.py file to use the same value of ell_max
6. Run python setup.py install
7. Use the code.

Also note that if you want high performance with really large ell values, you should probably look into other packages like spinsfast, and hack it to get Wigner D elements instead of just spin-weighted spherical harmonics (which involves calculating a complex phase).

[raziel2646]

Thanks a lot for your quick reply. Calculation of the Wigner-D elements with your package is incredibly faster than sympy, which is why I was hoping to be able to use it for very large ell (up to around 1000). I will definitely look into instructions above and spinsfast to see if I can figure it out.
Thanks again.

[moble]

It's faster than sympy because sympy uses symbolics to get exact values (which can then be evaluated to any precision you like). I doubt that you'll be patient enough to get to ell_max=1000 because evaluating the code I posted above will take forever — in particular, that last loop over ell, mp, and m. It could probably be made a lot faster, but I'm afraid I really don't have time to work on this right now.

So, I recommend using something like spinsfast, which was made to be fast and get to very large ell. You'll have to use Euler angles, and multiply by something like math.e**(1j*mp*alpha), but it will be very fast and very accurate.

[NiallMac]

Hi @moble, I'm also looking to get reduced Wigner-D elements for |m|,|m_p|<=2, up to ell~10^4. Would it be possible to get some pointers on how to do this with spinsfast? I took a look at the repo, but my untrained eye could not see which functions to use. Thanks!

[moble]

@NiallMac Sorry, I never saw your question, and just happen to be looking at this issue now, years later, when it's probably too late anyway.

For the record: I should clarify that spinsfast specializes in converting functions to and from a mode representation, and it does this on a grid. So basically, if you're looking for a 𝔇 element for an arbitrary rotation, it's not going to help at all. Plus, for the case |m|,|m_p|<=2, there's going to be a lot of wasted cycles and memory.

If you actually still care about this, I've started the successor of this package, the spherical package, which does more nearly what you want. I've only tested a few things up to ell=10^3, but I expect it would go to 10^4. It's also possible to limit the calculation to |m_p|<=2, though I haven't implemented limiting |m| at the same time. But given those limitations, you could calculate sYlm values (spin-weighted spherical harmonics), which are related to 𝔇 values by a simple multiplicative factor. If it's still relevant, take a look, and open an issue over there if you have questions.

[moble]

I won't be updating this package to compute ell>29 because, as mentioned above, this new package has been designed for that purpose.