% Banded radiosity equation solver 
function [qtot_out, q_out, J_out] = SolveBandRadEqs(Tin,F,e,varargin)
  N=size(F,2);                              % number of surfaces
  bands=[0 Inf];                            % default bands
  out=reshape([1:N],size(Tin));             % default out
  T=Tin(:);                                 % make T a column vector
  for i=1:length(varargin)                  % modify defaults
    if varargin{i}(end) == Inf
      bands=varargin{i};
    else
      out=varargin{i};
    end
  end
  B=size(e,2);                              % number of bands
  if B ~= length(bands)-1
    error('SolveBandRadEqs: size(e,2) ~= bands');
  end
  if N ~= size(e,1) || N ~= size(F,1) || N ~= length(T)
    error('SolveBandRadEqs: size(F) ~= [length(T) size(e,1)]');
  end
  J=zeros(N,B);                             % banded radiosities
  q=zeros(N,B);                             % banded heat fluxes
  T=T(:);                                   % T needs to be a column vector
  Fb=FlTb(bands,T);                         % banded fractional functions
  for b=1:B,                                % loop of all the bands
    gray= find(e(:,b) < 1);                 % gray band surface indices
    black=find(e(:,b) == 1);                % black band surface indices
    Eb=Eblk(T).*Fb(:,b);                    % blackbody emission in band
    re=e(gray,b)./(1-e(gray,b));            % e/(1-e) for gray band
    C=diag(1+re)-F(gray,gray);              % coefficient matrix
    D=re.*Eb(gray)+F(gray,black)*Eb(black); % resultant vector
    Jg=C\D;                                 % solve for unknown radiosities
    J(gray,b) = Jg;                         % gray band radiosities
    J(black,b) = Eb(black);                 % black band radiosities
    q(:,b) = J(:,b)-F*J(:,b);               % banded heat fluxes
  end
  J_out=J(out(:),:);                        % select requested output
  q_out=q(out(:),:);                        % select requested output
  qtot_out=sum(q_out,2);                    % total heat fluxes
  qtot_out=reshape(qtot_out,size(out));     % match format 
end
