| [l1,l2,l3,m1,m2min,m2max,thrcof,ndim,ier]=drc3jm(l1,l2,l3,m1,m2min,m2max,thrcof,ndim,ier); |
function [l1,l2,l3,m1,m2min,m2max,thrcof,ndim,ier]=drc3jm(l1,l2,l3,m1,m2min,m2max,thrcof,ndim,ier);
%***BEGIN PROLOGUE DRC3JM
%***PURPOSE Evaluate the 3j symbol g(M2) = (L1 L2 L3 )
% (M1 M2 -M1-M2)
% for all allowed values of M2, the other parameters
% being held fixed.
%***LIBRARY SLATEC
%***CATEGORY C19
%***TYPE doubleprecision (RC3JM-S, DRC3JM-D)
%***KEYWORDS 3J COEFFICIENTS, 3J SYMBOLS, CLEBSCH-GORDAN COEFFICIENTS,
% RACAH COEFFICIENTS, VECTOR ADDITION COEFFICIENTS,
% WIGNER COEFFICIENTS
%***AUTHOR Gordon, R. G., Harvard University
% Schulten, K., Max Planck Institute
%***DESCRIPTION
%
% *Usage:
%
% doubleprecision L1, L2, L3, M1, M2MIN, M2MAX, THRCOF(NDIM)
% INTEGER NDIM, IER
%
% CALL DRC3JM (L1, L2, L3, M1, M2MIN, M2MAX, THRCOF, NDIM, IER)
%
% *Arguments:
%
% L1 :IN Parameter in 3j symbol.
%
% L2 :IN Parameter in 3j symbol.
%
% L3 :IN Parameter in 3j symbol.
%
% M1 :IN Parameter in 3j symbol.
%
% M2MIN :OUT Smallest allowable M2 in 3j symbol.
%
% M2MAX :OUT Largest allowable M2 in 3j symbol.
%
% THRCOF :OUT Set of 3j coefficients generated by evaluating the
% 3j symbol for all allowed values of M2. THRCOF(I)
% will contain g(M2MIN+I-1), I=1,2,...,M2MAX-M2MIN+1.
%
% NDIM :IN Declared length of THRCOF in calling program.
%
% IER :OUT Error flag.
% IER=0 No errors.
% IER=1 Either L1.LT.ABS(M1) or L1+ABS(M1) non-integer.
% IER=2 ABS(L1-L2).LE.L3.LE.L1+L2 not satisfied.
% IER=3 L1+L2+L3 not an integer.
% IER=4 M2MAX-M2MIN not an integer.
% IER=5 M2MAX less than M2MIN.
% IER=6 NDIM less than M2MAX-M2MIN+1.
%
% *Description:
%
% Although conventionally the parameters of the vector addition
% coefficients satisfy certain restrictions, such as being integers
% or integers plus 1/2, the restrictions imposed on input to this
% subroutine are somewhat weaker. See, for example, Section 27.9 of
% Abramowitz and Stegun or Appendix C of Volume II of A. Messiah.
% The restrictions imposed by this subroutine are
% 1. L1.GE.ABS(M1) and L1+ABS(M1) must be an integer;
% 2. ABS(L1-L2).LE.L3.LE.L1+L2;
% 3. L1+L2+L3 must be an integer;
% 4. M2MAX-M2MIN must be an integer, where
% M2MAX=MIN(L2,L3-M1) and M2MIN=MAX(-L2,-L3-M1).
% If the conventional restrictions are satisfied, then these
% restrictions are met.
%
% The user should be cautious in using input parameters that do
% not satisfy the conventional restrictions. For example, the
% the subroutine produces values of
% g(M2) = (0.75 1.50 1.75 )
% (0.25 M2 -0.25-M2)
% for M2=-1.5,-0.5,0.5,1.5 but none of the symmetry properties of the
% 3j symbol, set forth on page 1056 of Messiah, is satisfied.
%
% The subroutine generates g(M2MIN), g(M2MIN+1), ..., g(M2MAX)
% where M2MIN and M2MAX are defined above. The sequence g(M2) 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 3j 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. Abramowitz, M., and Stegun, I. A., Eds., Handbook
% of Mathematical Functions with Formulas, Graphs
% and Mathematical Tables, NBS Applied Mathematics
% Series 55, June 1964 and subsequent printings.
% 2. Messiah, Albert., Quantum Mechanics, Volume II,
% North-Holland Publishing Company, 1963.
% 3. 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.
% 4. 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.
% 5. Schulten, Klaus and Gordon, Roy G., Recursive
% evaluation of 3j and 6j coefficients, Computer
% Phys Comm, v 11, 1976, pp. 269-278.
%***ROUTINES CALLED D1MACH, XERMSG
%***REVISION HISTORY (YYMMDD)
% 750101 DATE WRITTEN
% 880515 SLATEC prologue added by G. C. Nielson, NBS; parameters
% HUGE and TINY revised to depend on D1MACH.
% 891229 Prologue description rewritten; other prologue sections
% revised; MMATCH (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 THRCOF expanded. These changes were done by
% D. W. Lozier.
%***end PROLOGUE DRC3JM
%
persistent a1 a1s c1 c1old c2 cnorm dv eps firstCall gt hugemlv i indexmlv lstep m2 m3 n newfac nfin nfinp1 nfinp2 nfinp3 nlim nstep2 oldfac one ratio sign1 sign2 srhuge srtiny sum1 sum2 sumbac sumfor sumuni 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(c1), c1=0; end;
if isempty(c1old), c1old=0; end;
if isempty(c2), c2=0; end;
if isempty(cnorm), cnorm=0; end;
if isempty(dv), dv=0; end;
if isempty(eps), eps=0; end;
if isempty(hugemlv), hugemlv=0; end;
if isempty(m2), m2=0; end;
if isempty(m3), m3=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(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.0d0]; end;
if firstCall, eps =[0.01d0]; end;
if firstCall, one =[1.0d0]; end;
if firstCall, two=[2.0d0]; end;
firstCall=0;
%
%***FIRST EXECUTABLE STATEMENT DRC3JM
gt(:)=0;
ier = 0;
% HUGE is the square root of one twentieth of the largest floating
% point number, approximately.
hugemlv = sqrt(d1mach(2)./20.0d0);
srhuge = sqrt(hugemlv);
tinymlv = 1.0d0./hugemlv;
srtiny = 1.0d0./srhuge;
%
% MMATCH = ZERO
%
%
% Check error conditions 1, 2, and 3.
if((l1-abs(m1)+eps<zero) ||(rem(l1+abs(m1)+eps,one)>=eps+eps))
ier = 1;
[dumvar1,dumvar2,dumvar3,ier]=xermsg('SLATEC','DRC3JM',['L1-ABS(M1) less than zero or ','L1+ABS(M1) not integer.'],ier,1);
elseif((l1+l2-l3<-eps) ||(l1-l2+l3<-eps) ||(-l1+l2+l3<-eps) ) ;
ier = 2;
[dumvar1,dumvar2,dumvar3,ier]=xermsg('SLATEC','DRC3JM',['L1, L2, L3 do not satisfy ','triangular condition.'],ier,1);
elseif( rem(l1+l2+l3+eps,one)>=eps+eps ) ;
ier = 3;
[dumvar1,dumvar2,dumvar3,ier]=xermsg('SLATEC','DRC3JM','L1+L2+L3 not integer.',ier,1);
else;
%
%
% Limits for M2
m2min = max(-l2,-l3-m1);
m2max = min(l2,l3-m1);
%
% Check error condition 4.
if( rem(m2max-m2min+eps,one)>=eps+eps )
ier = 4;
[dumvar1,dumvar2,dumvar3,ier]=xermsg('SLATEC','DRC3JM','M2MAX-M2MIN not integer.',ier,1);
elseif( m2min<m2max-eps ) ;
%
% This is reached in case that M1 and M2 take more than one value.
% MSCALE = 0
nfin = fix(fix(m2max-m2min+one+eps));
if( ndim<nfin )
%
% Check error condition 6.
ier = 6;
[dumvar1,dumvar2,dumvar3,ier]=xermsg('SLATEC','DRC3JM',['Dimension of result array for ','3j coefficients too small.'],ier,1);
else;
%
%
%
% Start of forward recursion from M2 = M2MIN
%
m2 = m2min;
thrcof(1) = srtiny;
newfac = 0.0d0;
c1 = 0.0d0;
sum1 = tinymlv;
%
%
lstep = 1;
gt(:)=0;
while( true );
lstep = fix(lstep + 1);
m2 = m2 + one;
m3 = -m1 - m2;
%
%
oldfac = newfac;
a1 =(l2-m2+one).*(l2+m2).*(l3+m3+one).*(l3-m3);
newfac = sqrt(a1);
%
%
dv =(l1+l2+l3+one).*(l2+l3-l1) -(l2-m2+one).*(l3+m3+one)-(l2+m2-one).*(l3-m3-one);
%
%
if( lstep>2 )
c1old = abs(c1);
end;
c1 = -dv./newfac;
%
if( lstep>2 )
%
%
c2 = -oldfac./newfac;
%
% Recursion to the next 3j coefficient
x = c1.*thrcof(lstep-1) + c2.*thrcof(lstep-2);
thrcof(lstep) = x;
sumfor = sum1;
sum1 = sum1 + x.*x;
if( lstep==nfin )
gt(1)=1;
break;
end;
%
% See if last unnormalized 3j coefficient exceeds SRHUGE
%
if( abs(x)>=srhuge )
%
% This is reached if last 3j coefficient larger than SRHUGE,
% so that the recursion series THRCOF(1), ... , THRCOF(LSTEP)
% has to be rescaled to prevent overflow
%
% MSCALE = MSCALE + 1
for i = 1 : lstep;
if( abs(thrcof(i))<srtiny )
thrcof(i) = zero;
end;
thrcof(i) = thrcof(i)./srhuge;
end; i = fix(lstep+1);
sum1 = sum1./hugemlv;
sumfor = sumfor./hugemlv;
x = x./srhuge;
end;
%
%
% As long as ABS(C1) is decreasing, the recursion proceeds towards
% increasing 3j 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 M2 = M2MIN + 1, the third term in the recursion equation vanishes,
% hence
%
x = srtiny.*c1;
thrcof(2) = x;
sum1 = sum1 + tinymlv.*c1.*c1;
if( lstep==nfin )
break;
end;
end;
end;
%
if(gt(1)==0)
sumuni = sum1;
gt(2)=1;
end;
gt(1)=0;
%
%
% Keep three 3j coefficients around MMATCH for comparison later
% with backward recursion values.
%
% MMATCH = M2 - 1
if(gt(2)==0)
nstep2 = fix(nfin - lstep + 3);
x1 = x;
x2 = thrcof(lstep-1);
x3 = thrcof(lstep-2);
%
% Starting backward recursion from M2MAX taking NSTEP2 steps, so
% that forwards and backwards recursion overlap at the three points
% M2 = MMATCH+1, MMATCH, MMATCH-1.
%
nfinp1 = fix(nfin + 1);
nfinp2 = fix(nfin + 2);
nfinp3 = fix(nfin + 3);
thrcof(nfin) = srtiny;
sum2 = tinymlv;
%
%
%
m2 = m2max + two;
lstep = 1;
while( true );
lstep = fix(lstep + 1);
m2 = m2 - one;
m3 = -m1 - m2;
oldfac = newfac;
a1s =(l2-m2+two).*(l2+m2-one).*(l3+m3+two).*(l3-m3-one);
newfac = sqrt(a1s);
dv =(l1+l2+l3+one).*(l2+l3-l1) -(l2-m2+one).*(l3+m3+one)-(l2+m2-one).*(l3-m3-one);
c1 = -dv./newfac;
if( lstep>2 )
%
c2 = -oldfac./newfac;
%
% Recursion to the next 3j coefficient
%
y = c1.*thrcof(nfinp2-lstep) + c2.*thrcof(nfinp3-lstep);
%
if( lstep==nstep2 )
break;
end;
%
thrcof(nfinp1-lstep) = y;
sumbac = sum2;
sum2 = sum2 + y.*y;
%
%
% See if last 3j coefficient exceeds SRHUGE
%
if( abs(y)>=srhuge )
%
% This is reached if last 3j coefficient larger than SRHUGE,
% so that the recursion series THRCOF(NFIN), ... , THRCOF(NFIN-LSTEP+1)
% has to be rescaled to prevent overflow.
%
% MSCALE = MSCALE + 1
for i = 1 : lstep;
indexmlv = nfin - i + 1;
if( abs(thrcof(indexmlv))<srtiny )
thrcof(indexmlv) = zero;
end;
thrcof(indexmlv) = thrcof(indexmlv)./srhuge;
end; i = fix(lstep+1);
sum2 = sum2./hugemlv;
%
sumbac = sumbac./hugemlv;
end;
else;
%
% If M2 = M2MAX + 1 the third term in the recursion equation vanishes
%
y = srtiny.*c1;
thrcof(nfin-1) = y;
if( lstep==nstep2 )
break;
end;
sumbac = sum2;
sum2 = sum2 + y.*y;
end;
end;
%
%
%
% The forward recursion 3j coefficients X1, X2, X3 are to be matched
% with the corresponding backward recursion values Y1, Y2, Y3.
%
y3 = y;
y2 = thrcof(nfinp2-lstep);
y1 = thrcof(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;
thrcof(n) = ratio.*thrcof(n);
end; n = fix(nfin+1);
sumuni = sumfor + ratio.*ratio.*sumbac;
else;
%
for n = 1 : nlim;
thrcof(n) = ratio.*thrcof(n);
end; n = fix(nlim+1);
sumuni = ratio.*ratio.*sumfor + sumbac;
end;
end;
gt(2)=0;
%
%
% Normalize 3j coefficients
%
cnorm = one./sqrt((l1+l1+one).*sumuni);
%
% Sign convention for last 3j coefficient determines overall phase
%
sign1 = (abs(one).*sign(thrcof(nfin)));
sign2 =(-one).^fix(abs(l2-l3-m1)+eps);
if( sign1.*sign2<=0 )
cnorm = -cnorm;
end;
%
if( abs(cnorm)<one )
%
thresh = tinymlv./abs(cnorm);
for n = 1 : nfin;
if( abs(thrcof(n))<thresh )
thrcof(n) = zero;
end;
thrcof(n) = cnorm.*thrcof(n);
end; n = fix(nfin+1);
else;
%
for n = 1 : nfin;
thrcof(n) = cnorm.*thrcof(n);
end; n = fix(nfin+1);
end;
end;
elseif( m2min<m2max+eps ) ;
%
%
% This is reached in case that M2 and M3 can take only one value.
% MSCALE = 0
thrcof(1) =(-one).^fix(abs(l2-l3-m1)+eps)./sqrt(l1+l2+l3+one);
else;
%
% Check error condition 5.
ier = 5;
[dumvar1,dumvar2,dumvar3,ier]=xermsg('SLATEC','DRC3JM','M2MIN greater than M2MAX.',ier,1);
end;
end;
%
%
%
end %subroutine drc3jm
%DECK DRC6J
|
|
