Angular Emission Function of a City and Skyglow Modeling: A Critical Perspective

Light emissions from ground-based sources have been intensively studied over the last decades, aiming to better understand the mechanisms of how the light pollution spreads in a nocturnal environment while causing an overillumination of living spaces, an excessive exposure of light over time, or just an illumination of otherwise dark places. Artificial lights such as municipal lighting installations, lights from buildings, or lights from cars can have different spectra, intensity, and angular emission function, but also spatial distribution, and can even be characterized by different peak hours (Elvidge et al. 2010). Most of street lighting luminaires are designed to redirect the light beams downward to ensure good viewing conditions in cities or towns (Novák et al. 2014), while the remaining emissions are considerably influenced by buildings, trees, and other obstacles that block off part of the light directed to low elevation angles (Duriscoe et al. 2014; Bará et al. 2015b). The illumination scene inside a city or town may strongly depend on the spectral composition of light sources surrounding a measuring site. However, the diffuse light of a night sky is mostly due to the superposition of many optical signals from all city lights. In principle, each elementary area on the ground inside of a city can contribute to the skyglow directly (thanks to direct emissions upwards) or indirectly (through diffuse reflection).

The light intensity causing an illumination of an atmospheric elementary volume is a non-trivial sum of intensities of all beams emitted from different city-areas under different zenith angles (Kocifaj 2007, 2014b, and Bará et al. 2015a). Hence, the upward light is due to superimposed emission distributions. Reliable predictions of sky radiance patterns can only be made when the bulk radiant intensity as a function of zenith angle is generally known for a city, town, or any other form of settlement that emits the light.

The angular emission function is one of the key properties that predetermines the variability of the diffuse light of a night sky (Aubé 2015), and there is a substantial demand for obtaining this function for many cities worldwide. While some progress has been made in the past several years (Kyba et al. 2013), the city emission function (CEF) still remains poorly quantified. The first works on the CEF date back to Garstang (1986), who developed an empirical model based on observations in Denver, Colorado. Garstang's formula appeared to be suitable for modeling purposes, but most likely the Garstang's emission function (GEF) does not satisfy light emissions from many other cities with different lighting installations, urbanistic design, orography, etc. Undoubtedly, the GEF also overestimates the emissions to low-elevation angles; Luginbuhl et al. (2009b) have treated large-scale blocking by terrain features, and introduced a "blocking" extinction term to overcome this shortcoming. The extinction term is important for distant light sources, but may not necessarily work properly if used as a scaling factor for both uplight and reflected components, especially at low zenith angles.

The CEF is a common input to all radiative transfer (RT) models that simulate artificial light in a nighttime environment, in particular the diffuse light of a night sky. Irrespective of what method is used for solving the standard RT problems, e.g., an adding or doubling method, a method of successive orders, Ambartsumian's method, a Monte Carlo or various other discretization methods (van de Hulst 1980), a numerical computation is the only effective way to give realistic predictions for spectral and broadband radiances and irradiances under arbitrary meteorological conditions. Solving the radiative transfer equation (RTE) can significantly contribute to our understanding of the extent to which the artificial light can affect living organisms and the environment (Longcore & Rich 2004; Gaston et al. 2014). Because the RTE is a fundamental evaluation tool, it makes a quantification of illumination levels easily possible and allows for the numerical comparison of different lighting technologies.

However, any misuse of the CEF may cause the computational results appear inaccurate or even partially or fully incorrect. Unfortunately, the CEF is used non-uniformly over the scientific literature, and such ambiguity is always a potential source of errors. Among the typical misconceptions is that what is placed into boundary conditions for the RTE is radiant intensity rather than radiance. It is therefore the goal of this paper to provide a solid foundation of the CEF and point out systematic errors that originate from incorrect use of CEF in skyglow modeling.

The bulk emission function still remains largely unknown for large agglomerations; however, it is often characterized by a radiant intensity function independent of the size of a light-emitting area,

$$\begin{eqnarray}&&{I}_{{up},\lambda }(z)=\displaystyle \frac{d{{\rm{\Phi }}}_{\lambda }}{d{\rm{\Omega }}}.\end{eqnarray} \tag{ 1 }$$

Here, dΦ_λ is an elementary spectral radiant flux passing through a small solid angle, dΩ, so the spectral radiant intensity I_up,λ(z) is measured in W·nm⁻¹·sr⁻¹ (or in W·sr⁻¹, taking into account the broadband equivalent of I_up,λ). The above approach is well suited for a small light source compared with any characteristic length in the system, e.g., when compared with the distance between a city and an observer. In fact, cities or towns are large emitting surfaces, so the optical signals recorded by a detector would change with a geometric configuration. Let us consider a flat surface, dA, of uniform light in a simple thought experiment (see Figure 1).

Zoom In Zoom Out Reset image size

