function coef=pn2div(l,m)
% return the real spherical harmonics difference coefficients.
% Reference: Nucl. Sci. Eng., 154, 134 (2006).
% function coef=pn2div(l,m)
% (c) 2010 Alain Hebert, Ecole Polytechnique de Montreal
  if l < 0 ; error('pn2div: invalid first parameter.') ; end
  if (m < -l) || (m > l) ; error('pn2div: invalid second parameter.') ; end
  coef=zeros(1,8) ;
  coef(1)=sqrt((l+m)*(l+m-1)) ;
  coef(2)=-sqrt((l-m+1)*(l-m+2)) ;
  coef(3)=-sqrt((l-m)*(l-m-1)) ;
  coef(4)=sqrt((l+m+1)*(l+m+2)) ;
  coef(5)=-sqrt((l+m)*(l+m-1)) ;
  coef(6)=sqrt((l-m+1)*(l-m+2)) ;
  coef(7)=-sqrt((l-m)*(l-m-1)) ;
  coef(8)=sqrt((l+m+1)*(l+m+2)) ;
  if m == -1
    coef(3)=0 ; coef(4)=0 ;
    coef(7)=coef(7)*sqrt(2) ;
    coef(8)=coef(8)*sqrt(2) ;
  elseif m == 0
    coef(1)=0 ; coef(2)=0 ;
    coef(3)=coef(3)*sqrt(2) ;
    coef(4)=coef(4)*sqrt(2) ;
    coef(5)=0 ; coef(6)=0 ;
    coef(7)=coef(7)*sqrt(2) ;
    coef(8)=coef(8)*sqrt(2) ;
  elseif m == 1
    coef(1)=coef(1)*sqrt(2) ;
    coef(2)=coef(2)*sqrt(2) ;
    coef(5)=0 ; coef(6)=0 ;
  end
  coef(1:8)=coef(1:8)./(2*(2*l+1)) ;
  if m < 0 ; coef(1:8)=-coef(1:8) ; end