function [Jmua_n, Jmua_e, Jd_n, Jd_e]=rbjac(sd, phi, deldotdel, felem, evol, varargin)
%
% [Jmua_n, Jmua_e, Jd_n, Jd_e]=rbjac(sd, phi, deldotdel, felem, evol, isnodal)
%
% Building the Jacobian matrices using adjoint method and native matlab code
%
% author: Qianqian Fang (q.fang <at> neu.edu)
%
% input:
%     sd: the source-detector mapping table
%     phi: the forward solutions at all sources
%     deldotdel: grad*phi dot product with grad phi, computed as part of the computation
%     felem: forward mesh element list
%     evol: forward mesh element volume
%     iselem (optional): 1 if elem-based Jacobian; 0 (default) if is nodal-based 
%
% output:
%     Jmua_n: the nodal Jacobian for absorption coeff. mua
%     Jmua_e: the element-wise Jacobian for absorption coeff. mua
%     Jd_n: (optional) the nodal Jacobian for diffusion coeff D
%     Jd_e: (optional) the element-wise Jacobian for diffusion coeff D
%
% license:
%     GPL version 3, see LICENSE_GPLv3.txt files for details 
%
% -- this function is part of Redbird-m toolbox
%

if(nargin<5 || isempty(sd) || isempty(phi) || isempty(deldotdel)|| isempty(evol))
    error('you must give at least 5 inputs and they must not be empty');
end
if (isempty(varargin))
    iselem=0;
    rfcw = 1;
elseif (length(varargin) == 1)
    iselem = varargin{1};
    rfcw = 1;
else
    opt=varargin2struct(varargin{:});
    rfcw = jsonopt('rfcw',1,opt);
    iselem = jsonopt('iselem',0,opt);
end
nelem=size(felem,1);

wavelengths={''};

if isstruct(phi)
    if(isa(phi(1).phi,'containers.Map'))
        wavelengths = phi(1).phi.keys;
    else
        for ii = 1:length(phi)
            phi(ii).phi = containers.Map({''},{phi(ii).phi});
        end
    end        
else
    if(isa(phi,'containers.Map'))
        wavelengths=phi.keys;
        phitemp = phi; clear phi
        phi(rfcw).phi = phitemp;
        clear phitemp
    else
        phitemp = phi; clear phi
        phi(rfcw).phi = containers.Map({''},{phitemp});
        clear phitemp
    end
end
mode = 1:length(phi);

for ii = mode
    Jmua_n(ii).J = containers.Map();
    Jmua_e(ii).J = containers.Map();
    if(nargout>2)
        Jd_n(ii).J=containers.Map();
        Jd_e(ii).J=containers.Map();
    end
end

idx=[1 1 1 2 2 3 2 3 4 3 4 4
     2 3 4 3 4 4 1 1 1 2 2 3
     2 3 4 6 7 9 2 3 4 6 7 9];

for waveid=wavelengths
    wv=waveid{1};
    if(isa(sd,'containers.Map'))
        sdwv=sd(wv);
    else
        sdwv=sd;
    end
    for md = rfcw
        if (size(sdwv,2) < 4)
            sdwv(:,4) = md;
        end
        sdmd = sdwv(((sdwv(:,3) == 1) & (sdwv(:,4) == md | sdwv(:,4) == 3)),:);
        phiwv = phi(md).phi(wv);
        Jmua_node=zeros(size(sdmd,1),size(phiwv,1));
        Jmua_elem=zeros(size(sdmd,1),nelem);
        if(nargout>2)
            Jd_node=zeros(size(sdmd,1),size(phiwv,1));
            Jd_elem=zeros(size(sdmd,1),nelem);
        end
        for i=1:nelem
            phidotphi1=phiwv(felem(i,1:4),sdmd(:,1)).*phiwv(felem(i,1:4),sdmd(:,2));
            phidotphi2=phiwv(felem(i,idx(1,:)),sdmd(:,1)).*phiwv(felem(i,idx(2,:)),sdmd(:,2));  
            Jmua_elem(:,i)=-(sum(phidotphi1,1)+sum(phidotphi2,1)*0.5)*(0.1*evol(i));
            if(nargout>2)
                Jcol=phidotphi1.'*deldotdel(i,[1 5 8 10]).';
                Jcol=Jcol+phidotphi2.'*deldotdel(i,idx(3,:)).';
                Jd_elem(:,i)=-Jcol;
            end
            for j=1:4
                Jmua_node(:,felem(i,j))=Jmua_node(:,felem(i,j))+Jmua_elem(:,i);
            end
            if(nargout>2)
                for j=1:4
                    Jd_node(:,felem(i,j))=Jd_node(:,felem(i,j))-Jcol;
                end
            end
        end
        Jmua_node=Jmua_node*0.25;
        Jmua_n(md).J(wv)=Jmua_node;
        Jmua_e(md).J(wv)=Jmua_elem;
        if(nargout>2)
            Jd_node=Jd_node*0.25;
            Jd_n(md).J(wv)=Jd_node;
            Jd_e(md).J(wv)=Jd_elem;
        end
    end
end

% if only a single wavelength is required, return regular arrays instead of a map
if(length(wavelengths)==1)
    for md = mode
        Jmua_n(md).J=Jmua_n(md).J(wavelengths{1});
        Jmua_e(md).J=Jmua_e(md).J(wavelengths{1});
        if(nargout>2)
            Jd_n(md).J=Jd_n(md).J(wavelengths{1});
            Jd_e(md).J=Jd_e(md).J(wavelengths{1});
        end
    end
end

if (length(rfcw) == 1)      %%(length(mode) == 1)
    Jmua_n = Jmua_n(rfcw).J;    %%Jmua_n(1).J;
    Jmua_e = Jmua_e(rfcw).J;    %%Jmua_e(1).J;
    if (nargout>2)
        Jd_n = Jd_n(rfcw).J;    %%Jd_n(1).J;
        Jd_e = Jd_e(rfcw).J;    %%Jd_e(1).J;
    end
end

if(iselem && nargout<3)
    Jmua_n=Jmua_e;
    if nargout>2
        Jmua_e=Jd_e;
    end
end