function f=pnmars(l,m)
% return the Marshak boundary coefficients in slab geometry. These
% coefficients are specific to the left boundary.
% function f=pnmars(l,m)
% (c) 2008 Alain Hebert, Ecole Polytechnique de Montreal
  if mod(l,2) == 0 
    error('odd first index expected')
  end
  zgksi=[ 0.00529953837  0.0277124941  0.0671843886  0.122297794  0.191061884 ...
  0.270991623  0.359198213  0.452493757  0.547506273  0.640801787  0.729008377 ...
  0.808938146  0.877702236  0.932815611  0.972287536  0.994700432 ] ;
  wgksi=[ 0.01357623  0.0311267618  0.047579255  0.0623144843  0.0747979954 ...
  0.0845782608  0.0913017094  0.0947253034  0.0947253034  0.0913017094 ...
  0.0845782608  0.0747979954  0.0623144843  0.047579255  0.0311267618 ...
  0.01357623 ] ;
  f=(2*m+1)*sum(wgksi(:).*plgndr(l,0,zgksi(:)).*plgndr(m,0,zgksi(:))) ;