Correcting Impulse Invariance Method
You might be surprised to hear that your code implements neither the conventional impulse invariance method (IIM) nor the corrected IIM. Instead it implements the matched Z-transform, which maps the poles and zeros of the analog prototype according to 

$$z_k=e^{s_kT}\tag{1}$$

where $T=1/f_s$ is the sampling interval.

The two versions of the IIM can be derived as follows. If we take as an example the RIAA de-emphasis transfer function

$$H(s)=\frac{sT_2+1}{(sT_1+1)(sT_3+1)}\tag{2}$$

with $T_1=3.18\,\text{ms}$, $T_2=75\,\mu\text{s}$, and $T_3=318\,\mu\text{s}$, we first have to do a partial fraction expansion of $(2)$ to get

$$H(s)=\frac{A_1}{s-p_1}+\frac{A_2}{s-p_2}\tag{3}$$

with $p_1=-1/T_1$, $p_2=-1/T_3$, and

$$A_1=\frac{1}{T_1}\frac{T_2-T_1}{T_3-T_1},\qquad A_2=\frac{1}{T_3}\frac{T_3-T_2}{T_3-T_1}$$

With the constants used in $(3)$, the transfer function of the transformed system can directly be written as

$$H(z)=\frac{TA_1}{1-e^{p_1T}z^{-1}}+\frac{TA_2}{1-e^{p_2T}z^{-1}}\tag{4}$$

From $(4)$ it can be seen that the poles are the same as for the matched Z-transform, but the zeros are different. Note that for this example there is only one zero for the matched Z-transform as well as for the conventional IIM.

The corrected IIM according to Mecklenbräuker and Jackson subtracts a constant term from the transfer function $(4)$ of the conventional IIM:

$$H(z)=\frac{TA_1}{1-e^{p_1T}z^{-1}}+\frac{TA_2}{1-e^{p_2T}z^{-1}}-\frac{T}{2}(A_1+A_2)\tag{5}$$

Note that this solution has one more zero than the other two methods. The poles of all three methods are identical.

The figure below shows the approximation of the RIAA de-emphasis filter given by $(2)$ according to the three aforementioned transformations (IIM in red, IIM corrected in green, matched Z in magenta). IIM performs worst, matched Z-transform and the corrected IIM are similar but their errors have different signs. This fact could be used to come up with a filter that combines the corrected IIM and the matched Z-transform. Since their poles are identical, we only need to combine their numerator coefficients. The result of the combined filter is shown in black. Its approximation error is less than 1 dB over the whole frequency range.