function [Ct,rv,rp] = schabel(p,t)
% Ct = schabel(p,t)
% implementation of the Schabel model in
% Schabel, A Unified Impulse Resposne Model for DCE-MRI, MRM, 68:(2012)
% this model has been named GCTT (Gamma Capillary Transit time)
% and has the following parameters
% p(1) = tc = capillary transit time [min]
% p(2) = alfa = heterogeneity of the distribution of transit time within a voxel
% p(3) = E = extraction fraction (0<E<1)
% p(4) = kep = rate constant from extra-vascular-space (EES) to plasma [mL/mL/min]
% p(5) = F = blood flow [mL/mL/min]
% 
% 
% t = time  [min] (vector)
%
% Ct(t) = F (Rv(t)+Rp(t)) ** Ca(t)
% where ** is te convolution operator and
% Rv(t) = gamma(tc/tau,t/tau)
% tau = tc / alfa
% Rp(t) = E exp(-kep t) [1 - gamma(tc/tau,(1/tau-kep) t)] / (1 - kep tau)^(tc/tau)
%
% gamma(alfa,z) = 1-gamminc(z,a)


tc   = p(1); % min
alfa = p(2);
E    = p(3); % 0<E<1
kep  = p(4); % min^-1
F    = p(5); % min^-1

tau = tc/alfa;

uno = ones(size(t));
uno(find(t<0)) = 0;
t(find(t<0))=0; % time must be positive


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

Ct(t<0)=0; % before injection Ct is 0

rp = E * exp(-kep * t) .* [ gammainc((1/tau-kep)*t,alfa)] ./ (1 - kep*tau)^(alfa);
rv = uno - gammainc(t/tau,alfa);



	function rc = convolution(x)
if (x>=0)
		Rp = E * exp(-kep * x) .* [ gammainc((1/tau-kep)*x,alfa)] ./ (1 - kep*tau)^(alfa);
		Rv = 1 - gammainc(x/tau,alfa);
else
Rp = 0;
Rv = 0;
endif
		%aif = weinmann(t(k)-x); 
		aif = parker(t(k)-x);
		aif(find(t(k)-x)<0) = 0;
		rc = aif .* (Rp+Rv) * F;
	endfunction

endfunction

