%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% illuminate_PreXion.m                                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Illuminate 3D dataset by PreXion system, simulated by MCX.             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function fluence = illuminate_PreXion(object_data,object_properties,volume_center,object_unitinmm)
clear cfg;

% define global MCX parameters
cfg.gpuid=1;
cfg.autopilot=1;
cfg.tstart=0;
cfg.tend=4e-10;
cfg.tstep=4e-10;
cfg.seed=0;
cfg.isreflect=1;
cfg.isnormalized = 1;
cfg.isgpuinfo=0;

totalnphoton = 1e5; % original: 62*3*30720

% prepare MCX data
cfg.vol = object_data;
dim = size(object_data);
% cfg.prop = zeros(size(object_properties,1),4);
cfg.prop = object_properties;
cfg.unitinmm = object_unitinmm;

% partition number of photons
cfg.nphoton=ceil(totalnphoton/2);
% change dimensions in order to agree with MCX standard
cfg.vol = uint8(permute(cfg.vol,[2,1,3]));
flux_sum = zeros(dim(2),dim(1),dim(3));
failed_simulations = 0;

% define specific PreXion source properties, in mm
cfg.srctype='planar';
zheight=5; % Distance from surface of sample
sep_LED=7; % Separation between LED sources
len_src=50; % Dimensions of LED arrays
wid_src=7;  % % % 

% Two planar sources focused on surface of phantom
for i=1:2
    sideStep=((-1)^i)*sep_LED/cfg.unitinmm; % Two LED arrays side by side around transducer, distance from transducer at center
    cfg.srcpos= [volume_center(2)-(wid_src/2)/cfg.unitinmm+sideStep; volume_center(1)-(len_src/2)/cfg.unitinmm; volume_center(3)+(15+zheight)/cfg.unitinmm]; 
    cfg.srcparam1=[0 len_src/cfg.unitinmm 0 0];                                      % 15+50 because 15 is surface of phantom, .5cm away from surface
    cfg.srcparam2=[wid_src/cfg.unitinmm  0 0 0];
    dir_vec=[-sideStep 0 -zheight/cfg.unitinmm];
    cfg.srcdir=dir_vec/norm(dir_vec);
    % Calculate flux, sum
    flux=mcxlab(cfg);
    if (sum(isnan(flux.data(:))+isinf(flux.data(:)))==0)
        cwf=sum(flux.data,4);
        flux_sum = flux_sum + cwf;
    else
        failed_simulations = failed_simulations + 1;
    end
end
disp(failed_simulations)
fluence = permute(abs(double(flux_sum)),[2,1,3]); %change back coordinates to match object_data

end