%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   MCX/MMC advanced features - Demo 1,2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Demo 1 - beam focal length

USE_MCXCL=0;

clear mcxfg mmcfg
close all

% define the basic parameters for an MCX simulation

mcxfg.nphoton=1e6;
mcxfg.issrcfrom0=1;

mcxfg.vol=ones(40,40,40);
mcxfg.vol(end,end,end)=0;  %%< make this corner more complex
mcxfg.tstart=0;
mcxfg.tend=2e-9;
mcxfg.tstep=2e-9;

% define a disk shaped source

mcxfg.srctype='disk';      %%<- any area source can specify focal length
mcxfg.srcpos=[30 30 0];
mcxfg.srcparam1=[20,0 0 0];
mcxfg.srcdir=[0,0,1];      %% no focal length specified: collimated beam

mcxfg.prop=[0 0 1 1; 
    0.005, 0.01, 0.01, 1.37];

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

mcxfg.savedetflag='dx';

%% Run mcxlab simulation

[flux, detps]=mcxlab(mcxfg);

%% Plot the photon exit positions, see the shapes
figure;
plotmesh(detps.p,'.')
mcxplotvol(log10(flux.data))

%% Define focal length as the 4th number of cfg.srcdir - positive means converting, negative means diverging

mcxfg.srcdir=[0,0,1,100];  %% focus at 100 mm above cfg.srcpos
mcxfg.isreflect=0;   %%<- Hands-on suggestions, disable this and see changes

%% rerun simulation

[flux, detps]=mcxlab(mcxfg);

%% plot the results
mcxplotvol(log10(flux.data))

%%%========================================================================
%%   MMC detector 
%%%========================================================================

%% create a basic MMC session
mmcfg=mcxfg;
mmcfg=rmfield(mmcfg,{'vol','savedetflag','issrcfrom0'});
density=3;
maxvol=20;
[mmcfg.node,fc, mmcfg.elem]=meshasphere([30,30,30],30,density,maxvol);
mmcfg.elemprop=ones(size(mmcfg.elem,1),1);
mmcfg.srcdir(4)=[];


mmcfg.nphoton=1e7;

mmcfg.srctype='pencil';
mmcfg.srcpos=[30 30 0.1];
mmcfg.srcparam1=0;
mmcfg.e0=tsearchn(mmcfg.node,mmcfg.elem,mmcfg.srcpos);
mmcfg.prop=[0 0 1 1; 
    0.005, 1, 0.01, 1.37];

mmcfg.detpos=[30,30,60 10];
mmcfg.method='grid';            %% enable DMMC output as a volume

mmcfg.issaveexit=1;             %% use this flag to save exit position
mmcfg.debuglevel='P';

%%
[mmcflux, mmcdetps]=mmclab(mmcfg);  %% run MMCLAB

%%
mcxplotvol(log10(mmcflux.data))
figure;
plotmesh(mmcdetps.p,'.')      %% plot mmc exit positions
