% Demo 1: you need to fill in some important parameters to get you familar
% with MCXLAB

clear all
close all

rand_seed = randi([1 2^31-1], 1, 1);

% Create matrix that saves the data for a specific photon count.
tic
clear cfg

% TODO: the photon count (no more than 1e9)
cfg.nphoton=1e6;

% spherical volume design
x = 60;
y = 60;
z = 60;
[yy xx zz] = meshgrid(1:x,1:y,1:z);
S1 = double(sqrt((xx-30).^2+(yy-30).^2+(zz-20).^2)<=15);
S2 = double(sqrt((xx-30).^2+(yy-30).^2+(zz-20).^2)<=10);
S = S1+S2;
S = S+1;
volume = uint8(S);
cfg.vol=volume;

% optical properties
cfg.prop=[0 0 1 1;
          0.0005 1 0 1.37;
          0.0005 1 0 1.37*2;
          0.3 1 0 1.37];

% TODO: define the source position
cfg.srcpos=[30 30 1];

% TODO: define the source direction
cfg.srcdir=[0 0 1];

% different source type
    %%% Gaussian %%%
    cfg.srctype = 'gaussian';
    cfg.srcparam1 = 10;
    
    %%% Pencil %%%
%     cfg.srctype = 'pencil';
    
    %%% Disk %%%
%     cfg.srctype = 'disk';
%     cfg.srcparam1 = 10;

    %%% Isotropic %%%
%     cfg.srctype = 'isotropic';

    %%% Cone %%%
%     cfg.srctype = 'cone';
%     cfg.srcparam1 = pi/6;

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

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;

% TODO: RUN mcxlab
[flux,detpos]=mcxlab(cfg);
t = toc

caxis = [2 8];
figure,imagesc(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 spherical inclusion'],'FontSize', 15)