function tpsf = mcxdettpsf(detp, detnum, prop, time)
%
% tpsf=mcxdettpsf(detp,detnum,prop,time)
%
% Calculate the temporal point spread function curve of a specified detector
% given the partial path data, optical properties, and distribution of time bins
%
% author: Qianqian Fang (q.fang <at> neu.edu)
%     Ruoyang Yao (yaor <at> rpi.edu)
%
% input:
%     detp: the 2nd output from mcxlab. detp must be a struct with detid and ppath subfields
%     detnum: specified detector number
%     prop: optical property list, as defined in the cfg.prop field of mcxlab's input
%     time: distribution of time bins, a 1*3 vector [tstart tend tstep]
%           could be defined different from in the cfg of mcxlab's input
% output:
%     tpsf: caculated temporal point spread function curve of the specified detector
%
% this file is copied from Mesh-based Monte Carlo (MMC)
%
% License: GPLv3, see http://mcx.space/ for details
%

% select the photon data of the specified detector
detp.ppath = detp.ppath(detp.detid == detnum, :);
detp.detid = detp.detid(detp.detid == detnum);

% calculate the detected photon weight and arrival time
replayweight = mcxdetweight(detp, prop);
replaytime = mcxdettime(detp, prop);

% define temporal point spread function vector
nTG = round((time(2) - time(1)) / time(3));     % maximum time gate number
tpsf = zeros(nTG, 1);

% calculate the time bin, make sure not to exceed the boundary
ntg = ceil((replaytime - time(1)) / time(3));
ntg(ntg < 1) = 1;
ntg(ntg > nTG) = nTG;

% add each photon weight to corresponding time bin
for i = 1:length(replayweight)
    tpsf(ntg(i)) = tpsf(ntg(i)) + replayweight(i);
end

end