%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Demo 1: a basic MMCLAB example
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear cfg;
close all;
%== add path to the softwares
if ispc
    addpath(genpath('../../Software/win64/MCXStudio'));
elseif isunix
    addpath(genpath('../../Software/linux64/MCXStudio'));
elseif ismac
    addpath(genpath('../../Software/osx64/MCXStudio'));
end

%% Create a two-layer tetrahedral mesh model, allocate tissue type label for each element

%== Generate a two-layer triangular surface mesh; calculate centroid coordinates
% for each layer 
[node,face,c0]=latticegrid([0 60],[0 60],[0 15 20]);
%--------------------------------------------------------------------------
%%% Hands-on exercise: change thickness along z-axis (third input parameter) 
%%% for each layer respectively and rerun
%--------------------------------------------------------------------------

%== Define maximum tetrahedral element size for bottom (label 1) and top (label 2) layers
c0(:,4)=[4;2];   
%--------------------------------------------------------------------------
%%% Hands-on exercise: change maxvol for each layer respectively and rerun
%--------------------------------------------------------------------------

%== Generate volumetric mesh from closed surfaces; use centroid as region seeds 
% to create label for each tissue type
[cfg.node,cfg.elem]=surf2mesh(node,face,[],[],1,[],c0);

%== Plot volumetric meshes
figure(1);plotmesh(cfg.node,cfg.elem);title('two-layer mesh');
%--------------------------------------------------------------------------
%%% Hands-on exercise: adjust maximum tetrahedral volume for each layer and  
%%% see how this impact mesh densities and simulation speed
%--------------------------------------------------------------------------

%== Last column of elem is media type
cfg.elemprop=cfg.elem(:,5);

%== First 4 columns of elem are node index for each element
cfg.elem=cfg.elem(:,1:4); 

%== Define optical properties for each tissue type
% simulate the optical properties of skull and gray-matter of a human brain model
% see http://mcx.sourceforge.net/cgi-bin/index.cgi?MMC/Colin27AtlasMesh
%         mua(1/mm) mus(1/mm)  g      n
cfg.prop=[   0          0      1      1;       %<- tissue label 0: background    
            0.02       9.0    0.89   1.37;     %<- tissue label 1: gray-matter
            0.019      7.8    0.89   1.37];    %<- tissue label 2: skull
%--------------------------------------------------------------------------
%%% Hands-on exercise: change the optical properties for tissue 1 or 2 and rerun
%--------------------------------------------------------------------------
        
%% source & detector information

%== define source type as pencil beam
cfg.srctype='pencil';       
cfg.srcpos=[30 30 20];   % source position (in length unit [unitinmm], not in mm unit!)
cfg.srcdir=[0 0 -1];     % pencil beam direction (a unitary vector represents photon initial direction)

%== detector position and radius (in length unit [unitinmm], not in mm unit!)
cfg.detpos=[30 30 0 2]; 

%% Define time-gate settings

cfg.tstart=0;       % starting time: 0
cfg.tend=5e-9;      % ending time: 5 ns
cfg.tstep=1e-10;    % time step: 0.1 ns, results in 50 time windows

%% Define other simulation parameters

%== Define simulated photon count
cfg.nphoton=1e6;        

%== Allocate random number seed (-1 for totally random)
cfg.seed=1648335518;

%== "T" to print timing information, "P" to show progress bar
cfg.debuglevel='TP';

%== Do reflection at internal & external boundaries (refractive index mismatch)
cfg.isreflect=1;        

%% Preview the domain
figure(2);mcxpreview(cfg);title('simulation domain preview');

%% Run mmc simulation
tic;
[layer,detp,~,~]=mmclab(cfg);
toc;
%--------------------------------------------------------------------------
%%% Hands-on exercise: 'Plucker' is the default ray-tracer of mmc, try
%%% other ray-tracers by changing cfg.method, rerun and compare; you can
%%% type 'help mmclab' to learn more about this 
%--------------------------------------------------------------------------

%% output data processing

%== extract time-resolved fluence rate
layer=layer.data;

%== Integrate time-resolved fluence rate to continuous-wave fluence
layercw=sum(layer,2)*cfg.tstep;   

%== CW fluence distribution inside media, along the plane y=30
figure(3);
qmeshcut(cfg.elem,cfg.node,log10(layercw),'y=30','linestyle','none'); view([0 -1 0]);
xlim([0 60]);   zlim([0 20]);   axis equal;
title('log10(CW fluence [J/mm^2]) along plane y=30'); colorbar;

%% Generate TPSF on the detector

%== Calculate remaining weight for detected photons
w = mmcdetweight(detp,cfg.prop);    

%== Calculate time of flight for detected photons
t = mmcdettime(detp,cfg.prop);      

%== number of time gates
nTG = round((cfg.tend-cfg.tstart)/cfg.tstep);

%== Accumulate remaining weight of detected photons to their corresponding time-gate 
tpsf = zeros(nTG,1);
ntg = ceil((t-cfg.tstart)/cfg.tstep);
for i=1:length(w)
    tpsf(ntg(i)) = tpsf(ntg(i))+w(i);
end

figure(4);
plot([1:nTG]*cfg.tstep/1e-9,tpsf/max(tpsf));
xlabel('time (ns)');    ylabel('normalized TPSF');
title('Normalized TPSF for detecter #1');
