function [Amat, deldotdel] = rbfemlhs(cfg, deldotdel, wavelength, mode)
%
% [Amat,deldotdel]=rbfemlhs(cfg)
%   or
% [Amat,deldotdel]=rbfemlhs(cfg, wavelength)
% [Amat,deldotdel]=rbfemlhs(cfg, deldotdel, wavelength)
%
% create the FEM stiffness matrix (left-hand-side) for solving the
% diffusion equation at a given wavelength
%
% author: Qianqian Fang (q.fang <at> neu.edu)
%
% input:
%     cfg: the initial simulation data structure
%     deldotdel (optional): precomputed operator on the mesh (del_phi dot del_phi)
%         where del represents the gradient; see help rbdeldotdel
%     wavelength (optional): a string or number denoting the wavelength
%
% output:
%     Amat: the left-hand-side matrix of the FEM equation - a sparse matrix
%          of dimension Nn x Nn, where Nn is the number of nodes of the
%          forward mesh
%     deldotdel: if the 2nd input is not given, this function compute
%          deldotdel and return as the 2nd output
%
% license:
%     GPL version 3, see LICENSE_GPLv3.txt files for details
%
% -- this function is part of Redbird-m toolbox
%

nn = size(cfg.node, 1);
ne = size(cfg.elem, 1);

R_C0 = (1 ./ 299792458000.);

% cfg.prop is updated from cfg.param and contains the updated mua musp.
% if cfg.param is node/elem based, cfg.prop is updated to have 4 columns
% with mua/musp being the first two columns
% if cfg.param is segmentation based, cfg.prop has the same format as mcx's
% prop, where the first row is label 0, and total length is Nseg+1

prop = cfg.prop;
cfgreff = cfg.reff;
if exist('mode', 'var')
    if mode == 1
        omega = cfg.omega;
    elseif mode == 2
        omega = 0;
    end
else
    omega = cfg.omega;
end

if (nargin == 2 && numel(deldotdel) == 1)
    wavelength = deldotdel;
end

if (isa(cfg.prop, 'containers.Map')) % if multiple wavelengths, take current
    if (nargin < 3)
        error('you must specify wavelength');
    end
    if (~ischar(wavelength))
        wavelength = sprintf('%g', wavelength);
    end
    prop = cfg.prop(wavelength);
    cfgreff = cfg.reff(wavelength);
    if (isa(omega, 'containers.Map'))
        omega = omega(wavelength);
    end
end

% if deldotdel is provided, call native code; otherwise, call mex

if (nargin >= 2 && numel(deldotdel) > 1)
    % get mua from prop(:,1) if cfg.prop has wavelengths
    if (size(prop, 1) == nn || size(prop, 1) == ne)
        mua = prop(:, 1);
        if (size(prop, 2) < 3)
            musp = prop(:, 2);
        else
            musp = prop(:, 2) .* (1 - prop(:, 3));
        end
    elseif (size(prop, 1) < min([nn ne])) % use segmentation based prop list
        mua = prop(cfg.seg + 1, 1);
        if (size(prop, 2) < 3)
            musp = prop(cfg.seg + 1, 2); % assume g is 0
        else
            musp = prop(cfg.seg + 1, 2) .* (1 - prop(cfg.seg + 1, 3));
        end
    end
    dcoeff = 1 ./ (3 * (mua + musp));

    if (isfield(cfg, 'bulk') && isfield(cfg.bulk, 'n'))
        nref = cfg.bulk.n;
    elseif (isfield(cfg, 'seg') && size(prop, 1) < min([nn, ne]))
        nref = rbgetbulk(cfg);
        if (isa(nref, 'containers.Map'))
            nref = nref(wavelength);
        end
        nref = nref(4);
    else
        nref = prop(:, 4);
    end
    Reff = cfgreff;

    edges = sort(meshedge(cfg.elem), 2);

    % what LHS matrix needs is dcoeff and mua, must be node or elem long
    if (length(mua) == size(cfg.elem, 1))  % element based property
        Aoffd = deldotdel(:, [2:4, 6:7, 9]) .* repmat(dcoeff(:), 1, 6) + repmat(0.05 * mua(:) .* cfg.evol(:), 1, 6);
        Adiag = deldotdel(:, [1, 5, 8, 10]) .* repmat(dcoeff(:), 1, 4) + repmat(0.10 * mua(:) .* cfg.evol(:), 1, 4);
        if (omega > 0)
            Aoffd = complex(Aoffd, repmat(0.05 * omega * R_C0 * nref(:) .* cfg.evol(:), 1, 6));
            Adiag = complex(Adiag, repmat(0.10 * omega * R_C0 * nref(:) .* cfg.evol(:), 1, 4));
        end
    else  % node based properties
        w1 = (1 / 120) * [2 2 1 1; 2 1 2 1; 2 1 1 2; 1 2 2 1; 1 2 1 2; 1 1 2 2]';
        w2 = (1 / 60) * (diag([2 2 2 2]) + 1);
        mua_e = reshape(mua(cfg.elem), size(cfg.elem));
        if (length(nref) == 1)
            nref_e = nref * ones(size(cfg.elem));
        else
            nref_e = reshape(nref(cfg.elem), size(cfg.elem));
        end
        dcoeff_e = mean(reshape(dcoeff(cfg.elem), size(cfg.elem)), 2);
        Aoffd = deldotdel(:, [2:4, 6:7, 9]) .* repmat(dcoeff_e, 1, 6) + (mua_e * w1) .* repmat(cfg.evol(:), 1, 6);
        Adiag = deldotdel(:, [1, 5, 8, 10]) .* repmat(dcoeff_e, 1, 4) + (mua_e * w2) .* repmat(cfg.evol(:), 1, 4);
        if (omega > 0)
            Aoffd = complex(Aoffd, (omega * R_C0) * (nref_e * w1) .* repmat(cfg.evol(:), 1, 6));
            Adiag = complex(Adiag, (omega * R_C0) * (nref_e * w2) .* repmat(cfg.evol(:), 1, 4));
        end
    end
    % add partial current boundary condition
    edgebc = sort(meshedge(cfg.face), 2);
    Adiagbc = cfg.area(:) * ((1 - Reff) / (12 * (1 + Reff)));
    Adiagbc = repmat(Adiagbc, 1, 3);
    Aoffdbc = Adiagbc * 0.5;

    Amat = sparse([edges(:, 1); edges(:, 2); cfg.elem(:); edgebc(:, 1); edgebc(:, 2); cfg.face(:)], ...
                  [edges(:, 2); edges(:, 1); cfg.elem(:); edgebc(:, 2); edgebc(:, 1); cfg.face(:)], ...
                  [Aoffd(:); Aoffd(:); Adiag(:); Aoffdbc(:); Aoffdbc(:); Adiagbc(:)]);
else
    if (size(cfg.elem, 2) > 4)
        cfg.elem(:, 5:end) = [];
    end
    cfg.prop = prop; % use property of the current wavelength
    [Adiag, Aoffd, deldotdel] = rbfemmatrix(cfg);
    Amat = sparse([cfg.rows, cfg.cols, (1:nn)], [cfg.cols, cfg.rows, (1:nn)], [Aoffd, Aoffd, Adiag], nn, nn);
    deldotdel = deldotdel';
end
