Hello,
thank you for your reply.
Then it seems that "N[HankelH1[16, 3], 25]" is the way to go. Nevertheless, the result from Python
scipy.special.hankel1(16,3) == (-1.0049113641502094e-08-736865544.0285703j)
where the real part significantly differs, leads to a more meaningful physical solution in my specific problem. Nevertheless, as I stated before, this might be only the case due to the fact that the problem was designed by a Python user.
In general, I trust more in Mathematica because there I feel more confident.
However, I wonder if there is really a bug in Python for such a basic function? It is hard to believe for me. Maybe there are different conventions? However, for small n everything seems to be fine which makes me really mistrustful.