% Function for evaluating equation (6) in Reiss (2016) -- Imhof's integral
% u is the integration variable
% arg is the upper limit of the original integral (Pr(q < arg))
% lambd = quadratic form weights
% multip = n_r (multiplicities)
% delt = noncentrality deltas (nroot by 1 here)
function integ=imhof_integrand(u,arg,lambd,multip,delt)
nroots=size(lambd(:),1);
integ=0;
theta=-arg.*u/2;
rho_u=1;
for j=1:nroots
    theta=theta+(multip(j)/2)*atan(lambd(j).*u)+lambd(j)*delt(j)*u./(2*(1+lambd(j)^2*u.^2));
    rho_u=(rho_u.*(1+lambd(j)^2*u.^2).^(multip(j)/4)).*exp(lambd(j)^2*delt(j)*u.^2./...
        (2.*(1+lambd(j)^2.*u.^2)));
end;
integ=integ+sin(theta)./(u.*rho_u);