function f=pncmar(l1,l2,m1,m2)
% return the Marshak boundary coefficients in 1D tube geometry.
% function f=pncmar(l1,l2,m1,m2)
% (c) 2008 Alain Hebert, Ecole Polytechnique de Montreal
  if mod(l1,2) == 0
     error('odd legendre order expected.')
  elseif mod(m1,2) == 0
     error('odd azimuthal order expected.')
  end
  zgksi=[  -0.989400923 -0.944575012 -0.865631223 -0.755404413 -0.617876232 ...
  -0.458016783 -0.281603545 -0.0950125083  0.0950125083  0.281603545 ...
  0.458016783  0.617876232  0.755404413  0.865631223  0.944575012  0.989400923 ] ;
  wgksi=[  0.0271524601  0.0622535236  0.0951585099  0.124628969  0.149595991 ...
  0.169156522  0.182603419  0.189450607  0.189450607  0.182603419  0.169156522 ...
  0.149595991  0.124628969  0.0951585099  0.0622535236  0.0271524601 ] ;
  coef=1.0 ;
  if m1 > 0 , coef=coef*sqrt(2.0*factorial(l1-m1)/factorial(l1+m1)) ; end
  if m2 > 0 , coef=coef*sqrt(2.0*factorial(l2-m2)/factorial(l2+m2)) ; end
  sum1=sum(wgksi(:).*plgndr(l1,m1,zgksi(:)).*plgndr(l2,m2,zgksi(:))) ;
  if m1 == m2
    sumphi=0.25 ;
  elseif mod(m1+m2,2) == 0
    sumphi=0.0 ;
  else
    sumphi=(-1.)^fix((m1-m2)/2)*sign(m1-m2)*(1./(m1-m2)+1./(m1+m2))/(2.*pi) ;
  end
  f=sum1*coef*sumphi*(2.*l2+1.) ;