I have used SageMath to symbolically integrate an expression. The result contains a gamma function with two input parameters:
gamma(-1, 2*((x - xp)^2 + (y - yp)^2)/s^2)
Apparently this is called an incomplete gamma function. Now I want to use this expression in a Python code. I have tracked down the incomplete gamma function to scipy.special.gammainc. Unfortunately, this function does not allow negative input parameters and I have to use -1 as the first input parameter. How can I work around this?
The lower incomplete gamma function can be defined in terms of an improper integral according to Wikipedia. This integral can be related to the generalised form of the exponential integral. Both pages give this relation between the two:
E_n(x) = x^(n-1)*gamma(1-n, x)
So for the case in OP we would have n=2, which corresponds to -1 as a first input parameter for the gamma function. I have numerically verified in SageMath that the above relation holds up. The corresponding functions in SageMath are:
1. gamma(n, x) == gamma_inc(n, x)
2. E_n(x) == exp_integral_e(n, x)
Which according to the Wikipedia relation gives us (aside from round-off errors):
exp_integral_e(n, x) == x^(n-1)*gamma_inc(1-n, x)
The corresponding Python functions are:
1. gamma(n, x) == scipy.special.gammainc(n, x)
2. E_n(x) == scipy.special.expn(n, x)
Which according to the Wikipedia relation gives us (aside from round-off errors):
expn(n, x) == x**(n-1)*gammainc(1-n, x)
There is one small caveat. The gamma_inc function from SageMath accepts a negative first input parameter, whereas the gammainc function from scipy.special does not. However, the expn function from scipy.special does not have this limitation as it can be evaluated for n>=2 corresponding to a negative first input parameter for gamma_inc.
So the answer to the OP is to use the Wikipedia relation to replace the lower incomplete gamma function with the generalised exponential integral and to use scipy.special.expn for evaluation in Python.