| [l2,l3,l4,l5,l6,l1min,l1max,sixcof,ndim,ier]=rc6j(l2,l3,l4,l5,l6,l1min,l1max,sixcof,ndim,ier); |
function [l2,l3,l4,l5,l6,l1min,l1max,sixcof,ndim,ier]=rc6j(l2,l3,l4,l5,l6,l1min,l1max,sixcof,ndim,ier);
%***BEGIN PROLOGUE RC6J
%***PURPOSE Evaluate the 6j symbol h(L1) = {L1 L2 L3}
% {L4 L5 L6}
% for all allowed values of L1, the other parameters
% being held fixed.
%***LIBRARY SLATEC
%***CATEGORY C19
%***TYPE SINGLE PRECISION (RC6J-S, DRC6J-D)
%***KEYWORDS 6J COEFFICIENTS, 6J SYMBOLS, CLEBSCH-GORDAN COEFFICIENTS,
% RACAH COEFFICIENTS, VECTOR ADDITION COEFFICIENTS,
% WIGNER COEFFICIENTS
%***AUTHOR Gordon, R. G., Harvard University
% Schulten, K., Max Planck Institute
%***DESCRIPTION
%
% *Usage:
%
% REAL L2, L3, L4, L5, L6, L1MIN, L1MAX, SIXCOF(NDIM)
% INTEGER NDIM, IER
%
% CALL RC6J(L2, L3, L4, L5, L6, L1MIN, L1MAX, SIXCOF, NDIM, IER)
%
% *Arguments:
%
% L2 :IN Parameter in 6j symbol.
%
% L3 :IN Parameter in 6j symbol.
%
% L4 :IN Parameter in 6j symbol.
%
% L5 :IN Parameter in 6j symbol.
%
% L6 :IN Parameter in 6j symbol.
%
% L1MIN :OUT Smallest allowable L1 in 6j symbol.
%
% L1MAX :OUT Largest allowable L1 in 6j symbol.
%
% SIXCOF :OUT Set of 6j coefficients generated by evaluating the
% 6j symbol for all allowed values of L1. SIXCOF(I)
% will contain h(L1MIN+I-1), I=1,2,...,L1MAX-L1MIN+1.
%
% NDIM :IN Declared length of SIXCOF in calling program.
%
% IER :OUT Error flag.
% IER=0 No errors.
% IER=1 L2+L3+L5+L6 or L4+L2+L6 not an integer.
% IER=2 L4, L2, L6 triangular condition not satisfied.
% IER=3 L4, L5, L3 triangular condition not satisfied.
% IER=4 L1MAX-L1MIN not an integer.
% IER=5 L1MAX less than L1MIN.
% IER=6 NDIM less than L1MAX-L1MIN+1.
%
% *Description:
%
% The definition and properties of 6j symbols can be found, for
% example, in Appendix C of Volume II of A. Messiah. Although the
% parameters of the vector addition coefficients satisfy certain
% conventional restrictions, the restriction that they be non-negative
% integers or non-negative integers plus 1/2 is not imposed on input
% to this subroutine. The restrictions imposed are
% 1. L2+L3+L5+L6 and L2+L4+L6 must be integers;
% 2. ABS(L2-L4).LE.L6.LE.L2+L4 must be satisfied;
% 3. ABS(L4-L5).LE.L3.LE.L4+L5 must be satisfied;
% 4. L1MAX-L1MIN must be a non-negative integer, where
% L1MAX=MIN(L2+L3,L5+L6) and L1MIN=MAX(ABS(L2-L3),ABS(L5-L6)).
% If all the conventional restrictions are satisfied, then these
% restrictions are met. Conversely, if input to this subroutine meets
% all of these restrictions and the conventional restriction stated
% above, then all the conventional restrictions are satisfied.
%
% The user should be cautious in using input parameters that do
% not satisfy the conventional restrictions. For example, the
% the subroutine produces values of
% h(L1) = { L1 2/3 1 }
% {2/3 2/3 2/3}
% for L1=1/3 and 4/3 but none of the symmetry properties of the 6j
% symbol, set forth on pages 1063 and 1064 of Messiah, is satisfied.
%
% The subroutine generates h(L1MIN), h(L1MIN+1), ..., h(L1MAX)
% where L1MIN and L1MAX are defined above. The sequence h(L1) is
% generated by a three-term recurrence algorithm with scaling to
% control overflow. Both backward and forward recurrence are used to
% maintain numerical stability. The two recurrence sequences are
% matched at an interior point and are normalized from the unitary
% property of 6j coefficients and Wigner's phase convention.
%
% The algorithm is suited to applications in which large quantum
% numbers arise, such as in molecular dynamics.
%
%***REFERENCES 1. Messiah, Albert., Quantum Mechanics, Volume II,
% North-Holland Publishing Company, 1963.
% 2. Schulten, Klaus and Gordon, Roy G., Exact recursive
% evaluation of 3j and 6j coefficients for quantum-
% mechanical coupling of angular momenta, J Math
% Phys, v 16, no. 10, October 1975, pp. 1961-1970.
% 3. Schulten, Klaus and Gordon, Roy G., Semiclassical
% approximations to 3j and 6j coefficients for
% quantum-mechanical coupling of angular momenta,
% J Math Phys, v 16, no. 10, October 1975,
% pp. 1971-1988.
% 4. Schulten, Klaus and Gordon, Roy G., Recursive
% evaluation of 3j and 6j coefficients, Computer
% Phys Comm, v 11, 1976, pp. 269-278.
%***ROUTINES CALLED R1MACH, XERMSG
%***REVISION HISTORY (YYMMDD)
% 750101 DATE WRITTEN
% 880515 SLATEC prologue added by G. C. Nielson, NBS; parameters
% HUGE and TINY revised to depend on R1MACH.
% 891229 Prologue description rewritten; other prologue sections
% revised; LMATCH (location of match point for recurrences)
% removed from argument list; argument IER changed to serve
% only as an error flag (previously, in cases without error,
% it returned the number of scalings); number of error codes
% increased to provide more precise error information;
% program comments revised; SLATEC error handler calls
% introduced to enable printing of error messages to meet
% SLATEC standards. These changes were done by D. W. Lozier,
% M. A. McClain and J. M. Smith of the National Institute
% of Standards and Technology, formerly NBS.
% 910415 Mixed type expressions eliminated; variable C1 initialized;
% description of SIXCOF expanded. These changes were done by
% D. W. Lozier.
%***end PROLOGUE RC6J
%
persistent a1 a1s a2 a2s c1 c1old c2 cnorm denom dv eps firstCall gt hugemlv i indexmlv l1 lstep n newfac nfin nfinp1 nfinp2 nfinp3 nlim nstep2 oldfac one ratio sign1 sign2 srhuge srtiny sum1 sum2 sumbac sumfor sumuni three thresh tinymlv two x x1 x2 x3 y y1 y2 y3 zero ; if isempty(firstCall),firstCall=1;end;
if isempty(gt), gt=zeros(1,2); end;
%
if isempty(i), i=0; end;
if isempty(indexmlv), indexmlv=0; end;
if isempty(lstep), lstep=0; end;
if isempty(n), n=0; end;
if isempty(nfin), nfin=0; end;
if isempty(nfinp1), nfinp1=0; end;
if isempty(nfinp2), nfinp2=0; end;
if isempty(nfinp3), nfinp3=0; end;
if isempty(nlim), nlim=0; end;
if isempty(nstep2), nstep2=0; end;
if isempty(a1), a1=0; end;
if isempty(a1s), a1s=0; end;
if isempty(a2), a2=0; end;
if isempty(a2s), a2s=0; end;
if isempty(c1), c1=0; end;
if isempty(c1old), c1old=0; end;
if isempty(c2), c2=0; end;
if isempty(cnorm), cnorm=0; end;
if isempty(denom), denom=0; end;
if isempty(dv), dv=0; end;
if isempty(eps), eps=0; end;
if isempty(hugemlv), hugemlv=0; end;
if isempty(l1), l1=0; end;
if isempty(newfac), newfac=0; end;
if isempty(oldfac), oldfac=0; end;
if isempty(one), one=0; end;
if isempty(ratio), ratio=0; end;
if isempty(sign1), sign1=0; end;
if isempty(sign2), sign2=0; end;
if isempty(srhuge), srhuge=0; end;
if isempty(srtiny), srtiny=0; end;
if isempty(sum1), sum1=0; end;
if isempty(sum2), sum2=0; end;
if isempty(sumbac), sumbac=0; end;
if isempty(sumfor), sumfor=0; end;
if isempty(sumuni), sumuni=0; end;
if isempty(three), three=0; end;
if isempty(thresh), thresh=0; end;
if isempty(tinymlv), tinymlv=0; end;
if isempty(two), two=0; end;
if isempty(x), x=0; end;
if isempty(x1), x1=0; end;
if isempty(x2), x2=0; end;
if isempty(x3), x3=0; end;
if isempty(y), y=0; end;
if isempty(y1), y1=0; end;
if isempty(y2), y2=0; end;
if isempty(y3), y3=0; end;
if isempty(zero), zero=0; end;
%
if firstCall, zero =[0.0]; end;
if firstCall, eps =[0.01]; end;
if firstCall, one =[1.0]; end;
if firstCall, two =[2.0]; end;
if firstCall, three=[3.0]; end;
firstCall=0;
%
%***FIRST EXECUTABLE STATEMENT RC6J
ier = 0;
% HUGE is the square root of one twentieth of the largest floating
% point number, approximately.
hugemlv = sqrt(r1mach(2)./20.0);
srhuge = sqrt(hugemlv);
tinymlv = 1.0./hugemlv;
srtiny = 1.0./srhuge;
%
% LMATCH = ZERO
%
% Check error conditions 1, 2, and 3.
gt(:)=0;
if((rem(l2+l3+l5+l6+eps,one)>=eps+eps) ||(rem(l4+l2+l6+eps,one)>=eps+eps) )
ier = 1;
[dumvar1,dumvar2,dumvar3,ier]=xermsg('SLATEC','RC6J',['L2+L3+L5+L6 or L4+L2+L6 not ','integer.'],ier,1);
elseif((l4+l2-l6<zero) ||(l4-l2+l6<zero) ||(-l4+l2+l6<zero) ) ;
ier = 2;
[dumvar1,dumvar2,dumvar3,ier]=xermsg('SLATEC','RC6J',['L4, L2, L6 triangular ','condition not satisfied.'],ier,1);
elseif((l4-l5+l3<zero) ||(l4+l5-l3<zero) ||(-l4+l5+l3<zero) ) ;
ier = 3;
[dumvar1,dumvar2,dumvar3,ier]=xermsg('SLATEC','RC6J',['L4, L5, L3 triangular ','condition not satisfied.'],ier,1);
else;
%
% Limits for L1
%
l1min = max(abs(l2-l3),abs(l5-l6));
l1max = min(l2+l3,l5+l6);
%
% Check error condition 4.
if( rem(l1max-l1min+eps,one)>=eps+eps )
ier = 4;
[dumvar1,dumvar2,dumvar3,ier]=xermsg('SLATEC','RC6J','L1MAX-L1MIN not integer.',ier,1);
elseif( l1min<l1max-eps ) ;
%
%
% This is reached in case that L1 can take more than one value.
%
% LSCALE = 0
nfin = fix(fix(l1max-l1min+one+eps));
if( ndim<nfin )
%
% Check error condition 6.
ier = 6;
[dumvar1,dumvar2,dumvar3,ier]=xermsg('SLATEC','RC6J',['Dimension of result array for 6j ','coefficients too small.'],ier,1);
else;
%
%
% Start of forward recursion
%
l1 = l1min;
newfac = 0.0;
c1 = 0.0;
sixcof(1) = srtiny;
sum1 =(l1+l1+one).*tinymlv;
%
lstep = 1;
while( true );
gt(:)=0;
lstep = fix(lstep + 1);
l1 = l1 + one;
%
oldfac = newfac;
a1 =(l1+l2+l3+one).*(l1-l2+l3).*(l1+l2-l3).*(-l1+l2+l3+one);
a2 =(l1+l5+l6+one).*(l1-l5+l6).*(l1+l5-l6).*(-l1+l5+l6+one);
newfac = sqrt(a1.*a2);
%
if( l1<one+eps )
%
% If L1 = 1, (L1 - 1) has to be factored out of DV, hence
%
c1 = -two.*(l2.*(l2+one)+l5.*(l5+one)-l4.*(l4+one))./newfac;
else;
%
dv = two.*(l2.*(l2+one).*l5.*(l5+one)+l3.*(l3+one).*l6.*(l6+one)-l1.*(l1-one).*l4.*(l4+one))-(l2.*(l2+one)+l3.*(l3+one)-l1.*(l1-one)).*(l5.*(l5+one)+l6.*(l6+one)-l1.*(l1-one));
%
denom =(l1-one).*newfac;
%
%
if( lstep>2 )
c1old = abs(c1);
end;
c1 = -(l1+l1-one).*dv./denom;
end;
%
if( lstep>2 )
%
%
c2 = -l1.*oldfac./denom;
%
% Recursion to the next 6j coefficient X
%
x = c1.*sixcof(lstep-1) + c2.*sixcof(lstep-2);
sixcof(lstep) = x;
%
sumfor = sum1;
sum1 = sum1 +(l1+l1+one).*x.*x;
if( lstep==nfin )
gt(1)=1;
break;
end;
%
% See if last unnormalized 6j coefficient exceeds SRHUGE
%
if( abs(x)>=srhuge )
%
% This is reached if last 6j coefficient larger than SRHUGE,
% so that the recursion series SIXCOF(1), ... ,SIXCOF(LSTEP)
% has to be rescaled to prevent overflow
%
% LSCALE = LSCALE + 1
for i = 1 : lstep;
if( abs(sixcof(i))<srtiny )
sixcof(i) = zero;
end;
sixcof(i) = sixcof(i)./srhuge;
end; i = fix(lstep+1);
sum1 = sum1./hugemlv;
sumfor = sumfor./hugemlv;
x = x./srhuge;
end;
%
%
% As long as the coefficient ABS(C1) is decreasing, the recursion
% proceeds towards increasing 6j values and, hence, is numerically
% stable. Once an increase of ABS(C1) is detected, the recursion
% direction is reversed.
%
if( c1old<=abs(c1) )
gt(1)=1;
break;
end;
else;
%
% If L1 = L1MIN + 1, the third term in recursion equation vanishes
%
x = srtiny.*c1;
sixcof(2) = x;
sum1 = sum1 + tinymlv.*(l1+l1+one).*c1.*c1;
%
if( lstep==nfin )
break;
end;
end;
end;
%
if(gt(1)==0)
sumuni = sum1;
gt(2)=1;
end;
%
%
% Keep three 6j coefficients around LMATCH for comparison later
% with backward recursion.
%
% LMATCH = L1 - 1
if(gt(2)==0)
x1 = x;
x2 = sixcof(lstep-1);
x3 = sixcof(lstep-2);
%
%
%
% Starting backward recursion from L1MAX taking NSTEP2 steps, so
% that forward and backward recursion overlap at the three points
% L1 = LMATCH+1, LMATCH, LMATCH-1.
%
nfinp1 = fix(nfin + 1);
nfinp2 = fix(nfin + 2);
nfinp3 = fix(nfin + 3);
nstep2 = fix(nfin - lstep + 3);
l1 = l1max;
%
sixcof(nfin) = srtiny;
sum2 =(l1+l1+one).*tinymlv;
%
%
l1 = l1 + two;
lstep = 1;
while( true );
lstep = fix(lstep + 1);
l1 = l1 - one;
%
oldfac = newfac;
a1s =(l1+l2+l3).*(l1-l2+l3-one).*(l1+l2-l3-one).*(-l1+l2+l3+two);
a2s =(l1+l5+l6).*(l1-l5+l6-one).*(l1+l5-l6-one).*(-l1+l5+l6+two);
newfac = sqrt(a1s.*a2s);
%
dv = two.*(l2.*(l2+one).*l5.*(l5+one)+l3.*(l3+one).*l6.*(l6+one)-l1.*(l1-one).*l4.*(l4+one))-(l2.*(l2+one)+l3.*(l3+one)-l1.*(l1-one)).*(l5.*(l5+one)+l6.*(l6+one)-l1.*(l1-one));
%
denom = l1.*newfac;
c1 = -(l1+l1-one).*dv./denom;
if( lstep>2 )
%
%
c2 = -(l1-one).*oldfac./denom;
%
% Recursion to the next 6j coefficient Y
%
y = c1.*sixcof(nfinp2-lstep) + c2.*sixcof(nfinp3-lstep);
if( lstep==nstep2 )
break;
end;
sixcof(nfinp1-lstep) = y;
sumbac = sum2;
sum2 = sum2 +(l1+l1-three).*y.*y;
%
% See if last unnormalized 6j coefficient exceeds SRHUGE
%
if( abs(y)>=srhuge )
%
% This is reached if last 6j coefficient larger than SRHUGE,
% so that the recursion series SIXCOF(NFIN), ... ,SIXCOF(NFIN-LSTEP+1)
% has to be rescaled to prevent overflow
%
% LSCALE = LSCALE + 1
for i = 1 : lstep;
indexmlv = nfin - i + 1;
if( abs(sixcof(indexmlv))<srtiny )
sixcof(indexmlv) = zero;
end;
sixcof(indexmlv) = sixcof(indexmlv)./srhuge;
end; i = fix(lstep+1);
sumbac = sumbac./hugemlv;
%
sum2 = sum2./hugemlv;
end;
else;
%
% If L1 = L1MAX + 1 the third term in the recursion equation vanishes
%
y = srtiny.*c1;
sixcof(nfin-1) = y;
if( lstep==nstep2 )
break;
end;
sumbac = sum2;
sum2 = sum2 +(l1+l1-three).*c1.*c1.*tinymlv;
end;
end;
%
%
% The forward recursion 6j coefficients X1, X2, X3 are to be matched
% with the corresponding backward recursion values Y1, Y2, Y3.
%
y3 = y;
y2 = sixcof(nfinp2-lstep);
y1 = sixcof(nfinp3-lstep);
%
%
% Determine now RATIO such that YI = RATIO * XI (I=1,2,3) holds
% with minimal error.
%
ratio =(x1.*y1+x2.*y2+x3.*y3)./(x1.*x1+x2.*x2+x3.*x3);
nlim = fix(nfin - nstep2 + 1);
%
if( abs(ratio)<one )
%
nlim = fix(nlim + 1);
ratio = one./ratio;
for n = nlim : nfin;
sixcof(n) = ratio.*sixcof(n);
end; n = fix(nfin+1);
sumuni = sumfor + ratio.*ratio.*sumbac;
else;
%
for n = 1 : nlim;
sixcof(n) = ratio.*sixcof(n);
end; n = fix(nlim+1);
sumuni = ratio.*ratio.*sumfor + sumbac;
end;
end;
%
%
% Normalize 6j coefficients
%
cnorm = one./sqrt((l4+l4+one).*sumuni);
%
% Sign convention for last 6j coefficient determines overall phase
%
sign1 = (abs(one).*sign(sixcof(nfin)));
sign2 =(-one).^fix(l2+l3+l5+l6+eps);
if( sign1.*sign2<=0 )
cnorm = -cnorm;
end;
%
if( abs(cnorm)<one )
%
thresh = tinymlv./abs(cnorm);
for n = 1 : nfin;
if( abs(sixcof(n))<thresh )
sixcof(n) = zero;
end;
sixcof(n) = cnorm.*sixcof(n);
end; n = fix(nfin+1);
else;
%
for n = 1 : nfin;
sixcof(n) = cnorm.*sixcof(n);
end; n = fix(nfin+1);
end;
end;
elseif( l1min<l1max+eps ) ;
%
%
% This is reached in case that L1 can take only one value
%
% LSCALE = 0
sixcof(1) =(-one).^fix(l2+l3+l5+l6+eps)./sqrt((l1min+l1min+one).*(l4+l4+one));
else;
%
% Check error condition 5.
ier = 5;
[dumvar1,dumvar2,dumvar3,ier]=xermsg('SLATEC','RC6J','L1MIN greater than L1MAX.',ier,1);
end;
end;
%
end %subroutine rc6j
%DECK RC
|
|
