On Thu, Jun 16, 2011 at 8:34 PM, Ondrej Certik <[email protected]> wrote: > On Thu, Jun 16, 2011 at 6:38 PM, Ondrej Certik <[email protected]> > wrote: >> On Thu, Jun 16, 2011 at 12:04 AM, Sean Vig <[email protected]> wrote: >>> Hi all, >>> In working on some stuff with spin states, I ran into some problems with the >>> current implementation of the Wigner small-d matrix, Rotation.d in >>> sympy.physics.quantum.spin. I had written methods to change bases using the >>> Wigner D-function [0] and in testing decided to try >>>>>>qapply(JzBra(1,1)*JzKet(1,1).rewrite('Jx')) >>> which should be 1, as the spin ket is rewritten in the Jx basis and then >>> back to the Jz basis to apply the innerproduct, but I have that it gives 0. >>> I traced this back to a bug in the Rotation.d function, which currently has >>> an open issue [1]. For the Wigner D-function and the small d-matrix, the >>> conventions laid out in Varshalovich "Quantum Theory of Angular Momentum". >>> It seems the d-matrix is fine for positive values of angle argument, but >>> does not obey the symmetry d(j,m,mp,-beta)=(-1)**(m-mp)*d(j,m,mp,beta) and >>> does not agree with the tables in Varshalovich. Those terms that fail for >>> the j=1 case are in an XFAIL test in my branch. What is odd is that when I >> >> >> I think there is a bug in the code. See my comment here: >> >> http://code.google.com/p/sympy/issues/detail?id=2423#c3 >> >> The correct results for general beta are: >> >> d(1, 0, 1, beta) = sin(beta)/sqrt(2) >> d(1, 1, 0, beta) = -sin(beta)/sqrt(2) >> >> However, sympy gives: >> >>>>> Rotation.d(1,0,1,beta) >> ⎽⎽⎽ >> ╲╱ 2 ⋅(2⋅cos(β) + 2) >> ──────────────────── >> 4 >>>>> Rotation.d(1,1,0,beta) >> ⎽⎽⎽ >> -╲╱ 2 ⋅(cos(β) + 1) >> ─────────────────── >> 2 >> >> Which is wrong (it looks quite ok, that it changes sign, as it should, >> but something is wrong with the cos(beta) thing). The sympy code >> implements the "Eq. 7 in Section 4.3.2 of Varshalovich." >> >> >>> ran the equation used to define the d-matrix through Mathematica, I got >>> results that agreed with the sympy output, so the problem may be in the >>> equation and not a bug in the code. If anyone could take a look at that, I'd >>> appreciate it. >> >> Which *exact* equation did you run through Mathematica? Eq. 7 in >> section 4.3.2? What *exactly* did you get? Did you get an expression >> involving cos(beta), just like sympy above? Can you paste here the >> Mathematica code? I'll run it with my Mathematica, to verify, that we >> didn't make a mistake. >> >> >> Now we just need to systematically look at it, and nail it down. We >> need to get the correct expressions: >> >> d(1, 0, 1, beta) = sin(beta)/sqrt(2) >> d(1, 1, 0, beta) = -sin(beta)/sqrt(2) >> >> one way or the other, for general beta. Then things will start to work. >> >> Sean, let me know if you have any questions to the above. >> >> Once we fix this, we'll move on to the other problems you raised. >> >> I am CCing Brian, who implemented that code in SymPy. However, we >> should be able to fix this ourselves Sean ---- all we need to do is to >> take the eq. (7), and see what expression we get for J=1, M=1, M'=0, >> just put it there by hand (don't use mathematica) and see what you >> get. >> >> Post here your results, and I'll verify them with my independent >> calculation and we'll nail it down. > > Ok, I think that I have nailed it down. There are actually several problems: > > 1) First of all, this is the range for beta: > > 0 <= beta <= pi > > see Varshalovich, page 74. > > 2) See the attached screenshot of my calculation, which calculates > d(j, 0, 1, beta) from the equation (7) on page 75, and shows, that it > is equal to > sin(beta)/sqrt(2), consistent with Varshalovich result in the table 4.4. > > > As such, from 2) it follows, that the sympy result (see my previous > email) is *wrong*, as can be checked by substituting beta = pi/3 and > checking against sin(beta)/sqrt(2). For beta=pi/2, it happens to be > equal, but that is pure accident. > > From 1) it follows, that you can't call it for beta=-pi/2, you need to > adjust alpha and gamma instead.
Actually, I wasn't very clear here. You use the formula (7) to calculate the wigner d function for *any* J, M, M', for 0<=beta<=pi. In particular, you get the formulas: d(1, 1, 0, beta) = -sin(beta)/sqrt(2) d(1, 0, 1, beta) = +sin(beta)/sqrt(2) And I am stressing here that this is *only* valid for 0 <= beta <= pi. You can equivalently write the sin(beta) using cos as sin(beta) = sqrt(1-cos**2(beta)), which is valid in this whole domain. You actually only get expressions involving cos(beta) from (7), but since we are on this restricted domain, it doesn't matter. Now when you want to calculate the wigner d function for parameters outside of this domain, you need to use the symmetries, see the section 4.4. Example: calculate d(1, 1, 0, -pi/2), as you needed above. Then we use this symmetry: d(j, m, mp, -beta) = d(j, mp, m, beta), and write: d(1, 1, 0, -pi/2) = d(1, 0, 1, pi/2) = sin(pi/2)/sqrt(2) = 1/sqrt(2) and that's it. All should be clear now. Ondrej -- You received this message because you are subscribed to the Google Groups "sympy" group. To post to this group, send email to [email protected]. To unsubscribe from this group, send email to [email protected]. For more options, visit this group at http://groups.google.com/group/sympy?hl=en.
