function [lun,kprint,ipass]=qc36j(lun,kprint,ipass);
persistent diff first firstCall fmt fmt2 i ier ierjj ierjm iil1 iim2 indexmlv ipass1 ipass2 ipass3 ipass4 ipass5 jjval jmval l1 l1max l1min l2 l3 last m1 m2 m2max m2min m3 ndim nerr nsig r3jj r3jm r6j sixcof string thrcof tol x ; if isempty(firstCall),firstCall=1;end;
if isempty(iil1), iil1=0; end;
if isempty(iim2), iim2=0; end;
if isempty(string), string=repmat(' ',1,36); end;
if isempty(fmt), fmt=repmat(' ',1,30); end;
if isempty(fmt2), fmt2=repmat(' ',1,13); end;
if isempty(ipass1), ipass1=0; end;
if isempty(ipass2), ipass2=0; end;
if isempty(ipass3), ipass3=0; end;
if isempty(ipass4), ipass4=0; end;
if isempty(ipass5), ipass5=0; end;
if isempty(ier), ier=0; end;
if isempty(indexmlv), indexmlv=0; end;
if isempty(i), i=0; end;
if isempty(first), first=0; end;
if isempty(last), last=0; end;
if isempty(nsig), nsig=0; end;
if isempty(nerr), nerr=0; end;
if isempty(ierjj), ierjj=0; end;
if isempty(ierjm), ierjm=0; end;
if isempty(ndim), ndim=15 ; end;
if isempty(tol), tol=0; end;
if isempty(l1), l1=0; end;
if isempty(l2), l2=0; end;
if isempty(l3), l3=0; end;
if isempty(m1), m1=0; end;
if isempty(m2), m2=0; end;
if isempty(m3), m3=0; end;
if isempty(l1min), l1min=0; end;
if isempty(l1max), l1max=0; end;
if isempty(m2min), m2min=0; end;
if isempty(m2max), m2max=0; end;
if isempty(diff), diff=zeros(1,ndim); end;
if isempty(x), x=0; end;
if isempty(jjval), jjval=0; end;
if isempty(jmval), jmval=0; end;
if isempty(thrcof), thrcof=zeros(1,ndim); end;
if isempty(sixcof), sixcof=zeros(1,ndim); end;
if isempty(r3jj), r3jj=zeros(1,8); end;
if isempty(r3jm), r3jm=zeros(1,14); end;
if isempty(r6j), r6j=zeros(1,15); end;
if firstCall, r3jj=[2.78886675511358515993e-1,-9.53462589245592315447e-2,-6.74199862463242086246e-2,1.53311035167966641297e-1,-1.56446554693685969725e-1,1.09945041215655051079e-1,-5.53623569313171943334e-2,1.79983545113778583298e-2]; end;
if firstCall, r3jm=[2.09158973288615242614e-2,8.53756555321524722127e-2,9.08295370868692516943e-2,-3.89054377846499391700e-2,-6.63734970165680635691e-2,6.49524040528389395031e-2,2.15894310595403759392e-2,-7.78912711785239219992e-2,3.59764371059543401880e-2,5.47301500021263423079e-2,-7.59678665956761514629e-2,-2.19224445539892113776e-2,1.01167744280772202424e-1,7.34825726244719704696e-2]; end;
if firstCall, r6j=[3.49090513837329977746e-2,-3.74302503965979160859e-2,1.89086639095956018415e-2,7.34244825492864345709e-3,-2.35893518508179445858e-2,1.91347695521543652000e-2,1.28801739772417220844e-3,-1.93001836629052653977e-2,1.67730594938288876974e-2,5.50114727485094871674e-3,-2.13543979089683097421e-2,3.46036445143538730828e-3,2.52095005479558458604e-2,1.48399056122171330285e-2,2.70857768063318559724e-3]; end;
firstCall=0;
tol = 100.0.*r1mach(3);
if( kprint>=2 )
writef(lun,['%s \n'], [' THIS IS QC36J, A TEST PROGRAM FOR THE ','SINGLE PRECISION 3J6J PACKAGE.']);
writef(lun,['%s \n'], [' AN EXPLANATION OF THE VARIOUS ','TESTS CAN BE FOUND IN THE PROGRAM COMMENTS.']);
writef(lun,['%0.15g \n']);
end;
x = 1.0./3.0;
string=sprintf(['%35.25f' ' \n'], x);
%format(f35.25);
for i = 1 : 35;
if( strcmp(deblank(string([i:i])),deblank('3')) )
first = fix(i);
break;
end;
end;
for i = first : 35;
if( ~strcmp(deblank(string([i:i])),deblank('3')) )
last = fix(i - 1);
break;
end;
last = 36;
end;
nsig = fix(last - first + 1);
fmt([1:16]) = '(1X,F5.1,T8,G35.';
fmt([17:min(length(fmt),18)])=sprintf(['%2i'], nsig);
fmt([19:27]) = ',T45,G35.';
fmt([28:min(length(fmt),29)])=sprintf(['%2i'], nsig);
fmt([30:30]) = ')';
fmt2([1:10]) = '(1X,A,G35.';
fmt2([11:min(length(fmt2),12)])=sprintf(['%2i'], nsig);
fmt2([13:13]) = ')';
ipass1 = 1;
l2 = 4.5;
l3 = 3.5;
m2 = -3.5;
m3 = 2.5;
[l2,l3,m2,m3,l1min,l1max,thrcof,ndim,ier]=rc3jj(l2,l3,m2,m3,l1min,l1max,thrcof,ndim,ier);
if( ier~=0 )
ipass1 = 0;
else;
for iil1 = fix(l1min) : fix(l1max);
indexmlv = fix(iil1-l1min) + 1;
m1 = 1.0;
diff(indexmlv) = abs(thrcof(indexmlv)-r3jj(indexmlv));
if( diff(indexmlv)>abs(r3jj(indexmlv)).*tol )
ipass1 = 0;
end;
end; iil1 = fix(fix(l1max)+1);
l1 = iil1;
end;
if( kprint>=3 ||(kprint==2 && ipass1==0) )
writef(lun,['%s %s \n'], ' TEST 1, RECURRENCE IN L1, COMPARE VALUES OF 3J ', 'CALCULATED BY RC3JJ TO');
writef(lun,['%s %s \n'], ' VALUES CALCULATED BY EXPLICIT FORMULA FROM ' ,'MESSIAH''S QUANTUM MECHANICS');
writef(lun,[' L2 = ','%5.1f',' L3 = ','%5.1f' ' \n'], l2 , l3);
writef(lun,[' M1 = ','%5.1f',' M2 = ','%5.1f',' M3 = ','%5.1f' ' \n'], m1 , m2 , m3);
if( ier~=0 )
writef(lun,['%s %s %0.15g \n'], ' ERROR RETURNED FROM SUBROUTINE ' ,'RC3JJ: IER =' , ier);
else;
writef(lun,[' L1','%31f',' RC3JJ VALUE','%67f','FORMULA VALUE' ' \n']);
%format (' L1',t31,' RC3JJ VALUE',t67,'FORMULA VALUE');
for iil1 = fix(l1min) : fix(l1max);
indexmlv = fix(iil1-l1min) + 1;
disp({ iil1 , thrcof(indexmlv) , r3jj(indexmlv)});
if( diff(indexmlv)>abs(r3jj(indexmlv)).*tol )
writef(lun,[repmat(' ',1,1),'%s','%5.1f'], ['DIFFERENCE EXCEEDS ERROR ','TOLERANCE FOR iiL1 ='] , iil1);
end;
end; iil1 = fix(fix(l1max)+1);
l1 = iil1;
end;
end;
if( ipass1==0 )
if( kprint>=1 )
writef(lun,['%s \n'], ' TEST 1 FAILED ');
writef(lun,['%0.15g \n']);
end;
elseif( kprint>=2 ) ;
writef(lun,['%s \n'], ' TEST 1 PASSED ');
writef(lun,['%0.15g \n']);
end;
ipass2 = 1;
l1 = 8.0;
l2 = 7.5;
l3 = 6.5;
m1 = 1.0;
[l1,l2,l3,m1,m2min,m2max,thrcof,ndim,ier]=rc3jm(l1,l2,l3,m1,m2min,m2max,thrcof,ndim,ier);
if( ier~=0 )
ipass2 = 0;
else;
for iim2 = fix(m2min) : fix(m2max);
indexmlv = fix(iim2-m2min) + 1;
m3 = -m1 - iim2;
diff(indexmlv) = abs(thrcof(indexmlv)-r3jm(indexmlv));
if( diff(indexmlv)>abs(r3jm(indexmlv)).*tol )
ipass2 = 0;
end;
end; iim2 = fix(fix(m2max)+1);
m2 = iim2;
end;
if( kprint>=3 ||(kprint==2 && ipass2==0) )
writef(lun,['%s %s \n'], ' TEST 2, RECURRENCE IN M2, COMPARE VALUES OF 3J ', 'CALCULATED BY RC3JM TO');
writef(lun,['%s %s \n'], ' VALUES CALCULATED BY EXPLICIT FORMULA FROM ' ,'MESSIAH''S QUANTUM MECHANICS');
writef(lun,[' L1 = ','%5.1f',' L2 = ','%5.1f',' L3 = ','%5.1f' ' \n'], l1 , l2 , l3);
%format (' L1 = ',f5.1,' L2 = ',f5.1,' L3 = ',f5.1);
writef(lun,[' M1 = ','%5.1f',' M3 = -(M1+M2)' ' \n'], m1);
%format (' M1 = ',f5.1,' M3 = -(M1+M2)');
if( ier~=0 )
writef(lun,['%s %s %0.15g \n'], ' ERROR RETURNED FROM SUBROUTINE ' ,'RC3JM: IER =' , ier);
else;
writef(lun,[' M2','%31f',' RC3JM VALUE','%67f','FORMULA VALUE' ' \n']);
%format (' M2',t31,' RC3JM VALUE',t67,'FORMULA VALUE');
for iim2 = fix(m2min) : fix(m2max);
indexmlv = fix(iim2-m2min) + 1;
disp({ iim2 , thrcof(indexmlv) , r3jm(indexmlv)});
if( diff(indexmlv)>abs(r3jm(indexmlv)).*tol )
writef(lun,[repmat(' ',1,1),'%s','%5.1f'], ['DIFFERENCE EXCEEDS ERROR ','TOLERANCE FOR iiM2 ='] , iim2);
end;
end; iim2 = fix(fix(m2max)+1);
m2 = iim2;
end;
end;
if( ipass2==0 )
if( kprint>=1 )
writef(lun,['%s \n'], ' TEST 2 FAILED ');
writef(lun,['%0.15g \n']);
end;
elseif( kprint>=2 ) ;
writef(lun,['%s \n'], ' TEST 2 PASSED ');
writef(lun,['%0.15g \n']);
end;
ipass3 = 1;
l1 = 100.0;
l2 = 2.0;
l3 = 100.0;
m1 = -10.0;
m2 = 0.0;
m3 = 10.0;
[l2,l3,m2,m3,l1min,l1max,thrcof,ndim,ierjj]=rc3jj(l2,l3,m2,m3,l1min,l1max,thrcof,ndim,ierjj);
jjval = thrcof(3);
[l1,l2,l3,m1,m2min,m2max,thrcof,ndim,ierjm]=rc3jm(l1,l2,l3,m1,m2min,m2max,thrcof,ndim,ierjm);
jmval = thrcof(3);
if( ierjj~=0 || ierjm~=0 )
ipass3 = 0;
else;
diff(1) = abs(jjval-jmval);
if( diff(1)>0.5.*abs(jjval+jmval).*tol )
ipass3 = 0;
end;
end;
if( kprint>=3 ||(kprint==2 && ipass3==0) )
writef(lun,['%s %s \n'], ' TEST 3, COMPARE A COMMON VALUE CALCULATED BY ' ,'BOTH RC3JJ AND RC3JM');
writef(lun,['%s \n'], ' L1 = 100.0 L2 = 2.0 L3 = 100.0');
writef(lun,['%s \n'], ' M1 = -10.0 M2 = 0.0 M3 = 10.0');
if( ierjj~=0 )
writef(lun,['%s %s %0.15g \n'], ' ERROR RETURNED FROM SUBROUTINE ' ,'RC3JJ: IER =' , ierjj);
elseif( ierjm~=0 ) ;
writef(lun,['%s %s %0.15g \n'], ' ERROR RETURNED FROM SUBROUTINE ' ,'RC3JM: IER =' , ierjm);
else;
disp({ 'RC3JJ VALUE =' , jjval});
disp({ 'RC3JM VALUE =' , jmval});
if( diff(1)>0.5.*abs(jjval+jmval).*tol )
writef(lun,[repmat(' ',1,1),'%s'],'DIFFERENCE EXCEEDS ERROR TOLERANCE');
end;
end;
end;
if( ipass3==0 )
if( kprint>=1 )
writef(lun,['%s \n'], ' TEST 3 FAILED ');
writef(lun,['%0.15g \n']);
end;
elseif( kprint>=2 ) ;
writef(lun,['%s \n'], ' TEST 3 PASSED ');
writef(lun,['%0.15g \n']);
end;
ipass4 = 1;
l2 = 8.0;
l3 = 7.0;
m1 = 6.5;
m2 = 7.5;
m3 = 7.5;
[l2,l3,m1,m2,m3,l1min,l1max,sixcof,ndim,ier]=rc6j(l2,l3,m1,m2,m3,l1min,l1max,sixcof,ndim,ier);
if( ier~=0 )
ipass4 = 0;
else;
for iil1 = fix(l1min) : fix(l1max);
indexmlv = fix(iil1-l1min) + 1;
diff(indexmlv) = abs(sixcof(indexmlv)-r6j(indexmlv));
if( diff(indexmlv)>abs(r6j(indexmlv)).*tol )
ipass4 = 0;
end;
end; iil1 = fix(fix(l1max)+1);
l1 = iil1;
end;
if( kprint>=3 ||(kprint==2 && ipass4==0) )
writef(lun,['%s %s \n'], ' TEST 4, RECURRENCE IN L1, COMPARE VALUES OF 6J ', 'CALCULATED BY RC6J TO');
writef(lun,['%s %s \n'], ' VALUES CALCULATED BY EXPLICIT FORMULA FROM ' ,'MESSIAH''S QUANTUM MECHANICS');
writef(lun,[' L2 = ','%5.1f',' L3 = ','%5.1f' ' \n'], l2 , l3);
writef(lun,[' M1 = ','%5.1f',' M2 = ','%5.1f',' M3 = ','%5.1f' ' \n'], m1 , m2 , m3);
if( ier~=0 )
writef(lun,['%s %s %0.15g \n'], ' ERROR RETURNED FROM SUBROUTINE ' ,'RC6J: IER =' , ier);
else;
writef(lun,[' L1','%32f',' RC6J VALUE','%67f','FORMULA VALUE' ' \n']);
%format (' L1',t32,' RC6J VALUE',t67,'FORMULA VALUE');
for iil1 = fix(l1min) : fix(l1max);
indexmlv = fix(iil1-l1min) + 1;
disp({ iil1 , sixcof(indexmlv) , r6j(indexmlv)});
if( diff(indexmlv)>abs(r6j(indexmlv)).*tol )
writef(lun,[repmat(' ',1,1),'%s','%5.1f'], ['DIFFERENCE EXCEEDS ERROR ','TOLERANCE FOR iiL1 ='] , iil1);
end;
end; iil1 = fix(fix(l1max)+1);
l1 = iil1;
end;
end;
if( ipass4==0 )
if( kprint>=1 )
writef(lun,['%s \n'], ' TEST 4 FAILED ');
writef(lun,['%0.15g \n']);
end;
elseif( kprint>=2 ) ;
writef(lun,['%s \n'], ' TEST 4 PASSED ');
writef(lun,['%0.15g \n']);
end;
ipass5 = 1;
if( kprint<=2 )
xsetf(0);
else;
xsetf(-1);
end;
if( kprint>=3 )
writef(lun,['%s %s \n'],' TEST 5, CHECK FOR PROPER HANDLING ', 'OF INVALID INPUT');
end;
l2 = 2.0;
l3 = 100.0;
m1 = -6.0;
m2 = -4.0;
m3 = 10.0;
if( kprint>=3 )
writef(lun,['%0.15g \n']);
end;
xerclr;
[l2,l3,m2,m3,l1min,l1max,thrcof,ndim,ier]=rc3jj(l2,l3,m2,m3,l1min,l1max,thrcof,ndim,ier);
if( numxer(nerr)~=ier )
ipass5 = 0;
end;
l2 = 2.0;
l3 = 99.5;
m1 = -10.0;
m2 = 0.0;
m3 = 10.0;
if( kprint>=3 )
writef(lun,['%0.15g \n']);
end;
xerclr;
[l2,l3,m2,m3,l1min,l1max,thrcof,ndim,ier]=rc3jj(l2,l3,m2,m3,l1min,l1max,thrcof,ndim,ier);
if( numxer(nerr)~=ier )
ipass5 = 0;
end;
l2 = 3.2;
l3 = 4.5;
m1 = -1.3;
m2 = 0.8;
m3 = 0.5;
if( kprint>=3 )
writef(lun,['%0.15g \n']);
end;
xerclr;
[l2,l3,m2,m3,l1min,l1max,thrcof,ndim,ier]=rc3jj(l2,l3,m2,m3,l1min,l1max,thrcof,ndim,ier);
if( numxer(nerr)~=ier )
ipass5 = 0;
end;
l2 = 10.0;
l3 = 150.0;
m1 = -10.0;
m2 = 0.0;
m3 = 10.0;
if( kprint>=3 )
writef(lun,['%0.15g \n']);
end;
xerclr;
[l2,l3,m2,m3,l1min,l1max,thrcof,ndim,ier]=rc3jj(l2,l3,m2,m3,l1min,l1max,thrcof,ndim,ier);
if( numxer(nerr)~=ier )
ipass5 = 0;
end;
l1 = 100.0;
l2 = 2.0;
l3 = 100.0;
m1 = 150.0;
if( kprint>=3 )
writef(lun,['%0.15g \n']);
end;
xerclr;
[l1,l2,l3,m1,m2min,m2max,thrcof,ndim,ier]=rc3jm(l1,l2,l3,m1,m2min,m2max,thrcof,ndim,ier);
if( numxer(nerr)~=ier )
ipass5 = 0;
end;
l1 = 20.0;
l2 = 5.0;
l3 = 10.0;
m1 = -10.0;
if( kprint>=3 )
writef(lun,['%0.15g \n']);
end;
xerclr;
[l1,l2,l3,m1,m2min,m2max,thrcof,ndim,ier]=rc3jm(l1,l2,l3,m1,m2min,m2max,thrcof,ndim,ier);
if( numxer(nerr)~=ier )
ipass5 = 0;
end;
l1 = 1.0;
l2 = 1.3;
l3 = 1.5;
m1 = 0.0;
if( kprint>=3 )
writef(lun,['%0.15g \n']);
end;
xerclr;
[l1,l2,l3,m1,m2min,m2max,thrcof,ndim,ier]=rc3jm(l1,l2,l3,m1,m2min,m2max,thrcof,ndim,ier);
if( numxer(nerr)~=ier )
ipass5 = 0;
end;
l1 = 1.0;
l2 = 1.3;
l3 = 1.7;
m1 = 0.0;
if( kprint>=3 )
writef(lun,['%0.15g \n']);
end;
xerclr;
[l1,l2,l3,m1,m2min,m2max,thrcof,ndim,ier]=rc3jm(l1,l2,l3,m1,m2min,m2max,thrcof,ndim,ier);
if( numxer(nerr)~=ier )
ipass5 = 0;
end;
l1 = 100.0;
l2 = 10.0;
l3 = 110.0;
m1 = -10.0;
if( kprint>=3 )
writef(lun,['%0.15g \n']);
end;
xerclr;
[l1,l2,l3,m1,m2min,m2max,thrcof,ndim,ier]=rc3jm(l1,l2,l3,m1,m2min,m2max,thrcof,ndim,ier);
if( numxer(nerr)~=ier )
ipass5 = 0;
end;
l2 = 0.5;
l3 = 1.0;
m1 = 0.5;
m2 = 2.0;
m3 = 3.0;
if( kprint>=3 )
writef(lun,['%0.15g \n']);
end;
xerclr;
[l2,l3,m1,m2,m3,l1min,l1max,sixcof,ndim,ier]=rc6j(l2,l3,m1,m2,m3,l1min,l1max,sixcof,ndim,ier);
if( numxer(nerr)~=ier )
ipass5 = 0;
end;
l2 = 1.0;
l3 = 3.0;
m1 = 5.0;
m2 = 6.0;
m3 = 2.0;
if( kprint>=3 )
writef(lun,['%0.15g \n']);
end;
xerclr;
[l2,l3,m1,m2,m3,l1min,l1max,sixcof,ndim,ier]=rc6j(l2,l3,m1,m2,m3,l1min,l1max,sixcof,ndim,ier);
if( numxer(nerr)~=ier )
ipass5 = 0;
end;
l2 = 4.0;
l3 = 1.0;
m1 = 5.0;
m2 = 3.0;
m3 = 2.0;
if( kprint>=3 )
writef(lun,['%0.15g \n']);
end;
xerclr;
[l2,l3,m1,m2,m3,l1min,l1max,sixcof,ndim,ier]=rc6j(l2,l3,m1,m2,m3,l1min,l1max,sixcof,ndim,ier);
if( numxer(nerr)~=ier )
ipass5 = 0;
end;
l2 = 0.9;
l3 = 0.5;
m1 = 0.9;
m2 = 0.4;
m3 = 0.2;
if( kprint>=3 )
writef(lun,['%0.15g \n']);
end;
xerclr;
[l2,l3,m1,m2,m3,l1min,l1max,sixcof,ndim,ier]=rc6j(l2,l3,m1,m2,m3,l1min,l1max,sixcof,ndim,ier);
if( numxer(nerr)~=ier )
ipass5 = 0;
end;
l2 = 50.0;
l3 = 25.0;
m1 = 15.0;
m2 = 30.0;
m3 = 40.0;
if( kprint>=3 )
writef(lun,['%0.15g \n']);
end;
xerclr;
[l2,l3,m1,m2,m3,l1min,l1max,sixcof,ndim,ier]=rc6j(l2,l3,m1,m2,m3,l1min,l1max,sixcof,ndim,ier);
if( numxer(nerr)~=ier )
ipass5 = 0;
end;
if( ipass5==0 )
if( kprint>=1 )
writef(lun,['%s \n'], ' TEST 5 FAILED ');
writef(lun,['%0.15g \n']);
end;
elseif( kprint>=2 ) ;
writef(lun,['%s \n'], ' TEST 5 PASSED ');
writef(lun,['%0.15g \n']);
end;
if((ipass1==0) ||(ipass2==0) ||(ipass3==0) ||(ipass4==0) ||(ipass5==0) )
ipass = 0;
if( kprint>=1 )
writef(lun,[' QC36J FAILED SOME TESTS ' ' \n']);
end;
%format (' QC36J FAILED SOME TESTS ');
else;
ipass = 1;
if( kprint>=2 )
writef(lun,[' QC36J PASSED ALL TESTS ' ' \n']);
end;
%format (' QC36J PASSED ALL TESTS ');
end;
%format (' L2 = ',f5.1,' L3 = ',f5.1);
%format (' M1 = ',f5.1,' M2 = ',f5.1,' M3 = ',f5.1);
end %subroutine qc36j
