function [keff,rekeff]=mcflx(myNode,bc,alb,k_cycle,nsrck,ikz,kct,xyzl,varargin)
% Power iteration with the Monte-Carlo method in 1D/2D/3D Cartesian geometry.
% function [keff,rekeff]=mcflx(myNode,bc,alb,nsrck,ikz,kct,xyzl,varargin)
% (c) 2008 Alain Hebert, Ecole Polytechnique de Montreal
  ndim=size(xyzl,2) ; iseed=1111 ; ibank1=nsrck ; nu_cycle=zeros(1,2*nsrck) ;
  ngen1=repmat(struct('mix',0,'isonbr',-1,'nu',1.0,'pos',zeros(1,ndim)),1,2*nsrck) ;
  ngen2=ngen1 ;
  %----
  %  initial uniform source
  %----
  for itrk=1:nsrck
    for idim=1:ndim
      [iseed,rand]=randf(iseed) ;
      ngen1(itrk).pos(idim)=xyzl(1,idim)+rand*(xyzl(2,idim)-xyzl(1,idim)) ;
    end
  end
  %----
  %  power iteration
  %----
  sum1=0.0 ; sum2=0.0 ;
  for icycle=1:kct
    ibank2=0 ;
    for itrk=1:ibank1
      nu=ngen1(itrk).nu ; nloop=max(1,fix(nu/k_cycle)) ;
      ngen1(itrk).nu=nu/(k_cycle*nloop) ;
      for iloop=1:nloop
        [ngen iseed]=mctrk(myNode,bc,alb,ngen1(itrk),iseed,xyzl,varargin{:}) ;
        if ngen.isonbr > 0
          ibank2 = ibank2+1 ;
          if ibank2 > 2*nsrck, error('too many neutron tracks being banked'), end
          ngen2(ibank2)=ngen ; nu_cycle(ibank2)=ngen.nu ;
        end
      end
    end
    k_cycle=sum(cat(1,ngen2(1:ibank2).nu))/nsrck ;
    %----
    %  Russian roulette
    %----
    if ibank2 > nsrck
      nu_cycle=sort(nu_cycle,'descend') ; nulimit=nu_cycle(nsrck) ; ibank1=0 ;
      for itrk=1:ibank2
        lkeep=(ngen2(itrk).nu >= nulimit) ;
        if ~lkeep
          [iseed,rand]=randf(iseed) ; lkeep=(rand <= ngen2(itrk).nu/nulimit) ;
          ngen2(itrk).nu = nulimit ;
        end
        if lkeep
          ibank1=ibank1+1 ; ngen1(ibank1)=ngen2(itrk) ;
        end
      end
    else
      ibank1=ibank2 ; ngen1=ngen2 ;
    end
    %----
    %  k-effective of the problem with the relative error
    %----
    if icycle > ikz
      sum1=sum1+k_cycle ; sum2=sum2+k_cycle^2 ;
    end
    fprintf(1,'cycle number=%d k-eff cycle=%e\n',icycle,k_cycle)
  end
  keff=sum1/(kct-ikz) ; rekeff=sqrt((sum2-(kct-ikz)*keff^2)/(kct-ikz-1)/(kct-ikz)) ;
  fprintf(1,'cycle number=%d k-eff cycle=%e k-eff expected= %e sigma=%e\n', ...
  kct,k_cycle,keff,rekeff)