function [Ct] = tofts(p,t)
% Tofts model
% Fusco et al. Journal of Biomedical Graphics and computing, 2(2);2012

ktrans = p(1); % min^-1
kep = p(2); % min^-1
vp = p(3); % volume fraction 0<vp<1

Ct = zeros(size(t));

for k = 1:length(t)
	[Ct(k)] = quad("convolution",0,t(k));
end

Ct = Ct + vp * weinmann(t);
Ct(t<0)=0; % before injection Ct is 0


	function rc = convolution(tc)
		rc = weinmann(tc) .* ktrans .* exp(-kep*(t(k)-tc));
	endfunction

endfunction




