function f=pij_slab(x,sig,bc,alb)
% compute the reduced collision probability matrix in 1D slab geometry
% function f=pij_slab(x,sig,bc,alb)
% (c) 2005 Alain Hebert, Polytechnique Montreal
  npij=size(x,2)-1 ;
  vol=zeros(1,npij) ;
  for i=1:npij
    vol(i)=x(i+1)-x(i) ;
  end
  f=zeros(npij) ;
  %
  switch bc
    case {'void'}
      if alb(1) ~= 0.0
        error('left albedo equal to zero expected')
      elseif alb(2) ~= 0.0
        error('right albedo equal to zero expected')
      end
      for ip=1:npij
        f(ip,ip)=rii_f(1.0,0.0,sig(ip),vol(ip)) ;
        tau0=0.0 ;
        for jp=ip+1:npij
          f(ip,jp)=rij_f(1.0,tau0,sig(ip),sig(jp),vol(ip),vol(jp)) ;
          tau0=tau0+vol(jp)*sig(jp) ;
        end
      end
      for i=1:npij
        for j=i:npij
          f(j,i)=f(i,j) ;
        end
      end
    case {'tran'}
      if alb(1) ~= alb(2)
        error('left albedo equal to right albedo expected')
      end
%     complete this case !!!!!!!!!!
    case {'albe'}
      error('not requested for ENE6101')
    otherwise
      error('unknown type of boundary condition')
  end
  f=inv(diag(vol))*f ;
%
  if alb(1)==1 && alb(2)==1
    conservation=max(abs(ones(npij,1)-sum(f*diag(sig),2))) ;
    fprintf('conservation accuracy=%e\n',conservation)
  end
  pijsym=diag(vol)*f ; reciprocity=max(max(abs(pijsym-pijsym'))) ;
  fprintf('reciprocity accuracy=%e\n',reciprocity)