function phi=sn_slab(phi,q,vol,sigt,sigw,ngauss,bc)
% single SN scattering iteration with diamond differencing scheme
% function phi=sn_slab(phi,q,vol,sigt,sigw,ngauss,bc)
% (c) 2010 Alain Hebert, Ecole Polytechnique de Montreal
  nsct=size(sigw,1) ; nreg=size(vol,2) ; nsurf=2 ;
  %----
  %  set the volumic (qext) sources
  %----
  qext=zeros(nsct,nreg) ;
  iof=nsurf ;
  for ii3=1:nreg
    for il=0:nsct-1
      iof=iof+1 ;
      qext(il+1,ii3)=q(iof)+(2*il+1)*phi(iof)*sigw(il+1,ii3) ;
    end
  end
  %----
  %  set Gauss-Legendre angles.
  %----
  [u, w]=lgwt(ngauss, -1.0, 1.0) ; uu=u(ngauss:-1:1) ; u=uu ;
  pl=zeros(nsct,ngauss) ;
  for iangle=1:ngauss
    for il=1:nsct
      if il==1
        pl(il,iangle)=1 ;
      elseif il==2
        pl(il,iangle)=u(iangle) ;
      elseif il==3
        pl(il,iangle)=0.5*(3.0*u(iangle)^2-1.0) ;
      else
        error('nsct too big')
      end
    end
  end
  %----
  %  outer loop over mu levels.
  %----
  flux=zeros(nsct,nreg) ; curr=zeros(2,1) ;
  for ip=1:ngauss/2
    afbp=q(1) ; afbm=q(2) ;
    switch bc
      case {'void'}
      case {'tran'}
        %----
        %  shooting method.
        %----
        afba1=0 ; afba2=0 ; afbb1=1 ; afbb2=1 ;
        for i=nreg:-1:1
          qel=sum(qext(:,i).*pl(:,ip))/2 ;
          qsola=(qel*vol(i)-2*u(ip)*afba1)/(sigt(i)*vol(i)-2*u(ip));
          qsolb=(qel*vol(i)-2*u(ip)*afbb1)/(sigt(i)*vol(i)-2*u(ip));
          afba1=2*qsola-afba1 ; afbb1=2*qsolb-afbb1 ;
        end
        for i=1:nreg
          qel=sum(qext(:,i).*pl(:,ngauss-ip+1))/2 ;
          qsola=(qel*vol(i)+2*u(ngauss-ip+1)*afba2)/(sigt(i)*vol(i)+2*u(ngauss-ip+1));
          qsolb=(qel*vol(i)+2*u(ngauss-ip+1)*afbb2)/(sigt(i)*vol(i)+2*u(ngauss-ip+1));
          afba2=2*qsola-afba2 ; afbb2=2*qsolb-afbb2 ;
        end
        afbm=afba1/(1+afba1-afbb1) ; afbp=afba2/(1+afba2-afbb2) ;
      otherwise
        error('unknown type of boundary condition')
    end
    for i=nreg:-1:1
      qel=sum(qext(:,i).*pl(:,ip))/2 ;
      qsol=(qel*vol(i)-2*u(ip)*afbm)/(sigt(i)*vol(i)-2*u(ip));
      afbm=2*qsol-afbm ;
      flux(:,i)=flux(:,i)+w(ip)*qsol*pl(:,ip) ;
    end
    curr(1)=curr(1)+w(ip)*u(ip)*afbm ;
    for i=1:nreg
      qel=sum(qext(:,i).*pl(:,ngauss-ip+1))/2 ;
      qsol=(qel*vol(i)+2*u(ngauss-ip+1)*afbp)/(sigt(i)*vol(i)+2*u(ngauss-ip+1));
      afbp=2*qsol-afbp ;
      flux(:,i)=flux(:,i)+w(ngauss-ip+1)*qsol*pl(:,ngauss-ip+1) ;
    end
    curr(2)=curr(2)+w(ngauss-ip+1)*u(ngauss-ip+1)*afbp ;
  end
  %----
  %  reformat the flux
  %----
  phi(1:nsurf)=curr ;
  iof=nsurf ;
  for ii3=1:nreg
    for ii1=1:nsct
      iof=iof+1 ;
      phi(iof)=flux(ii1,ii3) ;
    end
  end