Non-continuous kink in photon energy flux vs. temperature

[jrsmith3]

Disclaimer: Its been awhile since I wrote the original code examples for this issue.

Problem

There's a non-continuous kink in the photon energy flux calculation vs. temperature (cf. this ipython notebook). The problem is that the photon energy flux is calculated via a finite sum of polylogarithm functions and the first order polylogarithm function returns values of zero for sufficiently small arguments. At some value of argument, the return value suddenly jumps to a nonzero value and yields the kink seen in the ipython notebook.

In ibei, the polylogarithm is implemented by the mpmath package via sympy. The first order polylogarithm can be expressed as a natural log.

Li_{1}(z) = -log(1 - z)

mpmath's implementation of the first order polylogarithm tests for the case of order unity and calls the natural log function in this case. The problem is that the argument to the natural log function is unity minus the polylogarithm argument. If the original polylogarithm argument is very small, 1 - z will be rounded to unity yielding zero for the final result (log(1) = 0).

>>> import mpmath
>>> z = 5e-20
>>> mpmath.polylog(1, z)
0

Solution

The solution is simply to increase the precision of the mpmath calculation to the order of the polylogarithm argument.

>>> import mpmath
>>> z = 5e-20
>>> mpmath.mp.dps = 22
>>> mpmath.polylog(1, z)
mpf('5.000000487562466251661797e-20')

This modification can be done within the uibei method by evaluating the size of the real_arg value.