%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MCXLAB - Monte Carlo eXtreme for MATLAB/Octave by Qianqina Fang
%
% In this example, we validate surface diffuse reflectance extracted from MC
% detected photon profiles against analytical solutions computed from
% diffusion equations (semi-infinite domain).
%
% [Yao2018] Ruoyang Yao, Xavier Intes, and Qianqian Fang, "Direct approach
% to compute Jacobians for diffuse optical tomography using perturbation Monte
% Carlo-based photon 'replay'," Biomed. Opt. Express 9, 4588-4603 (2018)
%
% [Kienle1997] Alwin Kienle and Michael S. Patterson, "Improved solutions
% of the steady-state and the time-resolved diffusion equations for reflectance
% from a semi-infinite turbid medium," J. Opt. Soc. Am. A 14, 246-254 (1997)
%
% This file is part of Monte Carlo eXtreme (MCX) URL:https://mcx.space
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% only clear cfg to avoid accidentally clearing other useful data
clear cfg;

%% Monte Carlo
cfg.nphoton = 1e8; % photon number
cfg.vol = uint8(ones(60, 60, 60)); % size in mm (length,width,height)

cfg.prop = [0.00  0.0  1.0 1.0   % first row: ambient medium by default
            0.01  10.0 0.9 1.4]; % media optical properties

cfg.srcpos = [30 30 0]; % (x,y,z) co-ordinates of light source
cfg.srcdir = [0 0 1];   % Direction - positive z-direction
cfg.issrcfrom0 = 1;

cfg.gpuid = 1;
cfg.autopilot = 1;

% time window: for time-solved results
cfg.tstart = 0;
cfg.tend = 5e-9;
cfg.tstep = 5e-9;

% save positions & directions of escaping photons
cfg.issaveexit = 1;

% set up detectors: non-overlapping fiber detectors along x = 30
detpy = 31:1:55; % y changes from 31 to 55
r = 0.25; % radius 0.25mm
for i = 1:size(detpy, 2)
    cfg.detpos(i, :) = [30, detpy(i), 0, r];
end

% visualize domain settings
figure(1);
mcxpreview(cfg);

% launch simulations
[fluence, detp] = mcxlab(cfg);

% compute diffuse reflectance at detectors
drefmc = mcxcwdref(detp, cfg);

%% Diffusion
% flux at detectors
fluxdiffusion = cwfluxdiffusion(cfg.prop(2, 1), cfg.prop(2, 2) * (1 - cfg.prop(2, 3)), ...
                                rbgetreff(1.4, 1.0), cfg.srcpos, cfg.detpos(:, 1:3));

% fluence at detectors
fluencediffusion = cwfluencediffusion(cfg.prop(2, 1), cfg.prop(2, 2) * (1 - cfg.prop(2, 3)), ...
                                      rbgetreff(1.4, 1.0), cfg.srcpos, cfg.detpos(:, 1:3));

% surface diffuse reflectance
drefdiffusion = 0.118 * fluencediffusion + 0.306 * fluxdiffusion; % Eq. 8 of Kienle1997

%% compare diffuse reflectance
figure(2);
h = semilogy(1:length(detpy), drefdiffusion, 'r--', 1:length(detpy), drefmc, 'b+');
legend("Diffusion", "MC");
title("log10(Diffuse reflectance) [1/mm^{2}]");
xlabel("Source-Detector seperation [mm]");
