function [ngen,iseed]=mctrk(myNode,bc,alb,ngen,iseed,xyzl,varargin)
% Random walk of a single particle from source to death.
% function [ngen,iseed]=mctrk(myNode,bc,alb,ngen,iseed,xyzl,varargin)
% (c) 2008 Alain Hebert, Ecole Polytechnique de Montreal
% based on a Fortran program written by Benoit Arsenault
  length=0.0 ; active=true ; idso=0 ;
  nreg=size(myNode,2) ; ng=size(myNode(1).sigt,2) ; nfiss=size(myNode(1).nusigf,1) ;
  %----
  % set the initial energy and direction of the particle
  %----
  if ngen.mix == 0
    ig=1 ;
  else
    [iseed,rand]=randf(iseed) ; xprob=0.0 ; ig=0 ;
    while (rand >= xprob) && (ig < ng)
      ig=ig+1 ;
      xprob=xprob+myNode(ngen.mix).chi(ngen.isonbr,ig) ;
    end
  end
  xsm=0.0 ;
  for ireg=1:nreg
    xsm=max(xsm,myNode(ireg).sigt(ig)) ;
  end
  [iseed,rand]=randf(iseed) ; mu=2.0*rand-1.0 ; mus=sqrt(1.0-mu*mu) ;
  [iseed,rand]=randf(iseed) ; phi=2.0*pi*rand ;
  vdir=[mus*cos(phi),mus*sin(phi),mu] ;
  %----
  % random walk of the particle
  %----
  while active
    if idso == 0
      [iseed,rand]=randf(iseed) ; length=-log(rand)/xsm ;
    end
    [ngen.pos,ireg,length,idso]=mctptr(length,vdir,ngen.pos,xyzl,varargin{:}) ;
    if ireg > 0
%     detect if an interaction is a virtual or a real collision
      ngen.mix=ireg ;
      [iseed,rand]=randf(iseed) ; virtual=rand <= ((xsm-myNode(ireg).sigt(ig))/xsm) ;
      if ~virtual
%       determine the type of reaction
        [iseed,rand]=randf(iseed) ; sigs0=sum(myNode(ireg).sigs(1:ng,ig)) ;
        pscat=sigs0/myNode(ireg).sigt(ig) ;
        pn2n=myNode(ireg).sn2n(ig)/myNode(ireg).sigt(ig) ;
        if rand <= pscat+pn2n
%         isotropic scattering or (n,2n) event
          if rand > pscat, ngen.nu =2.0*ngen.nu ; end
          [iseed,rand]=randf(iseed) ; xx=0.0 ; igroup=ng ;
          while (rand >= xx) && (igroup >= 1)
            xx=xx+myNode(ireg).sigs(igroup,ig)/sigs0 ;
            igroup=igroup-1 ;
          end
          ig=igroup+1 ; active=true ; xsm=0.0 ;
          for ireg=1:nreg
            xsm=max(xsm,myNode(ireg).sigt(ig)) ;
          end
        else
%         capture or fission event
          active=false ;
          nusigf0=sum(myNode(ireg).nusigf(1:nfiss,ig)) ;
          if nusigf0 > 0.0 && nfiss == 1
            ngen.nu =ngen.nu*myNode(ireg).nusigf(1,ig)/(myNode(ireg).sigt(ig)- ...
                sigs0-myNode(ireg).sn2n(ig)) ;
            ngen.isonbr =1 ;
          elseif nusigf0 > 0.0
            [iseed,rand]=randf(iseed) ; xprob=0.0 ; ngen.isonbr=0 ;
            while (rand >= xprob) && (ngen.isonbr < nfiss)
              ngen.isonbr=ngen.isonbr+1 ;
              xprob=xprob+myNode(ireg).nusigf(ngen.isonbr,ig)/nusigf0 ;
            end
            ngen.nu =ngen.nu*myNode(ireg).nusigf(iso,ig)/(myNode(ireg).sigt(ig)- ...
                sigs0-myNode(ireg).sn2n(ig)) ;
          else
            ngen.nu =0.0 ; ngen.isonbr =0 ;
          end
        end
        if active
          [iseed,rand]=randf(iseed) ; mu=2.0*rand-1.0 ; mus=sqrt(1.0-mu*mu) ;
          [iseed,rand]=randf(iseed) ; phi=2.0*pi*rand ;
          vdir=[mus*cos(phi),mus*sin(phi),mu] ;
        end
      end
    else
      ids=abs(idso) ; ngen.mix=0 ;
      switch bc
        case {'refl'}
%         reflective boundary conditions
          vdir(ids)=-vdir(ids) ;
        case {'void'}
          active=false ; ngen.nu =0.0 ; ngen.isonbr =0 ;
        case {'tran'}
%         periodic boundary conditions
          if idso > 0
            ngen.pos(ids)=ngen.pos(ids)-(xyzl(2,ids)-xyzl(1,ids)) ;
          else
            ngen.pos(ids)=ngen.pos(ids)+(xyzl(2,ids)-xyzl(1,ids)) ;
          end
        case {'albe'}
          if alb == 0.0
            active=false ; ngen.nu =0.0 ; ngen.isonbr =0 ;
          else
            [iseed,rand]=randf(iseed) ; active=(rand <= alb) ;
            if active
%             for isotropic bc, choose randomly the reentering direction
              indos=[ 2 3 1 ; 3 1 2 ] ;
              [iseed,rand]=randf(iseed) ; mu=2.0*rand-1.0 ; mus=sqrt(1.0-mu*mu) ;
              [iseed,rand]=randf(iseed) ; phi=pi*rand ;
              vdir(ids)=-sign(idso)*mus*sin(phi) ;
              vdir(indos(1,ids))=mus*cos(phi) ; vdir(indos(2,ids))=mu ;
              for idim=1:size(xyzl,2)
                if idim ~= ids
                  [iseed,rand]=randf(iseed) ;
                  ngen.pos(idim)=xyzl(1,idim)+rand*(xyzl(2,idim)-xyzl(1,idim)) ;
                end
              end
              idso=0 ;
            else
              ngen.nu =0.0 ; ngen.isonbr =0 ;
            end
          end
        otherwise
          error('unknown type of boundary condition')
      end
    end
  end