%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Demo 1: a basic MCXLAB example
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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));

%== 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 ...
%--------------------------------------------------------------------------
%%% 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


%% Preview the domain

mcxpreview(cfg);

%% Run simulation

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

%% Display fluence distributions

caxis = [2 8];
figure,imagesc(flipud(squeeze(log10(flux.data(:,30,:)))'),caxis),colormap jet,axis off;
axis tight
c = colorbar;
c.FontSize = 15;
xlabel('z axis'),ylabel('x axis')
title(['1e' num2str(log10(cfg.nphoton)) ' photons for a homogeneous cube'],'FontSize', 15)

%% Demo 1.2: adding a new tissue label

%== Define a new tissue label by modifying the cfg.vol 3D array
cfg.vol(30:end,:,:)=2;

%== Append a new row in cfg.prop to define the corresponding optical
%properties

cfg.prop(3,:)=[0.05 1 0 1.37]; 

%== Let's show a progress bar this time
cfg.debuglevel='P';

%% Preview the modified domain

mcxpreview(cfg);

%% Rerun the simulation

[flux2,detpos2]=mcxlab(cfg);

%% Plot the results for the 2-layer media
caxis = [0 8];
figure,imagesc(flipud(squeeze(log10(flux2.data(:,30,:)))'),caxis),colormap jet,axis off;
axis tight
c = colorbar;
c.FontSize = 15;
xlabel('z axis'),ylabel('x axis')
hold on;
plot([30 30],[0 60],'color','w','linestyle','--','linewidth',2)
title(['1e' num2str(log10(cfg.nphoton)) ' photons for a two-layered domain'],'FontSize', 15)