%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MCX19 Session 9  - Advanced Features 6
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Photon Replay
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% first run a baseline simulation

clear cfg cfgs
cfg.nphoton=1e5;
cfg.vol=uint8(ones(60,60,60));
cfg.srcpos=[30 30 0];
cfg.srcdir=[0 0 1];
cfg.gpuid=1;
% cfg.gpuid='11'; % use two GPUs together
cfg.autopilot=1;
cfg.issrcfrom0=1;
cfg.prop=[0 0 1 1;0.005 1 0 1.37];
cfg.tstart=0;
cfg.tend=5e-9;
cfg.tstep=5e-10;
cfg.issaveexit=1;
% calculate the flux distribution with the given config
cfg.detpos=[15 30 0 0.5];

%% save the detected photon seeds

[flux, detp, vol, seeds]=mcxlab(cfg);

%% pass the seeds tp cfg.seed
% pass the detps partial path to cfg.detphotons, then set output to jacobian

newcfg=cfg;
newcfg.seed=seeds.data;
newcfg.outputtype='jacobian';
newcfg.detphotons=detp.data;
newcfg.maxjumpdebug=1000000;
[flux2, detp2, vol2, seeds2, traj2]=mcxlab(newcfg);

%% obtain cw jacobian

jac=sum(flux2.data,4);
imagesc(log10(abs(squeeze(jac(:,30,:)))))

%% verify with trajectories

figure;
mcxpreview(cfg);
hold on;
newtraj=mcxplotphotons(traj2,'color','b');title('photon trajectories');
set(gca,'zlim',[0 10])

