%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Demo 3: detector configuration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear cfg
close all

% Comment the below line if you have an NVIDIA GPU
% USE_MCXCL=1;
rand_seed = randi([1 2^31-1], 1, 1);
% Create matrix that saves the data for a specific photon count.
tic
clear cfg

%% Define optical parameters

%== Define simulated photon count (suggested less than 1e8 to avoid lengthy simulations)
cfg.nphoton=1e6;

%== Define simulated domain volume (3D array of labels or floating point numbers)
cfg.vol=uint8(ones(60,60,60));
cfg.vol(30:end,:,:)=2;

%== Define optical properties for each tissue label
%         mua(1/mm) mus(1/mm)  g    n
cfg.prop=[   0         0       1    1;       %<- tissue label 0: background    
          0.0005       1       0    1.37     %<- tissue label 1
           0.05        1       0    1.37];   %<- tissue label 2 ...
%--------------------------------------------------------------------------
%%% Hands-on exercise: change the optical properties for tissu 1 and rerun
%--------------------------------------------------------------------------

%== Define source position (in grid-unit, not in mm unit!)
cfg.srcpos=[30 30 0];

%== Define source direction (a unitary vector)
cfg.srcdir=[0 0 1];

%== Define time-gate settings
cfg.tstart=0;   % starting time: 0
cfg.tend=1e-8;  % ending time: 10 ns
cfg.tstep=1e-8; % time step: 10 ns
cfg.seed = rand_seed;


%% Define GPU parameters

cfg.gpuid=1;         % use the first GPU
% cfg.gpuid='11';    % use two GPUs together
cfg.autopilot=1;     % assign the thread and block automatically

%% Define detector

%== Define detector position and radius
cfg.detpos=[30 25 0 2;
            30 20 0 2;
            30 15 0 2;
            30 10 0 2
            30 5 0 2];
% all options: 'dspmxvw'
% for details please type: help mcxlab

%== Define output structure
% d: detector id; s: scattering event count; p: partial path length
% x: exit position; v: exit direction
cfg.savedetflag = 'dspxv';

%% Preview the domain

figure,mcxpreview(cfg);

%% Run simulation

[flux,detpos]=mcxlab(cfg);
t = toc

%% Display deteced photon count

% plot detected photon number vs. detector index
det_pho = histc(detpos.detid,1:size(cfg.detpos,1));

figure,bar(1:5,det_pho,'LineWidth',2)
set(gca,'fontsize',15,'FontWeight','bold')
xlabel('detector index','FontSize',20,'FontWeight','bold')
ylabel('detected photon count','FontSize',20,'FontWeight','bold')