if (~exist('rbrun', 'file'))
    addpath(fullfile(pwd, '../matlab'));
end

clear cfg;

[cfg.node, cfg.face, cfg.elem] = meshabox([40 0 0], [160, 120, 60], 10);

nn = size(cfg.node, 1);

[xi, yi] = meshgrid(60:20:140, 20:20:100);
cfg.srcpos = [xi(:), yi(:), zeros(numel(yi), 1)];
cfg.detpos = [xi(:), yi(:), 60 * ones(numel(yi), 1)];
cfg.srcdir = [0 0 1];
cfg.detdir = [0 0 -1];

% label based properties for elements
cfg.prop = [
            0 0 1 1
            0.008 1 0 1.37
           ];
cfg.seg = ones(size(cfg.elem, 1), 1);

cfg.omega = 2 * pi * 70e6;

tic;
cfg = rbmeshprep(cfg);
fprintf('preparing mesh ... \t%f seconds\n', toc);

rbgetbulk(cfg);

%% RF forward using segmentation based properties

[detphi1, phi1] = rbrunforward(cfg);

%%  node-based properties

% cfg.seg=ones(size(cfg.node,1),1);
% cfg.prop=[
%     0 0 1 1
%     0.008 1 0 1.37
%     0.016 1 0 1.37
% ];
cfg.prop = repmat([0.008 1 0 1.37], size(cfg.node, 1), 1);

rbgetbulk(cfg);

[detphi2, phi2] = rbrunforward(cfg);

%% comparing the two solutions

rr = detphi1 ./ detphi2;
fprintf(1, 'max ratio=%f , min ratio=%f\n', max(abs(rr(:))), min(abs(rr(:))));
fprintf(1, 'max angle diff=%f , min angle diff=%f\n', max(angle(rr(:))), min(angle(rr(:))));
