function Jmua = rbjacmuafast(sd, phi, nvol, elem)
%
% Jmua=rbjacmuafast(sd, phi, nvol, elem)
%
% computing the nodal Jacobian matrix for mua using an approximated
% method as shownn in Qianqian's Thesis (nodal-adjoint formulation)
%
% author: Qianqian Fang (q.fang <at> neu.edu)
%
% input:
%     sd: the source-detector mapping table
%     phi: the forward solutions at all sources
%     nvol: the nodal volume computed by nodevolume from iso2mesh
%     elem: the mesh element list
%
% output:
%     Jmua: the approximated Jacobian for mua using nodal-adjoint formulation
%
% license:
%     GPL version 3, see LICENSE_GPLv3.txt files for details
%
% -- this function is part of Redbird-m toolbox
%

if (nargin < 3 || isempty(sd) || isempty(phi) || isempty(nvol))
    error('you must give at least the first 3 inputs and they must not be empty');
end

wavelengths = {''};

if (isa(phi, 'containers.Map'))
    wavelengths = phi.keys;
else
    phi = containers.Map({''}, {phi});
end

Jmua = containers.Map();

for waveid = wavelengths
    wv = waveid{1};
    phiwv = phi(wv);
    if (isa(sd, 'containers.Map'))
        sdwv = sd(wv);
    else
        sdwv = sd;
    end
    if (size(phiwv, 1) == length(nvol))
        Ja = zeros(size(sdwv, 1), size(phiwv, 1));
        for i = 1:size(sdwv, 1)
            Ja(i, :) = phiwv(:, sdwv(i, 1)) .* phiwv(:, sdwv(i, 2)) .* nvol(:);
        end
    elseif (nargin > 3 && size(nvol, 1) == size(elem, 1))
        Ja = zeros(size(elem, 1), size(sdwv, 1));
        for i = 1:size(sdwv, 1)
            for j = 1:4
                Ja(:, i) = Ja(:, i) + phiwv(elem(:, j), sdwv(i, 1)) .* phiwv(elem(:, j), sdwv(i, 2)) .* nvol(:);
            end
        end
        Ja = Ja.' * 0.25;
    else
        error('the row number of phi must be the same as the length of nvol or elem is needed');
    end
    Jmua(wv) = -Ja; % increasing mua, decreasing phi
end

% if only a single wavelength is required, return regular arrays instead of a map
if (length(wavelengths) == 1)
    Jmua = Jmua(wavelengths{1});
end
