%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  mcxyz skinvessel benchmark
%
%  must change mcxyz maketissue.m boundaryflag variable from 2 to 1 to get
%  comparable absorption fraction (40%), otherwise, mcxyz obtains slightly
%  higher absorption (~42%) with boundaryflag=2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear cfg flux

%% create the skin-vessel benchmark mesh
[no,fc]=latticegrid([0 200],[0 200],[20 32 200]); % create a 3-layer tissue
no(end,:)=no(end,:)+1e-5;

fc2=cell2mat(fc);
fc=[fc2(:,[1 2 3]); fc2(:,[1 3 4])];

[ncy,fcy]=meshacylinder([-1,99.5,99.5],[201,99.5,99.5],20,5); % add the vessel
[newnode,newelem]=surfboolean(no,fc,'first',ncy,fcy);  % merge the two domains

c0=[10,100,150,26]';
seeds=[ones(4,2)*100, c0];  % define the regions by index

%ISO2MESH_TETGENOPT='-Y -A'
[cfg.node,cfg.elem]=s2m(newnode,newelem(:,1:3),1,30,'tetgen',seeds,[]); % creating the merged mesh domain

cfg.unitinmm=0.005;
cfg.method='elem';

figure;
subplot(121);
plotmesh(cfg.node,cfg.elem);

cfg.elemprop=cfg.elem(:,5);
cfg.elem=cfg.elem(:,1:4);

%% define other properties

cfg.prop=[0.0000         0.0    1.0000    1
    3.5640e-05    1.0000    1.0000    1.3700
   23.0543    9.3985    0.9000    1.3700
    0.0458   35.6541    0.9000    1.3700
    1.6572   37.5940    0.9000    1.3700];

cfg.srcpos=[100 100 -1];
cfg.srcdir=[0 0 1];

cfg.tstart=0;
cfg.tend=5e-8;
cfg.tstep=5e-8;
% cfg.outputtype='energy'; %energy deposition in mmc varys with elem volume
cfg.outputtype='flux';
cfg.minenergy=0.01;

cfg.srctype='disk';
cfg.srcparam1=[0.3 0 0 0]/cfg.unitinmm; % in grid unit

%% define wide-field disk source by extending the mesh to the widefield src
srcdef=struct('srctype',cfg.srctype,'srcpos',cfg.srcpos,'srcdir',cfg.srcdir,...
    'srcparam1',cfg.srcparam1,'srcparam2',[]);

[cfg.node,cfg.elem] = mmcaddsrc(cfg.node,[cfg.elem cfg.elemprop],...
    mmcsrcdomain(srcdef,[min(cfg.node);max(cfg.node)]));


cfg.elemprop=cfg.elem(:,5);
cfg.elem=cfg.elem(:,1:4);

%% other simulation information

cfg.nphoton=1e7;
cfg.seed=1648335518;

cfg.debuglevel='TP';
cfg.isreflect=0;

%% mmc simulation

flux=mmclab(cfg);
flux=flux.data;
fluxcw=sum(flux,2)*cfg.tstep*100; % multiplying 100 converts 1/mm^2 to 1/cm^2 to match mcxyz output

%% plot simulated photon profiles

subplot(122);
hold on;
qmeshcut(cfg.elem(cfg.elemprop>0,1:4),cfg.node*cfg.unitinmm,log10(fluxcw),'x=0.5','linestyle','none');
view([1 0 0]);
set(gca,'zlim',[0 1],'ylim',[0 1],'zdir','reverse')

box on;
axis equal
title('MMC fluence rate (W/cm^2) per W simulated')
colorbar;
colormap(jet)
