Failed to apply double-float-buffer approach in #41 to `dref` cause loss of accuracy

[fangq]

This issue was brought up by Haohui Zhang in this mailing list post

https://groups.google.com/g/mcx-users/c/bbr2Yx_Qvws

running the below script returned conserved energy with absorb + reflect ~= 1 when nphoton is 1e6/1e7 or lower, but became less than 1 when nphoton increases to 1e8 (sum=0.85) or 1e9

% Bulid the shape
Vfinal = zeros(302,302,152);
Vfinal(2:301,2:301,2:151) = 1;
% Vfinal(:,:,end) = 0;
%% prepare cfg for MCX simulation
clear cfg
cfg.nphoton=1e8;
cfg.outputtype='energy';
cfg.seed = 77542;
% tissue labels:0-ambient air,1-tissue
cfg.vol=Vfinal;
cfg.prop = [0,0,1,1;0.01,2,0,1];
% light source
cfg.srctype = 'pencil';
cfg.srcpos=[149.5,149.5,0];
cfg.srcdir=[0,0,1]; 
cfg.issrcfrom0=1;
% time windows
cfg.tstart=0;
cfg.tend=1e-8;
cfg.tstep=1e-9;
% other simulation parameters
cfg.isspecular=1;
cfg.isreflect=1;
cfg.autopilot=1;
cfg.gpuid=1;
cfg.issaveref=1; % save diffuse reflectance
cfg.unitinmm = 1;
% Detector
cfg.detpos=[50,50,100,2];
cfg.debuglevel='P';
tic
%% run MCX simulation
[flux,detps]=mcxlab(cfg);
toc
%% post-simulation data processing and visualization
% convert time-resolved fluence to CW fluence
CWfluence=sum(flux.data,4); 
cwdref=sum(flux.dref,4);     % diffuse reflectance
absorb = sum(CWfluence,'all');
reflect = sum(cwdref(:,:,:),'all');
absorb + reflect

This issue was caused by the same floating-point accumulation error as described in #41 - i.e. adding a small energy-loss to a large accumulated single-precision floating point number cause a significant loss of accuracy.

In #41, this issue was solved by using a double float-buffer approach, but we only applied it for fluence accumulation. For dref, it still uses the single floating-point buffer for accumulation, thus cause the reflect value became lower than expected when running large number of photons.

Fixing this requiring applying the same strategy as we did in #41 but for dref.