Abstract: Recent methods to map artificial night-sky brightness and stellar visibility across large territories or their distribution over the entire sky at any site are based on computation of the propagation of light pollution with Garstang models, a simplified solution of the radiative transfer problem in the atmosphere that allows fast computation by reducing it to a ray-tracing approach. They are accurate for a clear atmosphere, when a two-scattering approximation is acceptable, which is the most common situation. We present here up-to-date extended Garstang models (EGM), which provide a more general numerical solution for the radiative transfer problem applied to the propagation of light pollution in the atmosphere. We also present the LPTRAN software package, an application of EGM to high-resolution Defense Meteorological Satellite Program (DMSP) Operational Linescan System (OLS) satellite measurements of artificial light emission and to GTOPO30 (Global 30 Arcsecond) digital elevation data, which provides an up-to-date method to predict the artificial brightness distribution of the night sky at any site in the world at any visible wavelength for a broad range of atmospheric situations and the artificial radiation density in the atmosphere across the territory. EGM account for (i) multiple scattering, (ii) wavelengths from 250 nm to infrared, (iii) the Earth's curvature and its screening effects, (iv) site and source elevation, (v) many kinds of atmosphere with the possibility of custom set-up (e.g. including thermal inversion layers), (vi) a mix of different boundary-layer aerosols and tropospheric aerosols, with the possibility of custom set-up, (vii) up to five aerosol layers in the upper atmosphere, including fresh and aged volcanic dust and meteoric dust, (viii) variations of the scattering phase function with elevation, (ix) continuum and line gas absorption from many species, ozone included, (x) up to five cloud layers, (xi) wavelength-dependent bidirectional reflectance of the ground surface from National Aeronautics and Space Administration (NASA) Moderate-Resolution Imaging Spectroradiometer (MODIS) satellite data, main models or custom data (snow included) and (xii) geographically variable upward light-emission function given as a three-parameter function or a Legendre polynomial series. Atmospheric scattering properties or light-pollution propagation functions from other sources can also be applied. A more general solution allows us to account also for (xiii) mountain screening, (xiv) geographical gradients of atmospheric conditions, including localized clouds and (xv) geographic distribution of ground surfaces, but suffers from too heavy computational requirements. Comparisons between predictions of classic Garstang models and EGM show close agreement for a US62 standard clear atmosphere and typical upward emission function.