function G = diffusionProb(r,t,z,mu,a,sigma)
%calculate Pr[response = r, RT <= t] for a diffusion model with parameters
% z = starting point
% mu = drift rate
% a = upper threshold (lower threshold 0)
% sigma^2 = diffusion rate

%get marginal probability of response r
if r==1
    P = (1-exp(-2*mu*z/sigma^2))/(1-exp(-2*mu*a/sigma^2));
else
    P = (exp(-2*mu*z/sigma^2)-exp(-2*mu*a/sigma^2))/(1-exp(-2*mu*a/sigma^2));
end

%get desired response-RT probability
A = pi*sigma^2/a^2;
B = exp((2*mu*[a-z;-z]-t*mu^2)/(2*sigma^2)); %one value for response 1 (a-z), one for response 2 (-z)
C = -pi^2*sigma^2*t/(2*a^2);
D = pi*[a-z;z]/a; %one value for response 1 (a-z), one for response 2 (z)
threshold = 1e-40; 
maxk = sqrt(1/C*log(threshold)); %exp(maxkn^2*C) = threshold
temp = 0;
for k=1:maxk
    temp = temp + 2*k*a^2*sigma^2/(k^2*pi^2*sigma^4+a^2*mu^2)*sin(k*D(r))*exp(k^2*C);
end
G = P - A*B(r)*temp;