The electromagnetic energy of all photons emitted from the ground area dA and reaching the detecting area dA₀ oriented normally to the beams of light is i_λcos(z)dΩdA, where i_λ is spectral radiance (W·m²·nm⁻¹·sr⁻¹) or, alternatively, its broadband equivalent (W·m²·sr⁻¹); z is the angle contained by the surface normal and the direction of the emission (actually the zenith angle); and dΩ equals to dA₀r⁻² with r being the distance between dA and dA₀. A measuring equipment with a sufficiently wide field of view, $d\tilde{{\rm{\Omega }}}\gt {{dAr}}^{-2}$ would record all photons emitted from dA toward dA₀. Therefore, the total energy crossing the detector r is i_λcos(z)dA₀r⁻²dA = i_λcos(z)dA₀ dΩ₀, where dΩ₀ = dAr⁻². Using a comparison of what is emitted with what is detected, one can conclude that radiance is conserved along lines through non-attenuating media; however, the number of photons and also the radiant energy are both proportional to the cosine of the angle between the direction of the emitted light and the surface normal (see, e.g., Thomas & Stamnes 2006). Note that surface normal is the vector that is normal to the light-emitting plane. However, radiance and radiant intensity function are sometimes imprecisely distinguished. Such vaguely specified inputs are potential sources of errors in the interpretation of the CEFs and their consequent skyglow effects.

For instance, the three-dimensional (3D) RTE should be solved subject to boundary conditions to obtain the diffuse radiation field that arose due to the interaction of direct beams of light with the atmosphere. A basic practice of the radiative transfer modeling is to minimize the complexity of computations, thus the rigorous vector RTE is traditionally approximated by its scalar one-dimensional (1D) counterpart assuming the incident light is monochromatic and unpolarized. With the aforementioned assumptions, the RTE would read:

$$\begin{eqnarray}&&\cos \xi \displaystyle \frac{d\,{i}_{\lambda }({\tau }_{\lambda },\xi ,\alpha )}{d{\tau }_{\lambda }}\\ &&\,=-\,{i}_{\lambda }({\tau }_{\lambda },\xi ,\alpha )\\ &&\,+\displaystyle \frac{{\varpi }_{\lambda }}{4\pi }{\displaystyle \int }_{-\pi /2}^{\pi /2}{\displaystyle \int }_{0}^{2\pi }{p}_{\lambda }(\xi ,\alpha ,\xi ^{\prime} ,\alpha ^{\prime} )\\ &&\,\times {i}_{\lambda }({\tau }_{\lambda },\xi ^{\prime} ,\alpha ^{\prime} )d\alpha ^{\prime} \,d\xi ^{\prime} \\ &&\,+\displaystyle \frac{{\varpi }_{\lambda }}{4\pi }{{\rm{\Phi }}}_{\lambda }{p}_{\lambda }(\xi ,\alpha ,z,A)\exp \left\{-\displaystyle \frac{{\tau }_{\lambda }}{\cos z}\right\},\end{eqnarray} \tag{ 2 }$$

where τ_λ is the optical depth at a certain altitude h, while ξ and α are observational zenith and azimuth angles. The angular distribution of light intensity scattered by an elementary volume of the atmosphere is characterized through p_λ(ξ, α, z, A)/(4 π) (sr⁻¹), where A is the azimuth angle of a light beam entering the atmosphere at its lower interface. Single scattering albedo, ϖ_λ, can be understood as the probability of the survival of a photon in any interaction of the electromagnetic wave with an atmospheric environment. In case of conservative scattering, there is no true absorption, thus ϖ_λ = 1. This approximation is sometimes applied in a characterization of optically thick media, such as clouds (e.g., Melnikova & Vasilyev 2005). In practically important cases, ϖ_λ is below 0.9, but scarcely smaller than 0.5.

The downward diffuse radiation entering the atmosphere at its upper interface is absent, so the only non-zero radiation in the system is generated by the municipality under consideration as outlined, e.g., by Kerola (2006), who expressed the incident radiation field through the radiant intensity I_up,λ(z). However, this can lead to confusion as the RTE (Equation (2)) describes the radiance field. I_up,λ(z) cannot replace i_λ in any part of Equation (2), including boundary conditions (see Liou 2002). The main concern with such a replacement is two-fold:

• 1D radiative transfer equation requires an atmospheric model with plane-parallel layers (Liou 2002; Liang 2004). In this concept, the ground surface is considered to be an infinite plane. However, a light-emitting city or town are finite-sized structures. Radiant intensity I_up,λ(z) of a finite-dimensional flat surface S necessarily approaches zero as the direction of rays escaping the surface becomes parallel to the ground. This is because I_up,λ(z) is a cosine-projected surface integral of i_λ:
  $$\begin{eqnarray}&&{I}_{{up},\lambda }(z)=\cos z{\int }_{S}{i}_{\lambda }(z)\,{dS},\end{eqnarray} \tag{ 3 }$$
  where the integration runs over the light-producing area of a city or town. It is not difficult to show that the second term (direct uplight) in the Garstang's emission function fails to fulfil the above condition (see Section 3 for more details). It is a direct consequence of Equation (3) that radiance (alternatively luminance) of a flat surface emitting in accordance with Garstang's model would exceed any limit. The image of such a surface becomes exceptionally bright when projected onto a cell of a CCD sensor and angle contained by direction normal to the emitting surface and direction toward a camera approaches 90°. Of course, both the saturation threshold and finite size of a photosensitive pixel preclude an uncontrolled increase of radiance impossible.
• The radiant intensity approach to boundary conditions for 1D RTE is also in conflict with the model of infinite emitting media. Let us assume that the optical system with a small finite field of view intercepts the light beams escaping from the ground. The area of sensing would increase proportionally to cos⁻¹ z. Assume S' is the scale of the problem, thus the integral (3) transforms into
  $$\begin{eqnarray}&&{I}_{{up},\lambda }(z)=\cos z{\int }_{0}^{S^{\prime} /\cos z}{i}_{\lambda }(z){dS},\end{eqnarray} \tag{ 4 }$$

