function U = upsilon(q, p)

U   = 1 + (p.sigma - 1)*exp(1/p.epsi)*(p.epsi^((p.sigma/p.epsi) - 1))*(uig(p.sigma/p.epsi, 1/p.epsi) - uig(p.sigma/p.epsi, (q.^(p.epsi/p.sigma))/p.epsi));


function Gamma_sx = uig(s,x)

% Gamma(s)   = \int_0^{\infty} t^(s-1) e^(-t) dt ; gamma function
% gamma(s,x) = \int_0^{   x  } t^(s-1) e^(-t) dt ; lower incomplete gamma
% Gamma(s,x) = \int_x^{\infty} t^(s-1) e^(-t) dt ; upper incomplete gamma

% Gamma(s)   = Gamma(s,0)
% Gamma(s,x) = Gamma(s) - gamma(s,x)


Gamma_s      = gamma(s);

gamma_sx     = (Gamma_s)*gammainc(x,s);          % Matlab uses the "regularized" lower incomplete gamma

Gamma_sx     = Gamma_s - gamma_sx;
