Permalink
Switch branches/tags
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
278 lines (266 sloc) 10.2 KB
function W=wigner3j(l1,l2,l3,m1,m2,m3,method,recon)
% DO NOT USE THIS FUNCTION: SEE WIGNER0J, WIGNER3JM, THREEJ.
% MOST ANNOYINGLY, THE LAST ORDER SIGN IS UNINTUITIVE. ALSO, IT'S A
% MIXTURE OF WIGNER AND GAUNT COEFFICIENTS AND NORMALIZATIONS. SEE DEMO6.
% Used only for initialization of WIGNERJM/WIGNER0J.
%
% W=wigner3j(l1,l2,l3,m1,m2,m3,method,recon)
%
% Calculates absolute values of Wigner 3j coefficients in the context of
% UN-normalized associated Legendre functions. One half the inner product
% over the entire interval of three standard, unnormalized, associated
% Legendre polynomials of degree l1,l2,l3 and positive order m1,m2,m3
% is the square of this Wigner 3j symbol, Dahlen & Tromp Eq. (C.226).
%
% INPUT:
%
% l1,l2,l3 Top row of the Wigner3j symbol; positive angular degree
% m1,m2,m3 Bottom row of the Wigner3j symbol
% method 'automatic' (default), or force:
% 'gl' by Gauss-Legendre integration
% 'recu' by recursion formula
% Note that, depending on what's possible, you may not
% always get what you asked for, it defaults back
% recon 1 Compute, maybe save to database [default]
% 0 Allow database reading
%
% OUTPUT:
%
% W The result, perhaps with all the recursions as a vector.
%
% EXAMPLE:
%
% wigner3j('demo1') % Compare against LEGENDREPRODINT
% wigner3j('demo2') % Compare the two calculation methods
% wigner3j('demo3') % Test various normalizations
% wigner3j('demo4') % Test against an exact arithmetic calculator
% wigner3j('demo5') % Test the special formulas against WIGNER3JM
% wigner3j('demo6') % Response to a user query
%
% See also: WIGNER0J, WIGNER3JM
%
% Last modified by fjsimons-at-alum.mit.edu, 05/31/2011
defval('m1',0)
defval('m2',0)
defval('m3',0)
defval('method','automatic')
defval('recon',1)
% What's coming behind here more properly fits in WIGNER3JM and uses a
% different normalization than what is being used behind
if m1<0 || m2<0 || m3<0
% Only special cases are allowed, with the Condon-Shortley phase
% Top row is of the form (l+a b l) with a=0,1,2, and b=1,2
% Bottom row is of the form (-m 0 m)
if l2==1 && m2==0 && m1==-m3 && l1==l3+1
l=l3; m=abs(m1);
W=(-1)^(l+m+1)*sqrt(2*((l+1)^2-m^2)/((2*l+3)*(2*l+2)*(2*l+1)));
disp('WIGNER3J by special formula 1')
return
end
if l2==2 && m2==0 && m1==-m3 && l1==l3+2
l=l3; m=abs(m1);
W=(-1)^(l+m)*sqrt(6*((l+1)^2-m^2)*((l+2)^2-m^2)/...
((2*l+5)*(2*l+4)*(2*l+3)*(2*l+2)*(2*l+1)));
disp('WIGNER3J by special formula 2')
return
end
if l2==2 && m2==0 && m1==-m3 && l1==l3+1
l=l3; m=abs(m1);
W=(-1)^(l+m+1)*2*m*sqrt(6*((l+1)^2-m^2)/...
((2*l+4)*(2*l+3)*(2*l+2)*(2*l+1)*2*l));
disp('WIGNER3J by special formula 3')
return
end
if l2==2 && m2==0 && m1==-m3 && l1==l3
l=l3; m=abs(m1);
W=(-1)^(l+m)*2*(3*m^2-l*(l+1))/...
sqrt((2*l-1)*(2*l)*(2*l+1)*(2*l+2)*(2*l+3));
disp('WIGNER3J by special formula 4')
return
else
error('Better use WIGNER3JM, WIGNERCYCLE, and THREEJ ')
end
end
% What's coming behind here has a very different normalization and no
% longer properly fits with the shortcuts inserted above. Now it treats
% the input orders as those of the triple Legendre product and constructs
% a Wigner-like symbol for it. Combined with the different normalization,
% this is not the best approach here, though it's still useful to check.
if ~isstr(l1)
fnpl=sprintf('%s/W3J-%i-%i-%i-%i-%i-%i.mat',...
fullfile(getenv('IFILES'),'WIGNER'),l1,l2,l3,m1,m2,m3);
if exist(fnpl,'file')==2 && recon==0
load(fnpl)
%disp(sprintf('%s loaded by WIGNER3J',fnpl))
else
% Negative m's are allowed unless the result is computed by integration
if abs(m1)>l1 || abs(m2)>l2 || abs(m3)>l3
error('Mismatch of angular degrees and orders')
end
% Conditions of Dahlen and Tromp Eq. (C.186) % BUT NO NEGATIVE M
if ~triangle(l1,l2,l3) || m1~=(m2+m3)
W=0;
return
end
% Conditions of Dahlen and Tromp Eq. (C.219)
if m1==0 && m2==0 && m3==0 && mod(l1+l2+l3,2)
W=0;
return
end
% Dahlen and Tromp Eqs. (B.49) and (C.226)
% Effectively it is now a double, not a triple product
if l3==0 && m1==m2 && ~strcmp(method,'recu') && ~strcmp(method,'gl')
% The triangle equality with l3=0 implies l1==l2
% The initial check on m implies that m3==0
if (abs(l1)+abs(m1))<21
disp('WIGNER3J by normalization formula Dahlen and Tromp B.49')
W=sqrt(factorial(l1+m1)/factorial(l1-m1)/(2*l1+1));
return
else
% If factorial fails need to use Gauss-Legendre
W=wigner3j(l1,l2,l3,m1,m2,m3,'gl');
return
end
end
% Via Gauss-Legendre integration, forced or by choice
if (l1~=l2 || m1~=0 || m2~=0 || m3~=0 || strcmp(method,'gl')) ...
&& ~strcmp(method,'recu')
% If the m is odd then the associated is not a polynomial
% and this is not the best way; would need to add new recursion
% relations to cover this case.
integrand=inline(sprintf(...
['rindeks(legendre(%i,x),%i).*',...
'rindeks(legendre(%i,x),%i).*',...
'rindeks(legendre(%i,x),%i)'],...
l1,m1+1,l2,m2+1,l3,m3+1));
% Sign information is lost, however.
% I suppose you can put this back in from Luscombe
% sign of last term is (-1)^(l2-l3+m2+m3), see WIGNER0J
% No wait, it is DT C220
% Dahlen and Tromp Eq. (C.226)
W=sqrt(gausslegendre([-1 1],integrand,l1+l2+l3)/2);
disp('WIGNER3J by Gauss-Legendre integration of triple product')
return
end
% All-zero bottom row; Dahlen and Tromp Eq. (C.220); sign positive
% Not by recursion
if m1==0 && m2==0 && m3==0 && (l1+l2+l3)<20 && ~strcmp(method,'recu')
S=(l1+l2+l3);
disp('WIGNER3J by summation formula Dahlen and Tromp C.220')
% PUT THE SIGN IN THERE DUDE
W=sqrt(1/factorial(S+1)*...
factorial(S-2*l1)*factorial(S-2*l2)*factorial(S-2*l3))*...
factorial(S/2)/...
factorial(S/2-l1)/factorial(S/2-l2)/factorial(S/2-l3);
return
end
% All-zero bottom row and l1==l2;
% Recursion http://functions.wolfram.com/07.39.17.0013.01
if m1==0 && m2==0 && m3==0 && l1==l2
% By upward recursion and saving all the terms
W=repmat(NaN,1,l3+1);
% Dahlen and Tromp Eq. (B.49)
W(1)=1/sqrt(2*l1+1); % Zeroth term
for j=2:2:l3
% All odd terms zero from triangle rule with l1==l2
W(j)=0;
% Keep sign positive
W(j+1)=(j-1)/j*...
sqrt((4*l1^2+4*l1-j^2+2*j)/(4*l1^2+4*l1-j^2+1))...
*W(j-1);
end
disp('WIGNER3J by recursion formula Mathematica')
end
end
% If you don't have anything by now, start all over with the defaults
defval('W',wigner3j(l1,l2,l3,m1,m2,m3))
elseif strcmp(l1,'demo1')
j=0; m=0;
for l=0:100
d(l+1)=wigner3j(l,l,j,m,m,m)^2-legendreprodint(l,m,l,m,-1,'gl')/2;
end
clf
plot(0:100,abs(d)/ep,'LineW',2)
xlabel('l')
ylabel(sprintf('W_{ll}^0 difference %s','\times'))
title('Difference between WIGNER3J and direct GL integration');
grid on
elseif strcmp(l1,'demo2')
m=0; l=round(rand*100);
for j=0:100
D=wigner3j(l,l,j,m,m,m);
d(j+1)=D(end)-wigner3j(l,l,j,m,m,m,'gl');
clear D
end
clf
plot(0:100,abs(d)/eps,'LineW',2)
xlabel('j')
ylabel(sprintf('W_{%i%i}^j difference %s eps',l,l,'\times'))
title('Difference between WIGNER3J using recursion and GL integration');
grid on
elseif strcmp(l1,'demo3')
l=round(100*rand(1));
disp(sprintf('Test normalization at degree %i',l))
difer(wigner3j(l,l,0,0,0,0)^2*(2*l+1)-1,[],[],NaN)
l=round(100*rand(1)); j=round(100*rand(1));
disp(sprintf('Test normalization at degrees %i and %i',l,j))
a=0; for i=0:l+j; a=a+wigner3j(l,j,i)^2*(2*i+1); end
difer(a-1,[],[],NaN)
elseif strcmp(l1,'demo4')
% Compare with the algorithm by Stone and Wood, which differs
% at most by the sign, e.g.
disp('See http://www-stone.ch.cam.ac.uk/cgi-bin/wigner.cgi')
difer(wigner3j(1,9,10)-sqrt(10/399),[],[],NaN)
difer(wigner3j(5,0,5,0,0,0)-sqrt(1/11),[],[],NaN)
elseif strcmp(l1,'demo5')
l=round(100*rand(1));
disp(sprintf('Compare with WIGNER3JM at degree %i',l))
difer(wigner3j(l+1,1,l,-1,0,1)-indeks(...
wigner3jm(l+1,1,l,-1,0,1),'end'),[],[],NaN)
l=round(100*rand(1));
disp(sprintf('Compare with WIGNER3JM at degree %i',l))
difer(wigner3j(l+2,2,l,-1,0,1)-indeks(...
wigner3jm(l+2,2,l,-1,0,1),'end'),[],[],NaN)
l=round(100*rand(1));
disp(sprintf('Compare with WIGNER3JM at degree %i',l))
difer(wigner3j(l+1,2,l,-1,0,1)-indeks(...
wigner3jm(l+1,2,l,-1,0,1),'end'),[],[],NaN)
l=round(100*rand(1));
disp(sprintf('Compare with WIGNER3JM at degree %i',l))
difer(wigner3j(l+0,2,l,-1,0,1)-indeks(...
wigner3jm(l+0,2,l,-1,0,1),'end'),[],[],NaN)
elseif strcmp(l1,'demo6')
% If you want something that is proportional to the triple product of
% spherical harmonics Y_{53},Y_{00},Y_{53}, you'll want to try
%l=5; m=3;
% But now let's do random values
maxl=48;
l=round(maxl*rand(1)); m=round(l*rand(1));
disp(sprintf(...
'Comparing with WIGNER3JM at degree %i and order %i %i',l,m))
W1=wigner3j(l,0,l,m,0,m);
% or else this which is just by a different method
W2=wigner3j(l,l,0,m,m,0);
% but then you also expect this to be proportional to
W3=indeks(wigner3jm(l,0,l,m,0,-m),'end');
% and you should be able to also use
W4=threej(l,0,l,m,0,-m);
% and the relation between the first and the last two, which pairs
% are identical, is due to an unfortunate mixture of normalization
% conventions and a deliberate confusion of Wigner and Gaunt
% coefficients. Bottom line is that this is why WIGNER3J is obsolete.
% Note that also the sign in WIGNER3J is not properly calibrated
format long
abs([W1/sqrt(factorial(l+m)/factorial(l-m)) ;
W2/sqrt(factorial(l+m)/factorial(l-m)) ; W3 ; W4])
% But the conclusion is that the various procedures are completely
% attuned to one another.
% Lastly, check your integrals over the sphere, this should be one
W1=indeks(wigner0j(l,0,l),'end')...
*indeks(wigner3jm(l,0,l,m,0,-m),'end');
W2=indeks(wigner3jm(l,0,l,0,0,0),'end')...
*indeks(wigner3jm(l,0,l,m,0,-m),'end');
W3=zeroj(l,0,l)*threej(l,0,l,m,0,-m);
(-1)^m*[W1 ; W2 ; W3]*sqrt((2*l+1)*(2*0+1)*(2*l+1))
format short
end