which becomes cosine-independent for homogeneously emitting surface, i.e.,

$$\begin{eqnarray}&&{I}_{{up},\lambda }(z)={i}_{\lambda }(z)\cos z{\int }_{0}^{S^{\prime} /\cos z}{dS}={i}_{\lambda }(z)S^{\prime} .\end{eqnarray} \tag{ 5 }$$

In other words, the radiant intensity and radiance of a semi-infinite medium tends to show the same functional dependence on zenith angle, except for a constant S' with dimension of area. Assume a lambertian flat surface extends to infinity in all directions and i_λ is its radiance. Then, it follows from Equation (5) that I_up,λ also becomes constant for any z, meaning that the first cosine-weighted term of the GEF cannot be a proper input to the 1D RTE (see Section 4). In addition, the bulk radiant intensity is measured in W·sr⁻¹, while the broadband radiance has dimension W·m⁻²·sr⁻¹ (the corresponding photometry quantities are cd and cd·m⁻², respectively). A use of radiance and radiant intensity must be clearly distinguished in any light-pollution model, otherwise it opens the door to a whole collection of errors (at least problems involving the dimensional analysis). Although I_up,λ and i_λ are incommensurable physical quantities, they have not received sufficient attention. For instance, both of these functions are characterized through the same variable in Kerola (2006, compare Equations (1) and (12) ibidem).

There is no doubt that natural terrain has a wavy profile rather than being an ideally flat surface. Light sources are also 3D objects with heterogeneous distributions in all directions, including vertical. This is why radiant intensity normally does not approach zero near the horizon. However, emissions toward the horizon are largely suppressed, so indeed the GEF is not the best tool to simulate skyglow from distant sources. There are different demands to the numerical and theoretical methods that are used in light-pollution research, so one of the primary goals is to develop more appropriate mathematical formulae for the CEF.

In general, any real-valued mathematical function that is multiple times differentiable can be expanded into convergent power series. A Taylor series represents a most known approach to the modeling of arbitrary function that satisfies the above conditions, but many other expansions are also used depending on the problem solved. For instance, Kolláth et al. (2016) formulated intensity distribution functions through associated Legendre polynomials, shifted to avoid negative intensities in the calculations. However, a CEF determined through the finite or infinite series is a pure mathematical construction that usually lacks a theoretical foundation. In addition, expansion coefficients suffer from indirect or completely missing correlation to measurable physical properties (such as ground reflectance, total emission upwards, and other parameters characterizing the system). Note that parametrization based on physical quantities rather than arbitrary weighting coefficients is extremely important for a correct understanding of underlying physics and also for subsequent modeling purposes with the aim to simulate the system behavior under various conditions. Such simulation can provide accurate predictions of night sky radiances for arbitrarily combined lighting technologies. The numerical modeling of diffuse light field is advantageous as it can be repeated for various set of input parameters with the intention of improving the planning, designing, and deployment of light-producing zones, taking into account a spread of light pollution into surrounding environment.

Among many others, an analytical formulation of the CEF is a very convenient approach to robust modeling when numerically fast code is absolutely necessary. Direct tasks are about the characterization of light-pollution propagation in a nocturnal environment based on known input data. Typically, it comprises numerical integration of elementary signals through all beam paths taking into account specific atmospheric and environmental conditions. The emission function can be expressed either in the form of an analytical formula or it can be an arbitrary data function. However, data function obtained for a respective locality (or city) does not allow for numerical experimentation with different light sources and their distribution in the territory of a respective city or town. Data function is a product of non-trivial superposition of all signals and characterizes the so-called bulk CEF. Experimentally, it is more or less impossible to identify contributions from different lamp types or to separate the direct emissions from diffusively reflected light. This is why the only experimental monitoring of night light is not sufficient to properly understand all aspects of how CEFs participate in forming the sky radiance and what the fundamentals of signal distortion are before it is recorded at the ground. Nevertheless, it is of great importance to recognize the role of CEF in all of its complexity, as it can help solve the inverse problems.

Inverse problems are about the reconstruction of a CEF using multiangle and/or multispectral optical data harvested at different measuring points. Inverse problems are basically formulated in terms of specific integral equations, where mapping between the CEF and the optical data measured is possible through a kernel of the integral equation. In other words, the inverse kernel operator acting on the measured data produces a solution to the CEF that is, however, very sensitive to errors in the data function. Small perturbations on the measured data may have quite a large response in the solution, thus a special mathematical treatment should be applied to determine sought function—e.g., a stabilizing functional is usually introduced to enforce the smoothness of the solution found, to penalize large norms, and, to penalize steep gradients of CEF. Retrieval of the CEF by solving an inverse radiative transfer problem is advantageous because (i) it is cheap enough, (ii) it provides information on the bulk emission function, and (iii) it is also applicable to routinely measured data such as radiance or luminance patterns (Kocifaj 2014a, 2014b). However, the capabilities of inversion techniques are quite constrained regarding accuracy and information content because data are noisy and sometimes incomplete, and also because the problem is ill-posed (Twomey 2002). This means the solution might become non-continuous, ambiguous, or non-physical (e.g., some values of CEF could be negative, or the CEF retrieved becomes a non-smooth and/or oscillatory function of zenith angle). The situation simplifies when city-emission-function is replaced by an analytical parametrizable formula. In such a case, the inverse problem reduces to a constrained optimization (e.g., Kocifaj & Solano Lamphar 2014). The scaling parameters are altered until the best theoretical match to the experimental radiance data is found. In this process, a set of CEFs is used to compute theoretical radiance under controlled conditions, while the data are compared against the measured ones.

