  Computation of electromagnetic fields due to point sources (Hertzian dipoles)
in cylindrically stratified media is a classical problem for which analytical
expressions of the associated tensor Green's function have been long known.
However, under finite-precision arithmetic, direct numerical computations based
on the application of such analytical (canonical) expressions invariably lead
to underflow and overflow problems related to the poor scaling of the
eigenfunctions (cylindrical Bessel and Hankel functions) for extreme arguments
and/or high-order, as well as convergence problems related to the numerical
integration over the spectral wavenumber and to the truncation of the infinite
series over the azimuth mode number. These problems are exacerbated when a
disparate range of values is to be considered for the layers' thicknesses and
material properties (resistivities, permittivities, and permeabilities), the
transverse and longitudinal distances between source and observation points, as
well as the source frequency. To overcome these challenges in a systematic
fashion, we introduce herein different sets of range-conditioned, modified
cylindrical functions (in lieu of standard cylindrical eigenfunctions), each
associated with non-overlapped subdomains of (numerical) evaluation to allow
for stable computations under any range of physical parameters. In addition
adaptively-chosen integration contours are employed in the complex spectral
wavenumber plane to ensure convergent numerical integration in all cases. We
illustrate the application of the algorithm to problems of geophysical interest
involving layer resistivities ranging from 1000 $\Omega \cdot$m to 10$^{-8}
\Omega \cdot$m, frequencies of operation ranging from 10 MHz down to the low
magnetotelluric range of 0.01 Hz, and for various combinations of layer
thicknesses.
