function [tpq,upq,vpq,wpq]=snqu02(nlf)
% set the type-2 level-symmetric quadrature base points and weights in
% one octant. See "Applied Reactor Physics", Table 3.4
% function [tpq,upq,vpq,wpq)]=snqu02(nlf)
% (c) 2010 Alain Hebert, Ecole Polytechnique de Montreal
   %----
   % set the unique quadrature values.
   %----
   m2=nlf/2 ; nw=fix(1+(nlf*(nlf+8)-1)/48) ;
   if nlf == 2
      zmu1=1/3 ;
   elseif nlf == 4
      zmu1=0.12251480 ;
   elseif nlf == 6
      zmu1=0.07109447 ;
   elseif nlf == 8
      zmu1=0.04761903 ;
   elseif nlf == 10
      zmu1=0.03584310 ;
   elseif nlf == 12
      zmu1=0.02796615 ;
   elseif nlf == 14
      zmu1=0.02310250 ;
   elseif nlf == 16
      zmu1=0.01931398 ;
   elseif nlf == 18
      zmu1=0.01692067 ;
   else
      error('order not available')
   end
   u=zeros(1,m2) ; u(1)=sqrt(zmu1) ;
   for i=2:m2
      u(i)=sqrt(zmu1+2*(i-1)*(1-3*zmu1)/(nlf-2)) ;
   end
   %----
   % compute the position of weights.
   %----
   ipr=0 ; inmax=0 ;
   jop=zeros(1,m2) ;
   npq=nlf*(nlf/2+1)/4 ; tpq=zeros(1,npq) ; upq=tpq ; vpq=tpq ; wpq=tpq ;
   inwei=zeros(1,nw) ;
   for ip=1:m2 ;
      jop(ip)=m2-ip+1 ;
      for iq=1:jop(ip)
         ipr=ipr+1 ;
         tpq(ipr)=u(ip) ; upq(ipr)=u(m2+2-ip-iq) ; vpq(ipr)=u(iq) ;
         is=min(min(ip,iq),m2+2-ip-iq) ;
         nw0=0 ;
         for ii=1:is-1
            nw0=nw0+fix((m2-3*(ii-1)+1)/2) ;
         end
         kk=ip-is+1 ; ll=iq-is+1 ;
         if kk == 1
            inwei(ipr)=nw0+min(ll,m2-3*(is-1)+1-ll) ;
         elseif ll == 1
            inwei(ipr)=nw0+min(kk,m2-3*(is-1)+1-kk) ;
         else
            inwei(ipr)=nw0+min(kk,ll) ;
         end
         inmax=max(inmax,inwei(ipr)) ;
      end
   end
   if inmax ~= nw , error('invalid value of nw') ; end ;
   if ipr ~= npq , error('invalid value of npq') ; end ;
   %----
   % set the rectangular system and solve it using the QR method.
   %----
   neq=0;
   for ipl=0:2:nlf
     for ipk=ipl:2:nlf-ipl
       if mod(ipl+ipk,2) == 1 , continue ; end ;
       neq=neq+1 ;
     end
   end
   zmat=zeros(neq,nw+1) ; ieq=0 ;
   for ipl=0:2:nlf
     for ipk=ipl:2:nlf-ipl
       if mod(ipl+ipk,2) == 1 , continue ; end ;
       ieq=ieq+1 ;
       zmat(ieq,1:nw)=0 ;
       for ipq=1:npq
          zmu=tpq(ipq) ; zeta=upq(ipq) ;
          iw=inwei(ipq) ;
          zmat(ieq,iw)=zmat(ieq,iw)+(zmu^ipk)*(zeta^ipl) ;
       end
       ref=1/(ipk+ipl+1) ;
       for i=1:2:ipl-1
          ref=ref*i/(ipk+i) ;
       end
       zmat(ieq,nw+1)=ref ;
     end
   end
   if ieq ~= neq , error('invalid value of neq') ; end ;
   [Q,R]=qr(zmat(:,1:nw),0) ; wei=R\Q'*zmat(:,nw+1) ;
   %----
   % set the level-symmetric quadratures.
   %----
   ipq=0 ;
   for ip=1:m2
      for iq=1:jop(ip)
         ipq=ipq+1 ;
         wpq(ipq)=wei(inwei(ipq))*pi/2 ;
      end
   end