Unfortunately, acceptable analytical parametrizable formula for the CEF is still missing in the light-pollution community.

The radiant intensity function is used to characterize a light source independently of its dimensions, while the RTE is about solving the differential equation for a radiative field that in general varies in space and time. Radiance i_λ is a scalar quantity and can be easily used to compute radiant intensity I_up,λ,  as shown in Equation (3). It has been proved in Section 2 that I_up,λ of a finite-sized light-emitting surface approaches zero at the horizon. This is why the Garstang's model does not work properly at low elevation angles. Garstang's emission function (GEF, Garstang 1986)

$$\begin{eqnarray}&&{I}_{{up},\lambda }(z)={GEF}(z)=\displaystyle \frac{{LP}}{2\pi }[2\,G(1-F)\cos z+0.554F\,{z}^{4}]\end{eqnarray} \tag{ 6 }$$

treats the contributions from reflected light and direct uplight separately: consult the first and second terms at the right-hand side (RHS) of Equation (6). Here, F is the fraction of the light emitted from light sources directly upward, while G characterizes the fraction of the light that is isotropically reflected from the ground. The function is normalized so that the integral of Equation (6) over upper hemisphere equals to

$$\begin{eqnarray}{\displaystyle \int }_{2\pi }{GEF}(z)d{\rm{\Omega }} & = & 2\pi {\displaystyle \int }_{2\pi }{GEF}(z)\,\sin (z){dz}\\ & = & {LP}[G(1-F)+F],\end{eqnarray} \tag{ 7 }$$

