%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Demo 2: source configuration and time-gate settings
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

%== 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];


%----------------- pencil beam -----------------
cfg.srctype = 'pencil';

%----------------- planar -----------------
% cfg.srctype = 'planar';
% % x/y/z of planar corner
% cfg.srcparam1 = [20 0 0 0]; % in grid-unit, not in mm unit!
% cfg.srcparam2 = [0 10 0 0];

%----------------- gaussian -----------------
% cfg.srctype = 'gaussian';
% cfg.srcparam1(1) = 5;   % gaussian waist radius (in grid-unit, not in mm unit!)
% cfg.srcpos=[30 30 0];

%== Define time-gate settings
cfg.tstart=0;   % starting time: 0
cfg.tend=10e-9;       % ending time: 10 ns
cfg.tstep=0.5e-9;      % time step: 0.5 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

figure,mcxpreview(cfg);

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

%% Display fluence distributions

figure,imagesc(flipud(squeeze(log10(sum(flux.data(:,30,:,1:end),4)))'),caxis),colormap jet,axis off;
axis tight
c = colorbar;
c.FontSize = 15;
xlabel('z axis'),ylabel('x axis')
xlabel('z axis'),ylabel('x axis')
title('CW solution')

%% Display fluence at each time-gate

caxis = [2 8];
figure,imagesc(flipud(squeeze(log10(flux.data(:,30,:,1)))'),caxis),colormap jet,axis off;
axis tight
c = colorbar;
c.FontSize = 15;
xlabel('z axis'),ylabel('x axis')
title([cfg.srctype ' at 0.5 ns'])

figure,imagesc(flipud(squeeze(log10(flux.data(:,30,:,2)))'),caxis),colormap jet,axis off;
axis tight
c = colorbar;
c.FontSize = 15;
xlabel('z axis'),ylabel('x axis')
title([cfg.srctype ' at 1 ns'])

figure,imagesc(flipud(squeeze(log10(flux.data(:,30,:,4)))'),caxis),colormap jet,axis off;
axis tight
c = colorbar;
c.FontSize = 15;
xlabel('z axis'),ylabel('x axis')
xlabel('z axis'),ylabel('x axis')
title([cfg.srctype ' at 2 ns'])

figure,plot(0.5:0.5:10,squeeze(flux.data(30,30,30,:)),'LineWidth',2)
set(gca,'fontsize',15,'FontWeight','bold')
xlabel('Time/$ns$','interpreter','latex','FontSize',20,'FontWeight','bold')
ylabel('Fluence rate/$J\cdot s^{-1}$','interpreter','latex','FontSize',20,'FontWeight','bold')
