moble / spherical_functions

Python/numba package for evaluating Wigner D matrices and spin-weighted spherical harmonics
MIT License
52 stars 11 forks source link

D matrices are conjugated relative to other standards #4

Closed moble closed 9 years ago

moble commented 9 years ago

It looks like I could change mine by changing how I define the "parts" (a and b) of the quaternions. In particular, I checked that conjugating the definitions of those parts preserves the key multiplication property, so the rest of the derivation should just follow through, except that the parts would get conjugated. So this would have to change in the quaternion package. Then, I would have to change my derivation of the D matrices, and remove the corresponding conjugation in the tests.

This would make my results agree exactly with the version given on Wikipedia. At the moment, I've implemented the Wikipedia version, but when comparing values, I have to apply a conjugation to the result of the Wikipedia version.

I should also consult with Barry. I think I may have misinformed him about the signs of the exponentials in the ugly Euler form. [Maybe SphericalFunctions actually comes out with the opposite signs, but I doubt it.]

moble commented 9 years ago

I think I'm gonna stick with my version, because rotations don't get weird conjugates / transpositions with my conventions.

ablech commented 3 years ago

Hi @moble! From what I see so far, this package seems to be a quite fast and very usable implementation of Wigner matrices. I like it alot!

However, I find your description of the convention that you are following somewhat misleading. In particular, the last formular in the section about Euler angles suggests, that you are following the same convention as Wikipedia - which is not the case as you described in this issue.

It would be nice, if you would consider adding a comment about this difference, to avoid having future users getting confused as well. ;-)

phockett commented 2 years ago

I tested this a while back vs. Zare's definitions (which I had in some old Matlab routine), and also caught this difference.

Details here in case it's useful to anyone else. (I checked the 3j definitions too, and those were consistent.)

To be honest, this stuff is always a bit of a mine-field, but I think anything that makes the rotations simpler is a good thing - one just needs to pay attention when comparing to other formalisms!

moble commented 2 years ago

[Sorry for not seeing the activity here sooner.]

the last formular in the section about Euler angles suggests, that you are following the same convention as Wikipedia - which is not the case as you described in this issue.

I believe I am actually using the same convention for Euler angles of a rotation as Wikipedia, and have been at least since this issue was closed. But there are endless other problems if you want to translate that into a formula for Wigner D as a function of Euler angles. Probably most importantly, is D considered a function of the rotation or the inverse rotation?

Note that Wigner's D matrix is not inherently a function of Euler angles. It is a function of a rotation, but even conventions for which function of a rotation differ.

IIRC, in Wigner's original text, he actually defined D in terms of the inverse of the rotation. To invert a rotation in terms of Euler angles, you flip the sign on each angle and reverse the order: (alpha, beta, gamma) -> (-gamma, -beta, -alpha). The effect on D is to take a conjugate and transpose (swap m and m'). But those same changes could be made by swapping other conventions — like what m and m' mean; or which order they occur in; or whether the field described / acted on by Wigner D is a function of the body or the frame coordinates; or whether your rotation is active or passive; or whether your Euler angles are active or passive (which is different from the other active or passive), or... any of a million different subtle choices.

It's a giant mess, and that's why I say

I find it nearly impossible to decipher the meaning behind Mathematica’s documentation (or anyone else’s for that matter), so the only comparison I’m willing to make is between actual results when using different codes.

The unfortunate reality is that you have to work through everything for yourself, decide what conventions mean to you, and then decide how you interpret any code you ever use. Don't rely on words; just rely on formulas and your own interpretation.

moble commented 2 years ago

@phockett Please also take note of the message at the top of the README for this repo. A quick perusal of the ePSproc docs led me to a quick perusal of Natalense and Lucchese (1999), in which I see talk of using up to ell=30. If that's a use case for your code, you should switch to the spherical package.

There's a fundamental numerical problem with calculating Wigner's D matrices and spherical harmonics from the usual formulas so that you shouldn't trust results beyond ell=25. I had to make a whole different code to get past that problem (by using fancier recurrence formulas).

phockett commented 2 years ago

It's a giant mess, and that's why I say

Having worked with this stuff for many years, I second that! I usually go with Zare's definitions, since that's what I'm used to, and they're relatively standard in my field... but there's still a lot of issues as @moble lists (and the excellent extended notes on this :smile:), so I usually end up just testing things numerically against my reference codes to confirm results. In practice this usually means code that I'm already using, which is at least self-consistent, and means I don't have to re-derive everything again, remember the fudges and subtleties, fix the typos in the literature etc. etc.

Re: spherical package - thanks for the reminder @moble (plus the hard work on the packages!). It's on my long to-do list already (https://github.com/phockett/ePSproc/issues/35), although hasn't been a serious issue for me as yet since I'm usually not beyond L=12 in practice^. (E.g. here's a typical computational example, with lmax=8).

^ In case you're interested in the reasons here vs. Natalense and Lucchese (1999) - they're discussing expanding full wavefunctions and solving these variationally for photoionization problems, I'm using the results of that type of calculation to calculate observables, which typically have a much-truncated expansion in L for reasonable numerical thresholds. See background notes here if you're interested.