dorefine2d(ibox,select,firstgenresidbox)

function [ibox,newidx] = dorefine2d(ibox,select,firstgenresidbox)
% Box refinement process

N = length(ibox);  N0 = N;
for i = find(select(:)).'
  % Coords of proposed box midpoints.
  dx = 0.5*ibox(i).dimen;

  if ibox(i).isboundary == false
    p = repmat(ibox(i).midpt,4,1) + 0.5*[1 1;-1 1;-1 -1;1 -1]*diag(dx);
  elseif ibox(i).isboundary == true
    dimen = ibox(i).dimen;
    if dimen(2) == 0
        p = repmat(ibox(i).midpt,2,1) + 0.5*[1 1;-1 1]*diag(dx);
    elseif dimen(1) == 0
        p = repmat(ibox(i).midpt,2,1) + 0.5*[1 1;1 -1]*diag(dx);
    end
  end

  n = size(p,1);
  p = num2cell(p,2);
  [ibox(N+(1:n)).midpt] = deal(p{:});
  [ibox(N+(1:n)).dimen] = deal(dx);
  [ibox(N+(1:n)).isboundary] = deal(ibox(i).isboundary);
  [ibox(N+(1:n)).child] = deal([]);
  [ibox(N+(1:n)).Nchild] = deal(0);
  [ibox(N+(1:n)).depth] = deal(1+ibox(i).depth);

  ibox(i).child = N+(1:n);
  ibox(i).Nchild = n;
  switch firstgenresidbox
      case 0
        ibox(i).depth = 2+ibox(i).depth;
  end
  N = N+n;
end
newidx = (N0+1:length(ibox))';