%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   MCX/MMC advanced features - Demo 3 diffuse reflectance
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

USE_MCXCL=0;

clear mcxfg
close all

% define the basic parameters for an MCX simulation

mcxfg.nphoton=1e7;
mcxfg.issrcfrom0=1;     %%<- this defines the lower-most corner as [0,0,0] instead of [1,1,1]

mcxfg.vol=ones(40,40,40);
mcxfg.tstart=0;
mcxfg.tend=2e-9;
mcxfg.tstep=2e-9;

% define a disk shaped source

mcxfg.srctype='disk';
mcxfg.srcpos=[30 30 0];
mcxfg.srcparam1=[20,0 0 0];
mcxfg.srcdir=[0,0,1];

mcxfg.prop=[0 0 1 1; 
    0.005, 0.1, 0.01, 1.37];    %%<- define a low scattering medium to see beam profile

mcxfg.detpos=[
    20,20,40 2; 
    40,40,40 2; 
    20,40,40 2
    ];

mcxfg.savedetflag='dx';  %%<- output detected photon exit position

mcxfg.vol(:,:,1)=0;      %% key: to output dref, you must pad the medium by a layer of 0s
mcxfg.issaveref=1;       %% key: you must set this flag

%% Run simulation to get flux.dref and plot

[flux, detps]=mcxlab(mcxfg);

mcxplotvol(flux.dref)

%% Demo 4: continuous media

[xx,yy,zz]=ndgrid(1:40,1:40,1:40); 
mua=xx/40;
mus=yy/40*10;

mcxfg.vol=single(reshape([mua(:) mus(:)]',[2 size(xx)]));  %%<- a 4-D array vol defines continuous medium

size(mcxfg.vol)

%%

[flux2, detps2]=mcxlab(mcxfg);

mcxplotvol(flux2.data)


%% Demo 5: specify bounary condition

mcxfg.vol=uint8(ones(40,40,40));

mcxfg.bc='aa_mmm';   %%<- 'a' for absorbing BC, 'm' for mirror BC, '_' to fallback to cfg.isreflect

%%

[flux3]=mcxlab(mcxfg);
mcxplotvol(flux3.data)