where P is the population of a city that produces an output L lumens per head of the population (note that Garstang's model was originally based on photometric quantities). The term (1 − F) characterizes the light emitted from all light sources downward, so the fraction of light directed upward after a diffuse reflection from the ground equals to G(1 − F). It is not difficult to show that placing Equation (6) into Equation (3) causes the radiance of a finite-sized light source to increase proportionally to cos⁻¹(z) as z approaches 90°. Luginbuhl et al. (2009b) introduced an exponential term to reduce the GEF near horizon:

$$\begin{eqnarray}&&{I}_{{up},\lambda }^{L.{et}.{all}}(z)\propto {GEF}(z){10}^{-0.4{E}_{b}/\cos (z)},\end{eqnarray} \tag{ 8 }$$

where E_b is a dimensionless steepness factor that works as the "blocking" extinction with peak effect on ${I}_{{up},\lambda }^{L.{et}.{all}}(z)$ at z = 90°. The above correction to the GEF can guarantee that the radiant intensity approaches zero at the horizon. In fact, there is also another theoretical reason for why I_up,λ(90°) = 0 even if the emissions from finite-sized ground surface are not blocked by obstacles (consult Section 2). Applying the exponential term to the GEF would unfortunately invalidate Lambert's cosine law, which dictates that the radiance of a Lambertian emitter (or reflecting surface) becomes constant at any direction. Namely, the combination of Equations (6) and (8) with Equation (3) generates the radiance

$$\begin{eqnarray}&&{i}^{L.{et}.{all}}=\displaystyle \frac{{LP}}{\pi }\displaystyle \frac{G(1-F)}{S}\,{10}^{-0.4{E}_{b}/\cos (z)},\end{eqnarray} \tag{ 9 }$$

where the second term at the RHS of Equation (6) was omitted because S is assumed to be a Lambertian surface. One can immediately find that ${i}^{L.{et}.{all}}$ is not a Lambertian radiance. Although Equation (8) can apparently predict reasonable levels of night sky radiances, there is a systematic problem with such a model as it violates one of fundamental laws of photometry.

Our current work aims to replace Equations (6) and (8) with a more convenient corrective function,

$$\begin{eqnarray}{I}_{{up}}^{\varphi }(z) & = & {{GEF}}^{\varphi }(z)=\displaystyle \frac{{LP}}{2\pi }\phantom{\rule{}{3ex}}[2G(1-F)\cos z\\ & & +\left.F\,\varphi (R,a){z}^{R}\,\exp \left\{-\displaystyle \frac{a}{\cos z}\right\}\right],\end{eqnarray} \tag{ 10 }$$

proposed to remedy the above-mentioned limitations. The first term at the RHS of Equation (10) models the Lambertian-component of upward emission, while the second term at the RHS of Equation (10) contains exponential correction ${e}^{-a/\cos z}$. In addition, the R-th power of zenith angle aims to simulate different radiant intensity gradations near the horizon. The normalization condition for the second term at the RHS of Equation (10) dictates

$$\begin{eqnarray}&&\tfrac{\varphi (R,a)}{2\pi }{\displaystyle \int }_{2\pi }{z}^{R}{e}^{-a/\cos z}d{\rm{\Omega }}\\ &&\quad =\varphi (R,a){\displaystyle \int }_{0}^{\pi /2}{z}^{R}{e}^{-a/\cos z}\sin z\,{dz}=1.\end{eqnarray} \tag{ 11 }$$

The formula

$$\begin{eqnarray}{I}_{{up}}^{\varphi }(z) & = & {{GEF}}^{\varphi }(z)=\displaystyle \frac{{LP}}{2\pi }\,\left[\phantom{{\frac{{\frac{{M}^{M}}{{M}^{M}}}^{M}}{{\frac{{M}^{M}}{{M}^{M}}}^{M}}}^{M}}\right.\,2G(1-F)\cos z\\ & & +\left.\displaystyle \frac{F\,{z}^{R}\,{e}^{-a/\cos z}}{{\displaystyle \int }_{0}^{\pi /2}{x}^{R}{x}^{-a/\cos z}\sin x\,{dx}}\right],\end{eqnarray} \tag{ 12 }$$

follows from Equation (11), where x is the integral variable. Unfortunately, the integral has no analytical solution; however, a real-time computations are now easily possible with high-performance computers (e.g., 1500 integrals made on a standard PC equipped with Intel I5 processor give the 30 × 50 matrix entries of φ(R, a) within some tens of milliseconds, where the integration step was kept as low as 0.2 degree). The normalization function φ(R, a) on a bounded interval of interest is depicted in Figure 2, showing a simple dependence of φ(R, a) on scaling parameters a and R.

Zoom In Zoom Out Reset image size

The φ (R, a) is increasing the steeper the higher the value of a, but the functional dependence of φ (R, a) on R is generally weak. Gradation along "a-axis" is similar to exponential function, but the analytical fit needs more attention. The function φ (R, a) offers the benefits of analytical differentiation, thus it is possible to adjust the scaling parameters on demands, e.g., by specifying the elevation angle ${z}_{{peak}}^{\varphi }$ for the peak emission. One needs to solve the identity

$$\begin{eqnarray}&&\left(\displaystyle \frac{{R}_{{peak}}}{{z}_{{peak}}\,\sin {z}_{{peak}}}-\displaystyle \frac{{a}_{{peak}}}{{\cos }^{2}{z}_{{peak}}}\right){({z}_{{peak}})}^{{R}_{{peak}}}\\ &&\quad \times \exp \left\{-\displaystyle \frac{{a}_{{peak}}}{\cos {z}_{{peak}}}\right\}=\displaystyle \frac{2\,G(1-F)}{F\,\varphi ({R}_{{peak}},\,{a}_{{peak}})},\end{eqnarray} \tag{ 13 }$$

to retrieve all combinations of [R, a] that satisfy the condition for peak emission. There are several approaches to solve the above nonlinear equation, e.g., bisection method. As a result, a_peak is determined assuming R_peak is selected, or vice versa. We know the GEF uses the value of R = 4 (see Equation (6)). Placing the same value into Equation (13), we can immediately find for F = G = 0.15 that peak emissions at z_peak ≈ 60° require the parameter a_peak to be about 1.0, while for z_peak ≈ 70° we only need a_peak as low as 0.4 (see the detail plots in Figures 3 and 4).

Figure 3 documents the angular behavior of radiant intensity function for uplight fraction F = 0.15 and ground reflectance G = 0.15. The same values were also considered by Garstang in his original work (Garstang 1986). However, the present data rather suggest the lower values of G, see e.g., Aubé et al. (2013), who used G ≈ 0.08 based on asphalt reflectance from the NASA's ASTER spectral library. The values below 0.1 are also applicable to F, specifically for Topolobampo city (Solano Lamphar & Kocifaj 2016), while the measurements performed in Bratislava city have shown that F is from 0.09 to 0.12 depending on spectral band analyzed (Kocifaj et al. 2015). There is actually observed data on the value of F for Flagstaff (Luginbuhl et al. 2009a), where an inventory of outdoor lighting resulted in an estimate of 0.083. To be consistent with these observations, we also made the computations for F = 0.08 in Figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

There is no doubt that the modified ${{GEF}}^{\varphi }$ can model enhanced emissions to low elevation angles and rapid intensity decay near the horizon. Peak intensity can be also controlled by scaling the parameters a and R. However, a few combinations of [a, R] should be avoided due to unrealistic radiant intensity distribution near zenith. In particular, it deals with R below 2–3.

We proceed with another way to tackle the problem by introducing a more convenient and smarter option for the radiant intensity function. The characteristic gradation of radiant intensity from "darker" zenith to "brighter" horizon is treated by the inverse cosine function

$$\begin{eqnarray}{I}_{{up}}^{{\rm{\Psi }}}(z) & = & {{GEF}}^{{\rm{\Psi }}}(z)=\displaystyle \frac{{LP}}{2\pi }\phantom{\rule{}{3ex}}[2G(1-F)\cos z\\ & & +\left.F\,\displaystyle \frac{{\rm{\Psi }}(R,a)}{{\cos }^{R}z}\,\exp \left\{-\displaystyle \frac{a}{\cos z}\right\}\right].\end{eqnarray} \tag{ 14 }$$

Equation (14) takes advantage of analytical formulation of the normalization condition for Ψ (R, a),

$$\begin{eqnarray}&&{\rm{\Psi }}(R,a)=\displaystyle \frac{\,1}{{\int }_{0}^{\pi /2}\tfrac{1}{{\cos }^{R}z}{e}^{-a/\cos z}\sin z\,{dz}}=\displaystyle \frac{{a}^{R-1}}{{\rm{\Gamma }}(R-1,a)},\end{eqnarray} \tag{ 15 }$$

where ${\rm{\Gamma }}(R-1,a)$ is an incomplete Gamma function. Note that the generating function Ψ (R, a) can be formulated by means of finite series in case {$R\in {\mathbb{N}}:R\gt 1$}, i.e., when $R=\{2,3,4,\ldots \}$. Then Equation (15) transforms into a new form:

$$\begin{eqnarray}&&{\rm{\Psi }}(R,a)=\displaystyle \frac{{a}^{R-1}\,{e}^{a}}{(R-2)!\sum _{k=0}^{R-2}\tfrac{{a}^{k}}{k!}}.\end{eqnarray} \tag{ 16 }$$

Analytically, formulae exist for any R ≥ 2, e.g., Ψ(2, a) = a e^a. The complete analytic extension of the normalization function for R = 1 is ${\rm{\Psi }}(1,a)=1/{E}_{1}(a)$, while ${\rm{\Psi }}(0,a)\ ={[{e}^{-a}-a{E}_{1}(a)]}^{-1}$. Here, ${E}_{1}(a)$ is an exponential integral that can be computed as follows:

$$\begin{eqnarray}&&{E}_{1}(a)=-\gamma -\mathrm{ln}(a)-\sum _{n=1}^{\infty }\displaystyle \frac{{(-1)}^{n}{a}^{n}}{n\,\cdot n!},\end{eqnarray} \tag{ 17 }$$

with $\gamma \cong 0.577$ being Euler-Mascheroni constant. Equation (17) is a divergent series for a = 0, but ${\mathrm{lim}}_{a\to 0}a\,{E}_{1}(a)=0$, so Ψ(0, 0) = 1.

The existence of an analytical solution for Ψ (R, a) is a good reason for why Equation (14) can be accepted as a preferred choice in modeling the radiant intensity distributions. Notice that for modeling purposes, it is often necessary to repeat the entire set of calculations with the aim to simulate the system behavior under various conditions—e.g., when light signals from different reflective surfaces are combined with direct emissions from various imperfectly shielded lamps, both being characterized by four independent parameters, namely G, F, R, and a. We recall that the parameter G represents the fraction of the light that is isotropically reflected from the ground, while F is the fraction of the light emitted from light sources directly upwards. Two scaling factors R and a are to control the form of emission function.

A structure of light emissions at angles near the zenith is smooth and lacks a pronounced minima when simulated by GEF^Ψ(z). However, this is not guaranteed with GEF^φ(z), at least for some combinations of [a, R]—compare Figure 2 against Figure 4, or Figure 3 against Figure 5. The generating function Ψ (R, a) even allows for modeling of relatively narrow emission lobes that might be typical for some light sources, e.g., partly or fully unshielded street lighting luminaires or light-emitting sources with non-horizontally oriented shielding.

Zoom In Zoom Out Reset image size

The parameters [a, R] satisfying the condition of peak emission for GEF^Ψ (z) can be found by solving Equation (18):

$$\begin{eqnarray}&&\left({R}_{{peak}}-\displaystyle \frac{{a}_{{peak}}}{\cos {z}_{{peak}}}\right)\exp \left\{-\displaystyle \frac{{a}_{{peak}}}{\cos {z}_{{peak}}}\right\}{\left(\displaystyle \frac{1}{{z}_{{peak}}}\right)}^{1+{R}_{{peak}}}\\ &&\quad =\displaystyle \frac{2\,G(1-F)}{F\,{\rm{\Psi }}({R}_{{peak}},{a}_{{peak}})}.\end{eqnarray} \tag{ 18 }$$

The bisection method is also suitable to solve Equation (18), meaning that R_peak can be determined as a function or a_peak, or vice versa. The peak emissions also depend on G and F. Figure 6 demonstrates that lowering the fraction of upwardly emitted light also suppresses the amplitude of peak emission (compare with Figure 5, where F exceeds that in Figure 6).

Zoom In Zoom Out Reset image size

Heterogeneity of radiating sources is one of the native properties of all ground emissions. This increases the complexity of radiant intensity function when formed as a sum of different signals. In fact, the radiant intensity function is additive. In many situations, it is convenient to expand such complex function in continuous or discrete orthogonal basis functions. Mathematically, it means that the CEF is sought in form of series with expansion coefficients. This approach does not always guarantee that the generated function is positive at all points on an interval of interest. However, one of the basic prerequisites for the CEF is that this is a non-negative distribution function. It is therefore necessary to diverge from a seemingly effective solution in the form of infinite series that might simulate almost everything (without any theoretical base) and search rather for reasonable analytical models.

The Introduction stated that a use of the CEF as an input to the radiative transfer equation is a typical misapprehension and consequently a potential source of systematic errors in the interpretation of computed data. Mass calculations based on the integration of elementary signals over the entire volume of the atmosphere, over a dense grid of azimuth and zenith angles, and over the whole spectrum generally makes the reverse traceability of all inputs difficult or rather impossible. This is why some inconsistencies in the results obtained usually stay unidentified. In addition, model computations also depend on other inputs, such as aerosol microphysics, cloud optics, ground albedo, elevation angle of obscured horizon, etc. It is therefore more than necessary to be careful in using a correct generating function that is either the CEF or the radiance depending on the model setup.

The RTE (Equation (2)) is about obtaining the radiative field generated by a two-component ground radiance: one originating from diffuse reflection, while other one having relation to direct emissions. Both of them are parts of boundary conditions, meaning that, e.g., surface radiance and not the CEF is an input to the RTE (we refer readers to fundamentals of radiative transfer e.g., Equations (3.41) and (2.53) in Zdunkowski et al. 2007). Because of the same reasons, the radiance i_λ(z) has been replaced by ${GEF}(z)/\cos (z)$ in (Kocifaj 2007, compare Equation (27) with Equation (24) taking into account the following convention: z = z₀, ${I}_{{up},\lambda }={GEF}=B$). The theoretical outcomes ibidem are then consistent with alternative derivations by Joseph et al. (1991), who dealt with the number of photons in order to characterize light output.

A wrong belief that the CEF is a generating function for the RTE can result in a number of adverse consequences. Let us illustrate one of them in a simple mathematical experiment. Let Φ_up be the flux emitted by a city upward. Employing Equation (6) in the model of city emission function one can immediately find

$$\begin{eqnarray}{{\rm{\Phi }}}_{{up}}^{{true}} & = & {\displaystyle \int }_{2\pi }{GEF}(z)d{\rm{\Omega }}={\displaystyle \int }_{A=0}^{2\pi }{\displaystyle \int }_{z=0}^{\pi /2}{GEF}(z)\sin z\,{dz}\,{dA}\\ & = & ({LP})\ \{G(1-F)+F\}.\end{eqnarray} \tag{ 19 }$$

We demonstrated in Equation (5) that both the CEF and the radiance show the same behavior when a flat light-emitting surface extends to infinity. In fact, it means the CEF imitates the angular pattern of the radiance. However, an opposite interpretation would lead to misleading results. Namely, if using the above-mentioned linear proportionality in the form of i(z) ∝ I_up,λ (z), one might tend to replace I_up,λ by GEF and then simulate the radiant intensity function mistakenly as a cosine-weighted radiance. Running the integration over an upper hemisphere, we obtain

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{{up}}^{{false}}=({LP})\ \left\{\tfrac{2}{3}G(1-F)+0.233\,F\right\}.\end{eqnarray} \tag{ 20 }$$

The radiative flux obtained from Equation (20) differs from what we found in Equation (19). Let us, e.g., assume that direct uplight is absent. A proper normalization would then result in zenith emissions increased by a factor of 3/2 (compare Equation (20) with Equation (19)). In contrast, the correction factor should be as large as 4.3 if G = 0. Both of the above examples might puzzle the light-pollution community. Although an intentional mistake made in our thought experiment seems impossible to overlook, it might not be true when solving complex radiative transfer problems subject to diverse boundary conditions.

Simulation errors caused by ill-conditioned boundary conditions are examined below using SkyGlow v.5c, a publicly available tool developed by Kocifaj & Kundracik (2015). More specifically, we computed sky luminance and horizontal illuminance at three discrete distances from the city of Bielsko-Biala: 5.4 km, 19.1 km and 40.5 km eastward from a city center (see Figure 7).

Zoom In Zoom Out Reset image size

The radiation emitted from other cities surrounding Bielsko-Biala was ignored to avoid possible bias of computational outcomes. Trial tests were performed for high-pressure sodium (HPS) spectra and the GEF as a generating function, while the fractions of the reflected and an upwardly emitted light were kept consistent with what was introduced by Garstang (1986), i.e., F = G = 0.15. The light output of 1.3 × 10⁸ lumens was obtained as a product of population P and per capita output L (≈750 lm/head). Computations were carried out for a turbid atmosphere characterized by Aerosol Optical Depth (AOD) AOD = 0.3 for the reference wavelength of 500 nm. We adopted Ångström coefficient υ = 1.3 to model AOD as a function of wavelength, i.e., AOD_λ ∝ λ^−υ. The parameter υ often ranges from −0.5 to more than 2 (Cachorro et al. 2000). The angular distribution of the scattered light is dictated by weighted contributions of Rayleigh and Mie scattering functions. The latter one is tuned through the aerosol asymmetry parameter (ASY), which theoretically ranges from −1 to 1. In the numerical runs, we used the value of ASY = 0.9 to characterize the background aerosols with single scattering albedo ϖ_A = SSA = 0.85. Keep in mind that in the light-scattering theories SSA is defined as a fraction of light scattered from the total light experiencing the extinction process. Typically, SSA > 0.5, but quite frequently it is found in the interval from 0.8 to 0.95 (Horvath et al. 2002).

Theoretical luminance patterns that satisfy both Equations (3) and (19) for finite-sized light sources are depicted in Figures 8 and 9 for clear sky and overcast conditions, respectively. In the latter case, a large field of altocumulus clouds with a base at the altitude of 1 km was taken into consideration. All plots are false color images shown on a logarithmic scale. A cosine distortion introduced in Equation (20) generally does not change the fashion in which the scattered light forms sky luminance patterns. This is documented in Figure 10, where overcast luminance distributions for well-posed and ill-posed boundary conditions are subtracted one from each other. The amplitudes of the above differences quickly decay as the distance to the city increases. For an observer situated at the city edge, the largest differences cumulate near zenith position while the horizon luminance data are almost the same. This is because F is as large as 0.15. Note that both parameters, F and G, rule the radiant intensity function, but that only F dominated the emissions to small elevation angles. Therefore, the cosine distortion effect is scarcely apparent in Figures 8 and 9. The largest discrepancies in horizontal illuminance due to the cosine distortion is about 15% for clear sky conditions and even 30% for overcast conditions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Errors induced by replacing the radiance with radiant intensity function in the RTE are more apparent when the uplight fraction is significantly lowered or even completely absent. Figures 11–12 are just to demonstrate the luminance distributions for F = 0, i.e., all photons emitted upwards originate solely from diffuse reflection. In contrast to the situation with F = 0.15 (Figures 8–10), the errors for F = 0 generally increase as the distance to the city increases. Such opposite behavior is due to the diffuse reflection that is now the exclusive source of radiation at all distances (see Figures 11–12). The larger the distance, the lower the elevation angles become important in forming scattered light of a night sky. However, a lowering of elevation angle translates into large errors because of cosine distortion that is maximum near the horizon; see that $\cos {(z)}_{z\to {90}^{^\circ }}\to 0$. This is why errors are almost equal for cloudless and overcast skies at large distances from a city. Notice that the errors obtained in Figures 11–12 are about 2–3 orders of magnitude larger than those in Figure 10, even if the luminance patterns obtained for correct and corrupted boundary conditions do not differ in their structure (compare the three Figures in the top panels of Figures 11–12 with those in the middle therein). The most apparent differences between the luminance maps are, therefore, related to their amplitudes rather than their spatial redistribution. This implies that incorrect computational results might potentially stay unnoticed, most likely due to missing information on total light output from the city. The light emissions still remain poorly quantified and often are estimated based on the population, even if the inventory of light sources is available for a respective city or town.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Here, it is worth noting that the error margin is only acceptable for zenith luminance and only when data are taken at the city edge or at its close surrounding. This is thanks to the cosine factor that approaches unity in zenith. The sky elements characterized by a reduced cosine distortion spreads over a wider area in case of overcast sky (compare the lowermost left diagrams in Figures 11 and 12). Such minima centered in zenith are due to diffuse reflection by cloud base.

Ground-based light sources are conventionally treated as light-emitting flat surfaces (Falchi & Cinzano 2012), while the effects of the Earth curvature on computational results are as low as 2% for distances below 50 km (Shirkey 2011). Radiant intensity function of isotropically emitting finite-sized flat surface decreases proportionally to the cosine of zenith angle, which necessarily results in rapid decline in the number of photons directed towards horizon. We have shown that the same applies to a non-uniformly emitting flat surface, otherwise the radiance of such surface would reach extreme values when the line of sight becomes a tangent to the Earth. It was quite easy to show that the second term of the Garstang's emission function (GEF) for direct uplight fails to fulfil the above requirement. This is why a use of GEF in radiative transfer modeling is at least questionable and requires more attention, especially when skyglow from distant light sources is analyzed.

We developed two alternative city emission functions (see Equations (10) and (14)) that tackle the above problem by introducing corrective functions and strengthening the model considerably: (i) an analytical formula for the CEF has the advantage of an exact mathematical expression throughout its natural domain, meaning that predictions based on such formula are definite and always guarantee a proper understanding of underlying physics; (ii) the CEFs allow for CPU non-intensive computations and thus become well suited for mass modeling studies; (iii) the CEF strengths also include parametrization based on theoretically sound physical quantities that have relation to the city or town under examination; (iv) the CEFs introduced in this paper fulfill the minimum requirements to ensure a successful solution of standard and inverse problems; and (v) Equations (10) and (14) satisfy normalization conditions when integrated over the upward hemisphere, thus guaranteeing the energy conservation law. The last one is sometimes overlooked in favor of mathematical constructions which, however, ignore fundamental principles of optics and rather prefer manipulation with arbitrary functions. For instance, the CEF is sometimes replaced by an infinite or finite polynomial series with scaling parameters that lack any physical meaning, thus potentially resulting in negative values of the CEF.

It is worthwhile to model CEFs by means of Equation (10) or (14), as the conditions for maximum radiant intensity can be adjusted according actual needs. By selecting elevation angle of peak emission all applicable combinations of scaling parameters [R, a] can be found, experimented with, and tested. For instance, by varying R one can reshape primary emission lobe while keeping a constant.

A use of radiant intensity function instead of radiance in the boundary conditions for the RTE can induce errors in computed luminance field. This error dramatically increases as the uplight fraction of emitted radiation approaches zero. However, the varying amplitude of the luminance and not its spatial distribution is the major implication of ill-conceived substitution of the above-mentioned functions. Assuming F = 0.15, the discrepancies between computed horizontal irradiances are only 15%–30% and are exclusively located near the city edges; however, the error margin increases by a factor of 100 or even 1000 for F = 0.0 and the peak errors are found at large distances from a radiating city or town.

Nowadays, a number of cities have improved public lighting in many aspects. Consequently, they achieved a better cleaning up of the night skies through the reduction of direct light emissions upwards. The values of F are now often below 0.1, meaning that the errors originated from the aforementioned misconception could be really large. While the errors in luminance data were virtually invisible for F = 0.15, now they may significantly influence numerical predictions for both cloudless and overcast conditions. There is no doubt that the light outputs from cities or towns are still poorly known. Therefore, the skyglow alterations that span the range of 1–3 orders of luminance magnitude can be mistakenly interpreted as the effect of incorrect estimate of lumens per head of the population rather than the effect of cosine distortion due to ill-posed inputs to the RTE.

Our results appear to be a sufficient motivation for theoreticians and modellers to revise the theoretical and computational models and be concerned with remediation if necessary.

This work was supported by the Slovak Research and Development Agency under contract No: APVV-14-0017. Computational work was supported by the Slovak National Grant Agency VEGA (grant No. 2/0016/16). H.A.S.L. received support from the National Council of Science and Technology of Mexico under the project Cátedras CONACYT # 2723.