function [rhs, loc, bary, optode] = rbfemrhs(cfg, sd, wv, md)
%
% newcfg=rbmeshprep(cfg)
%
% Create the right-hand-sides for the FEM system equation, here we solve
% forward systems for both source and detector locations in order to use
% the adjoint method to create Jacobians.
%
% author: Qianqian Fang (q.fang <at> neu.edu)
%
% input:
%     cfg: the initial simulation data structure
%
% output:
%     rhs: the right-hand-side of the FEM equation; if multiple sources are
%          used, the rhs is a matrix of dimension Nn x (Ns+Nd), where Nn is
%          the total number of nodes, and Ns the number of sources and
%          Nd is the total number of detectors
%     loc: the indices of the forward mesh element that encloses each
%          source or detector; Nan means the source is outside of the mesh
%          or a wide-field source/detector
%     bary: the barycentric coordinates of the source/detector if it is
%          enclosed by a tetrahedral element; Nan if outside of the mesh
%     optode: the full source+detector position list returned by
%          rbgetoptodes.m
%
% license:
%     GPL version 3, see LICENSE_GPLv3.txt files for details
%
% -- this function is part of Redbird-m toolbox
%

if (nargin > 2)
    sd = sd(wv);
    if (nargin > 3 && size(sd, 2) > 3)
        sd = sd((sd(:, 4) == md | sd(:, 4) == 3), :);
        if (isfield(cfg, 'rhsweight') && isstruct(cfg.rhsweight))
            rhsweight = cfg.rhsweight(md).weight;
        end
    end
    [optsrc, optdet, widesrc, widedet] = rbgetoptodes(cfg, wv);
else
    [optsrc, optdet, widesrc, widedet] = rbgetoptodes(cfg);
end

if exist('sd', 'var')
    if (size(sd, 2) > 3)
        tempsd = sd(sd(:, 3) == 1 & ismember(sd(:, 4), [1 2 3]), :);
    else
        tempsd = sd(sd(:, 3) == 1, :);
    end
    activeOpt = unique(tempsd(:, 1:2));
else
    activeOpt = 1:size([optsrc; optdet], 1) + size([widesrc; widedet], 1);
end

if ((size(optsrc, 1) < 1 && size(widesrc, 1) < 1) || (size(optdet, 1) < 1 && size(widedet, 1) < 1))
    error('you must provide at least one source or detector');
end

loc = [];
bary = [];

rhs = sparse(size(cfg.node, 1), size([widesrc; widedet], 1) + size([optsrc; optdet], 1));
if ~isempty(optsrc)
    [locsrc, barysrc] = tsearchn(cfg.node, cfg.elem, optsrc);
    rhsoptsrc = point2rhs(optsrc, locsrc, barysrc, cfg);
else
    locsrc = [];
    barysrc = [];
    rhsoptsrc = [];
end
if ~isempty(optdet)
    [locdet, barydet] = tsearchn(cfg.node, cfg.elem, optdet);
    rhsoptdet = point2rhs(optdet, locdet, barydet, cfg);
else
    locdet = [];
    barydet = [];
    rhsoptdet = [];
end

loc = [locsrc; nan * ones(size(widesrc, 1), 1); locdet; nan * ones(size(widedet, 1), 1)];
bary = [barysrc; nan * ones(size(widesrc, 1), 4); barydet; nan * ones(size(widesrc, 1), 4)];

rhs = [rhsoptsrc widesrc' rhsoptdet widedet'];

inactive = setdiff(1:size(rhs, 2), activeOpt);

%% Test rhs weight
if (isfield(cfg, 'rhsweight') && ~exist('rhsweight', 'var'))
    rhs = rhs .* cfg.rhsweight.';
elseif exist('rhsweight', 'var')
    rhs = rhs .* rhsweight.';
end
%%

rhs(:, inactive) = 0;

function optoderhs = point2rhs(srcs, locidx, baryc, cfg)
if (size(srcs, 1) ~= size(baryc, 1))
    error('barycentric indices should equal number of sources or detectors');
end

optoderhs = zeros(size(cfg.node, 1), size(srcs, 1));
for ii = 1:size(srcs, 1)
    if ~isnan(locidx(ii))
        optoderhs(cfg.elem(locidx(ii), :), ii) = baryc(ii, :)';
    end
end
