A Practical Analytic Model for Daylight
A. J. PreethamPeter Shirley

University of Utah

Brian Smits

www.cs.utah.edu

Figure 1: Left: A rendered image of an outdoor scene with a constant colored sky and no aerial perspective. Right: The same image with aphysically-based sky model and physically-based aerial perspective.

Abstract

1 Introduction

Sunlight and skylight are rarely rendered correctly in computer
graphics. A major reason for this is high computational expense.
Another is that precise atmospheric data is rarely available. We
present an inexpensive analytic model that approximates full spectrum daylight for various atmospheric conditions. These conditions
are parameterized using terms that users can either measure or estimate. We also present an inexpensive analytic model that approximates the effects of atmosphere (aerial perspective). These models
are fielded in a number of conditions and intermediate results verified against standard literature from atmospheric science. Thesemodels are analytic in the sense that they are simple formulas basedon fits to simulated data; no explicit simulation is required to use
them. Our goal is to achieve as much accuracy as possible without
sacrificing usability.

Most realistic rendering research has dealt with indoor scenes. Increased computing power and ubiquitous measured terrain data has
made it feasible to create increasingly realistic images of outdoor
scenes. However, rendering outdoor scenes is not just a matterof scaling up rendering technology originally developed for indoor
scenes. Outdoor scenes differ from indoor scenes in two important
aspects other than geometry: most of their illumination comes directly from the sun and sky; and the distances involved make theeffects of air visible. These effects are manifested as the desaturation and color shift of distant objects and is usually known as aerial
perspective. In this paper we present an efficient closed form approximation that captures the visually salient aspects of these phenomena and is easy to incorporate into a rendering system. Our
approach is motivated by the desire to generate images of real terrain, so we pay attention to maintaining physically-based radiometry, and use input parameters that are readily available to computer
graphics researchers. We feel that current approaches are either far
too general and expensive to be easily used, or include too many
simplifications to generate a sufficiently realistic appearance.
The importance of the phenomena modeled in this paper is emphasized in the psychology and art literature. Psychologists assertthat aerial perspective is a fundamental depth cue that humans use
to estimate distances, and the only absolute depth cue available for
distant unfamiliar objects [10]. Painters use aerial perspective andvariation in sky color in almost all landscape paintings. Da Vinci
devoted an entire chapter of his notebooks to the painting of theseeffects [6]. He stated several characteristics that painters shouldcapture including the whitening of the sky toward the horizon, theincreasing density of aerial perspective toward the ground, and thehue shift of distant objects. Any models of the sky and aerial perspective should make sure they capture these subjective effects. The

CR Categories: I.3.7 [Computing Methodologies ]: Computer
Graphics—Three-Dimensional Graphics and RealismColor, shading, shadowing, and texture
Keywords: sunlight, skylight, aerial perspective, illumination

Permission to make digital or hard copies of all or part of this work forpersonal or classroom use is granted without fee provided that copies
are not made or distributed for profit or commercial advantage and that
copies bear this notice and the full citation on the first page. To copy
otherwise, to republish, to post on servers or to redistribute to lists,requires prior specific permission and/or a fee.
SIGGRAPH 99, Los Angeles, CA USA
Copyright ACM 1999 0-201-48560-5/99/08 . . . $5.00

91

model presented in this paper does capture these effects and is careful to predict their magnitude and color as well. Rendered images
with and without these effects are shown in Figure 1. The image on
the right of Figure 1 was rendered using the techniques from thispaper.
To produce a realistic outdoor image, we need to model the aspects of atmosphere that produce the color of the sky and the effects of aerial perspective. To be most convenient and efficient for
rendering, two formulas are needed. The first should describe thespectral radiance of the sun and sky in a given direction. The second should describe how the spectral radiance of a distant object ischanged as it travels through air to the viewer. Although computer
graphics researchers have captured these effects by explicit modeling, there has so far been no such compact formulas that do not
introduce gross simplifications (e.g., the sky is a uniform color).
While it is possible to directly simulate the appearance of a particular sky given particular detailed conditions, this is inconvenient
because it is a complex and CPU-intensive task, and data for detailed conditions is generally not available. It would be more convenient to have a parameterized formula that takes input data that is
generally available, or is at least possible to estimate. While suchformulas exist for sky luminance, there have not been any for sky
spectral radiance. Given a sky spectral radiance formula, there havebeen no closed-form formulas that account for accurate aerial perspective. This paper presents such a set of formulas that are parameterized by geographic location, time and date, and atmospheric conditions. The formulas are for clear and overcast skies only. While
we do not present results for partially cloudy skies, our clear sky
results should be useful for developing such a model.
Our formulas are parametric fits to data from simulations of thescattering in the atmosphere. These models are analytic in the sensethat they are simple formulas based on fits to simulated data; no explicit simulation is required to use them. We will downplay themechanics of our simulation which is based largely on previous
work in computer graphics. Instead we emphasize a careful discussion on its underlying assumptions and accuracy as well as allmaterial needed to implement our model. In Section 2 we review
previous work on modeling the atmospheric phenomena that are
responsible for the appearance of the sky and objects under natural
illumination. In Section 3 we describe a new model for the spectral
radiance of the sun and sky. In Section 4 we extend that model toinclude the effects of aerial perspective. We present images createdusing these models in Section 5. We discuss limitations of the models and future work in Section 6. All needed formulas and data for
implementing the model are given in the appendices.

primary

secondaryprimary

Figure 2: The earth’s atmosphere receives almost parallel illumination from the sun. This light is scattered into the viewing direction
so that the sky appears to have an intrinsic color. Light may scatterseveral times on the way to the viewer, although primary scattering
typically dominates.

2.1 Atmospheric Phenomena
The visually rich appearance of the sky is due to sunlight scattered by a variety of mechanisms (Figure 2). These mechanisms
are described in detail in the classic book by Minnaert [22], andwith several extensions in the more recent book by Lynch and Livingston [19]. For a clear sky, various types of atmospheric particles
are responsible for the scattering. Because the scattering is not necessarily the same for all light wavelengths, the sky takes on varying
hues.
The details of the scattering depend on what types of particles arein the atmosphere. Rayleigh developed a theory for scattering by air
molecules less than 0:1 in diameter [28]. The crux of the theoryis that the monochromatic optical extinction coefficient varies approximately as ,4 , and this has been verified experimentally. Thismeans that blue light (400nm) is scattered approximately ten times
as much as red light (700nm), which is the usual explanation for
why the sky is blue. Because the short wavelengths in sunlight are
preferentially scattered by the same effect, sunlight tends to become
yellow or orange, especially when low in the sky because more atmosphere is traversed by the sunlight on the way to the viewer.
Although Rayleigh scattering does explain much of the sky’s appearance, scattering from haze is also important. The term hazerefers to an atmosphere that scatters more than molecules alone, butless than fog [21]. Haze is often referred to as a haze aerosol because the extra scattering is due to particles suspended in the molecular gas. These particles are typically much bigger than molecules,
and Mie scattering models the scattering behavior of these particles. Because the haze particles typically scatter more uniformly
than molecules for all wavelengths, haze causes a whitening of thesky. The actual particles come from many sources – volcanic eruptions, forest fires, cosmic bombardment, the oceans – and it is very
difficult to precisely characterize the haze of a given sky. Many researchers, starting with Angstrom, have attempted to describe hazeusing a single heuristic parameter. In the atmospheric sciences literature, the parameter turbidity is used [21].
Turbidity is a measure of the fraction of scattering due to haze as
opposed to molecules. This is a convenient quantity because it canbe estimated based on visibility of distant objects. More formally,turbidity T is the ratio of the optical thickness of the haze atmosphere (haze particles and molecules) to the optical thickness of theatmosphere with molecules alone:

2 Background
Many applications use estimates of energy levels of skylight andsunlight to aid in simulation. For this reason modeling skylight has
been studied in many fields over several decades. For rendering, weneed a function of the form:
sky : (direction, place, date/time, conditions) ! spectral radiance :Here “place” is the geographic coordinates (e.g., latitude/longitude)of the viewer. Such a formula would allow a renderer to query thesky for a color in a specific direction for either display or illumination computation. The spectral radiance of the sun should in principle be given by a compatible formula.
In this section we review previous approaches to generating formulas related to the sky function above. We will see that no efficient
formula of the form of sky has previously appeared, but that manytechniques are available that bring us close to that result.

T = (tm + th )=tm;

where tm is the vertical optical thickness of the molecular atmosphere, and th is the vertical optical thickness of the haze
R s atmosphere. Optical thickness for a given path is given by 0 (x)dx92pure air256

exceptionally clear
128
64

very clear

32
16
Rm
(km) 8

clear

θ

θs

light haze

4

haze
2

W

thin fog

1

γ

v
N

0.5

1

2

4

8
Turbidity

16

32

64

S

Figure 3: Meteorological range Rm for various turbidity values.
Values computed from source data in McCartney [21]φsEFigure 4: The coordinates for specifying the sun position and the
direction v on the sky dome.

where (x) is the scattering coefficient (fraction scattered per meter of length traveled) which may vary along the path. Several
other definitions of turbidity are used in various fields, so some caremust be taken when using reported turbidity values. Since turbidity
varies with wavelength, its value at 550nm is used for optical applications [21]. Turbidity can also be estimated using meteorologic
range, as is shown in Figure 3. Meteorological range Rm is the
distance under daylight conditions at which the apparent contrast
between a black target and its background (horizon sky) becomesequal to the threshold contrast ( = 0:02) of an observer, and it
roughly corresponds to the distance to most distant discernible geographic feature.
Although turbidity is a great simplification of the true natureof the atmosphere, atmospheric scientists have found it a practicalmeasure of great utility. Because it does not require complex instrumentation to estimate turbidity, it is particularly well-suited for
application in graphics, and we use it to characterize atmosphericconditions throughout the rest of this paper.

later adopted as a standard by the CIE [4]:
,





e,3 + 0:45 cos2 1 , e,0:32= cos 
Yz (0:91 + 10e,3s + 0:45 cos2  ) (1 , e,0:32) ; (1)swhere Yz is the luminance at the zenith, and the geometric terms
are defined in Figure 4. The zenith luminance Yz can be found:

0 91 + 10

in tables [17], or can be based on formulas parameterized by sun
position and turbidity [4, 16].
In computer graphics, the CIE luminance formula has been usedby several researchers (e.g., [23, 30]). To get spectral data for values returned by the CIE luminance formula, Takagi et al. inferredassociated color temperature with luminance levels using empirical
data for Japanese skies, and used this color temperature to generate a standard daylight spectrum [29]. In the Radiance system theluminance is multiplied by a unit luminance spectral curve that is
approximately the average sky color [31].
Moon and Spencer developed a formula for the luminance distribution of overcast skies which was later adopted by the CIE [4]:2.2 Atmospheric Measurements and Simulation
One way to develop a sky model is to use measured or simulated
data directly. The CIE organized the International Daylight Measurement Program (IDMP) to collect worldwide information on
daylight availability. Several other efforts have collected measured
data that can be used directly. The data sources do not include spectral radiance measurements, so they are not directly useful for our
purposes. Ineichen et al. surveyed these data sources and comparedthem to analytic sky luminance models [12]. They did find that existing sky luminance models are reasonably predictive for real skies
in a variety of locations around the world.
Various computer graphics researchers have simulated atmospheric effects. Blinn simulated scattering for clouds and dusty
surfaces to generate their appearance [1]. Klassen used a planar layer atmospheric model and single scattering to simulate skycolor [18]. Kaneda et al. employed a similar simulation using a
spherical atmosphere with air density changing exponentially withaltitude [15]. Nishita et al. extended this to multiple scattering
to display sky color [25] and also simulated atmospheric scattering to display earth and atmosphere from space [26]. All of thesemethods require a lengthy simulation for a given sky condition, butthey have the advantage of working with arbitrarily complex atmospheric conditions.

Yz (1 + 2 cos ) =3;(2)

There are various more complicated formulas for overcast sky luminance, but they vary only subtly from Equation 2 [17]. The zenith
values for luminance of overcast skies can be found from tables[17] or from analytic results adopted by CIE [4].
To gain efficiency over brute-force simulations, while retaining
the efficiency of the CIE representation, researchers have used basis
functions on the hemisphere to fit simulation data. Dobashi et al.
used a series of Legendre basis functions for specific sky data [7].
These basis functions can be used to fit any sky data, and does not
supply a specific analytic sky model. Rather, it provides a representation and a fitting methodology for some arbitrary data set. Thesebasis functions have the advantage of being orthogonal, but havethe associated property that care must be taken to keep the approximation nonnegative everywhere. Because these basis functions are
not tailored specifically for sky distributions, many terms might beneeded in practice. Our work differs from that of Dobashi et al.
in the choice of basis functions. More importantly, we supply theparameters resulting from our simulations, so the formulas in thispaper can be used directly.
Nimeroff et al. used steerable basis functions to fit various sky
luminance models including the CIE clear sky model [24]. They
demonstrated that the steerable property yielded great advantage in
rendering applications. They used approximately ten basis functions for their examples.
Brunger used the SKYSCAN data to devise a sky radiancemodel [2]. His model represented the sky radiance distribution as a

2.3 Analytic Sky Models
For simpler sky conditions, various researchers have proposed parametric models for the sky. Pokrowski proposed a formula for sky
luminance (no wavelength information) based on theory and sky
measurements. Kittler improved this luminance formula which was93sunlight
inscatter

light
from
objectresults for fog effects where density variation of fog was exponential [15] which is a special case of Max’s layered fog model. However, it is not possible to analytically solve for the extended case of
air combined with haze.
Several researchers have simulated aerial perspective using explicit modeling [15, 18]. This in fact is just a particular instance ofgeneral light scattering simulation. While such techniques have theadvantage of working on arbitrary atmospheric conditions, they arealso computationally expensive.
Ward-Larson has implemented a simpler version of aerial perspective in the Radiance system [30]. He assumes a constant ambient illumination that does not vary with viewing direction. This
produces an efficient global approximation to aerial perspective, butdoes not allow the changes in intensity and hue effects for changing
viewer or sun position.
In Ebert et al., the aerial perspective effect is modeled througha simulation of single Rayleigh scattering [8]. The color of distantmountains is a linear combination of the mountain color and skycolor whose weighting varies with distance. They include a sophisticated discussion of how to numerically integrate the resulting expressions. Although they restrict themselves to pure air (turbidity1), their techniques could easily be extended to include haze because they use numeric techniques. The only shortcoming of their
method is that the quadrature they perform is intrinsically costly,
although they minimize that cost as much as possible.

skylight
inscatter

skylight
inscatteroutscatter

Figure 5: The color of a distant object changes as the viewer moves
away from the object. Some light is removed by out-scattering, andsome is added by in-scattering.

composition of two components, one depending on viewing angle
from zenith and the other on scattering angle. An analytic radiancemodel is very useful for illumination engineers for energy calculations, but what the graphics community needs is a spectral radiancemodel and not a radiance model.
Perez et al. developed a five parameter model to describe the sky
luminance distribution [27]. Each parameter has a specific physicaleffect on the sky distribution. The parameters relate to (a) darkening or brightening of the horizon, (b) luminance gradient near thehorizon, (c) relative intensity of the circumsolar region, (d) widthof the circumsolar region and (e) relative backscattered light. Thesebasis functions can be fit to any data, and are designed to capture
the overall features of sky distributions without ringing or a data
explosion. Perez et al.’s model is given by:

3 Sunlight and Skylight
This section describes our formulas for the spectral radiance of thesun and the sky. The input to the formulas is sun position and turbidity. Sun position can be computed from latitude, longitude, time,and date using formulas given in the Appendix. We assume the U.S.
Standard Atmosphere for our simulations. We use Elterman’s data
for the density profile for haze up to 32km [9].

F (; ) = (1 + AeB=  )(1 + CeD + E cos ); (3)where A; B; C; D and E are the distribution coefficients and and are the angles shown in Figure 4. The luminance Y for sky incos

2

3.1 Sunlight
For sunlight we use the sun’s spectral radiance outside the earth’s
atmosphere, which is given in the Appendix. To determine howmuch light reaches the earth’s surface we need to compute the fraction removed by scattering and absorption in the atmosphere. Sunlight is scattered by molecular and dust particles and absorbed byozone, mixed gases and water vapor. In what order, this attenuation
takes place does not matter because attenuation is multiplicative
and thus commutative. Iqbal gives direct radiation attenuation coefficients for the various atmospheric constituents [13], so we cancompute the total attenuation coefficient if we know the accumulated densities along the illumination path.
The sun’s extraterrestrial spectral radiance is multiplied with thespectral attenuation due to each atmospheric constituents to give us
the sun’s spectral radiance at earth’s surface. Transmissivity due tothese constituents are given in the Appendix. To test our formulas
we compared the sunlight results at turbidity two with measured
values of sun chromaticity and luminance given by Wyszecki andStiles [32]. These numbers were with two percent of each other.

any viewing direction depends on the distribution function and thezenith luminance and is given by

Y = Yz F (; )=F (0; s):

(4)

The Perez model is similar to the CIE model, but has been found tobe slightly more accurate if the parameters A through E are chosen
wisely [12]. The Perez formula has been used in graphics withslight modification by Yu et al. [33].
What would be most convenient for computer graphics applications is a spectral radiance analog of Equation 1 that captures thehue variations suggested by real skies and full simulations. Such aform will be introduced in Section 3.

2.4 Aerial Perspective
A sky model is useful for both direct display and illuminating the
ground. However, it is not directly applicable to how the atmosphere changes the appearance of distant objects (Figure 5). Unlike
a sky model, aerial perspective effects cannot be stored in a simple
function or precomputed table because they vary with distance andorientation. A subtlety of aerial perspective is that it can cause color
shifts in any direction. Typically, it causes a blue-shift, but when theviewing direction is near the sun, it can cause a yellow-shift. It is
hard to predict such color shifts without a physically-based model.
Max presented an analytic single scattering model for light scattering through haze with uniform density and a generalized result
for the case of layered fog [20]. Kaneda et al. presented analytical

3.2 Skylight Model
Skylight is much more complicated to model than sunlight. Given amodel for the composition of the atmosphere, we can run a simulation using the methods of previous researchers. However, we would
then have the data for only one turbidity and sun position. What we
do is compute the sky spectral radiance function for a variety of sun
positions and turbidities, and then fit a parametric function. Basic issues that must be addressed are the assumptions used for thesimulation, and the parametric representation we use to fit the data.

94

L0

for the two gas constituents. This does not make the aerial perspective equations solvable analytically, but it does make them tractableenough to be approximated accurately. This approximation will be
described for the rest of this section. We assume that the earth isflat, which is a reasonable assumption for viewers on the ground.
Aerial perspective results when the light L0 from a distant object
is attenuated on the way to the viewer. In addition, light from thesun and sky can be scattered towards the viewer. This is shownin Figure 6. If  is the extinction factor as L0 travels a distances to reach the eye, and Lin is the in-scattered light, then L(s) =L0  + Lin.
Both the extinction factor  and the in-scattered light Lin are a
result of the scattering properties of the different particles in theatmosphere. Because the scattering coefficients of particles is proportional to the density of particles, the scattering coefficients also
decrease exponentially with height. Thus, (h) = 0 e, h , where0
is the value of scattering coefficient at earth’s surface and is
the exponential decay constant. In our case, h is a function of the
distance from the viewer, as shown in Figure 6, and can be represented as h(x) = h0 + x cos  . We can now write the expression
for asinscatter
ω
φ
L(s)

θ

x
h(x)
h0

Figure 6: The variables used to compute aerial perspective.

For the simulation we used the method of Nishita et al. [25]. The
earth was assumed flat for zenith angles less than seventy degrees
and spherical for other angles. This allowed several terms to beevaluated analytically for the smaller angles. Third and higher order scattering terms were ignored as their contribution to skylightis not significant. Reflectance of light from the earth’s surface wasalso ignored. This simulation was run for twelve sun positions andfive different turbidities (2 through 6). The spectral radiance wascomputed for 343 directions in a sky dome for each of these combinations. Because the amount of computation required was large(about 600 CPU hours in all) a number of careful optimizations
were employed to make the computation feasible such as an aggressive use of lookup tables and adaptive sampling of directions.
For our parametric formula for luminance we use Perez et al.’s
formulation (Equation 4). This formulation has been battle-tested
and has few enough variable that the optimization stage of the fitting process is likely to converge. We use this in preference to the
CIE model because it has a slightly more general form and can thus
capture more features of the simulated data. To account for spectral
variation, we also fit chromatic variables. We found Perez’s formulation to be a poor way to represent the CIE X and Z variables, but
the chromaticities x and y are well represented with this five parameter model. The functions were fit using Levenberg-Marquardt
non-linear least squares method in MATLab [11]. A linear function
in turbidity was obtained to describe the five parameters for Y , x
and y . The zenith values for Y , x and y were also fit across different
sun positions and turbidities.
Chromaticity values x and y are similarly behaved and are givenby the same model. Thus,hx( ( )) =

u(x);

(5)

where u(x) = e, (h0 +x cos ) is the ratio of density at point x to
the density at earth’s surface. The other scattering term we needmust describe the fraction of light scattered into the viewing direction (; ) from a solid angle ! . This is commonly denoted(!; ; ; h). Using the same trick as for (h), it can be rewrittenas!; ; ; h(x)) = 0(!; ; )u(x):((6)

4.1 Extinction Factor
The extinction factor  can be determined directly given our assumptions of an exponential density of particles. Attenuation oflight due to particles withover a disR s total scattering coefficienttance s is given by e, 0 dx . Using equation 5 and integrating, we
haveR = e, s u(x)dx = e, e, h
0

0

0

0 (1

,e, cos s )
cos



:

For convenience, we make the substitutions K = , cos  and
H = e, h0 , allowing the extinction factor to be neatly written
as e,K (H ,u(s)) for a single type of particle.
Atmosphere contains both molecules and haze, both of which
scatter light. The scattering properties of a particle is independentof the presence of other particles and therefore the total attenuation
due to the presence of two types of particles is equal to the productof the attenuation by each individual particles. This means the total
extinction due to both these particles is0x = xz FF(0(;;  )) ; and y = yz FF(0(;;  )) ;s

0s

where F is given by Equation 3 with different values of(A; B; C; D; E ) for x and y . The distribution coefficients andzenith values for luminance Y , and chromaticities x and y are givenin the Appendix. The luminance Y and chromaticities x and y canbe converted to spectral radiance on the fly using the CIE daylight
curve method described in the Appendix.

 = e,K (H ,u (s)) e,K (H ,u (s)) ;
1

1

1

2

2

2

(7)

where the subscript “1” denotes haze particles and the subscript “2”denotes molecular particles.

4 Aerial Perspective Model
Unlike the sun and sky, aerial perspective cannot be precomputed
for a given rendering. At every pixel it is a complex integral thatmust be evaluated numerically. Because we want to capture thesubjective hue and intensity effects of aerial perspective we must
preserve a reasonable degree of accuracy. But to make the problem
tractable we assume a slightly simpler atmospheric model than we
did for skylight: we approximate the density of the particles as exponential with respect to height. The rate of decrease is different

4.2 Light Scattered into Viewing Ray
At every point on the ray, light from the sun/sky is scattered into
the viewing direction. Let Ls (! ) denote the spectral radiance ofsun and sky in the direction ! . We can assume that the spectral
radiance from sun and sky, Ls (! ) does not depend on altitude because the viewer and source are close to the earth’s surface. Let
S (; ; x) be the term to denote the light scattered into the viewing95direction (; ) at point x. Using the angular scattering coefficientfrom equation 6, we can express the light scattered into the viewing
direction at x as

S (; ; x)Z

are determined such that g (v ) interpolates the position and slope ofthe endpoints of f (x). The resulting integral,

I1 = ,

Ls (!) (!; ; ; h)d!

=Z

Ls (!) 0(!; ; )u(x)d!
0= S (; )u(x);Rwhere S 0 (; ) = Ls (! ) 0 (!; ; )d! is the light scattered into
the viewing direction (; ) at ground level.
If we denote attenuation (equation 7) from 0 to x along viewing ray as  (0::x) then, the total light scattered into the viewing0

S (; ; x) (0::x)dx =

Z s
0

=

Z s

Our model was implemented in a C++ path tracer [14] that accepts
30m digital elevation data. All images are of a constant albedo terrain skin of approximately 4000km2 . The 30m resolution cells visible in the foreground of the images give an idea of scale. The implementation of the model was not carefully optimized, and slowed
down the program by approximately a factor of two on a MIPS
R10000 processor. The images are 1000 by 750 pixels and were
run with 16 samples per pixel.
Figure 7 shows the same landscape at different times of day andturbidities for a viewer looking west. Note that near sunset, there
is much warm light visible in the aerial perspective for the higher
turbidities. This is as expected because the high concentrations ofaerosols present at high turbidities tend to forward scatter the sunlight which has had much of the blue removed by the thick atmosphere for shallow sun angles.
We used a high value for turbidity for the last picture of Figure 7.
For these high values, we would typically expect an overcast sky forsuch high turbidities, and this is shown in the figure using the CIE
overcast sky luminance and a flat spectral curve. For intermediate
turbidities our model and the overcast model should be interpolated
between as recommended for the CIE luminance models. These
unusual conditions are the “hazy, hot, and humid” weather familiar
to the inland plains.
Figure 8 shows a comparison between the model used by WardLarson in the Radiance package and our model for a summer sky
a half hour before sunset with turbidity 6. Our implementation of
Ward-Larson’s model uses the correct luminance but the relative
spectral curve of the zenith. It correctly sets the attenuation at one
kilometer and uses an exponential interpolant elsewhere. For inscattering it uses the product of the zenith spectral radiance and thecomplement of the attenuation factor. This is our best estimate forsetting the “ambient” in-scattering term suggested by Ward-Larson.
We could certainly hand-tune this in-scattering term to produce better results for one view, but it would cause problems for other views
because Ward-Larson’s model does not take view direction into account. Note that for our model at sunset the east view always has
a blue-shift in the hue (because of backward Rayleigh scattering),
and a yellowish shift for west views depending upon turbidity. Thiseffect is not possible to achieve with a model that does not vary withdirection.

S 0(; )u(x) (0::x)dx:

S20(; )u2 (x) (0::x)dx000= S1 (; )I1 + S2 (; )I2 ;(8)
Rs
where Ii = 0 ui (x) (0::x)dx. A table of S10 (; ) and S20 (; )for different  and  can be precomputed thus avoiding expensive

computation for every pixel.
We show how to solve I1 in this paper; the solution for I2 isanalogous. First we expand  (0::x) and examine the results.

I1 =

Z s
0

u1 (x)e,K (H ,u (x)) e,K (H ,u (x)) dx:
1

1

1

2

2

2

(9)

If j s cos j  1 which would happen when the viewing ray isclose to horizon or the distances considered are small, the term

, x cos e,K(H ,u(x)) = e, H ,e
1

I1

=

Z s
0

=

Z s
0

=

cos



 e, Hx. Thus

u1 (x)e,K (H ,u (x)) e,K (H ,u (x)) dx
1

1

1

2

2

2

e,H e, x cos  e, H xe, H xdx
1

1

1

1

2

2

,( 1 cos + 1 H1 + 2 H2 )se,H 1 , ecos  + H + H :
1
1 1
2 2
1

(10)

Otherwise, two different approaches could be taken to solving
these integrals. The simplest and most accurate method of calculating the integrals I1 and I2 are by numerical integration techniques.
This is too expensive for the model to remain practical. We makeapproximations to the expressions above to present the results inclosed form. In Equation 9 we make the following substitution, v =u1(x) = e, 1 (h0 +x cos ). Therefore, dv = , 1 cos u1 (x)dx.
We now have

I1 = ,

1
1 cos

Z u1 (s) u (0) e
1(11)

5 Results

S10(; )u1 (x) (0::x)dx +

Z0 s

1

the Appendix.

Since there are two kinds of particles (haze and molecules), the total
light scattered into viewing direction is:

Lin,K1 (H1 ,v) g (v )dv;

 u (0) e

can be integrated by parts, leaving an analytic approximation for

direction for a single type of particle is:
Z s

1 cos

Z u1 (s)I1. This result and the coefficients for the polynomial are given in

=

Lin =

1

6 Conclusions and Future Work
We have presented a reasonably accurate analytic model of skylightthat is relatively easy to use. It captures the effects of different atmospheric conditions and times of day. In the same spirit, we havepresented a model for aerial perspective. The use of both models greatly enhances the realism of outdoor rendering with minimal
performance penalties, which may allow widespread use of theseeffects for rendering.
Our models use uniform (exponential or nearly exponential) density distributions of particles. These assumptions do not hold forcloudy (or partly cloudy) skies. They also do not hold for fog or,K1 (H1 ,v) e,K2 (H2 ,u2 (x)) dv:

We replace the term f (x) = e,K2 (H2 ,u2 (x)) with a Hermite cubic
polynomial g (v ) = av 3 + bv 2 + cv + d so that I1 is integrable inclosed form. The coefficients a, b, c and d for the cubic equivalent96

Figure 8: Top: the CIE clear sky model using constant chromaticity
coordinates and Ward-Larson’s aerial perspective approximationfor east viewing directions and the same viewpoint. Bottom: thenew model. Note the change in hue for different parts of the sky forthe new model.

the effects of localized pollution sources and inversion effects that
often occur near some large cities. In these cases the density distribution of particles is much more complicated than in our model.
In these cases, our model can be used as boundary conditions formore complex simulations.

7 Acknowledgments
Thanks to Sumant Pattanaik for the original CIE sky model implementation, to Chuck Hansen and the Los Alamos ACL lab forenabling computer runs, to Bill Thompson and Simon Premoze forsupplying digital terrain data, and to the Utah SGI Visual Supercomputing Center for thousands of hours of CPU time. This work
was supported by NSF awards 9731859 and 9720192, and the NSF
STC for Computer Graphics & Visualization.

A Appendix
Although much of the data in this appendix is available in the literature, it is not in sources readily accessible to most graphics professionals. The information here should allow users to implement ourmodel without sources other than this paper.

Figure 7: The new model looking west at different times/turbidities.
Top to bottom: morning, tubididity 2; evening, turbidity 2; evening,turbidity 6; overcast, turbidity 10.

97

n

A.1 Transmittance expressions for atmosphericconstituents
Simple results are given describing the attenuation of direct radiation by various atmospheric constituents using the data given by
Iqbal [13]. The formulas permit atmospheric parameters such asozone layer thickness, precipitable water vapor and turbidity to bevaried independently. These results are used in the computation ofsunlight received at earth’s surface.
Relative optical mass m is given by the following approximation,
where sun angle s is in radians:s Transmittance due to Rayleigh scattering of air molecules (r; ),Angstrom’s turbidity formula for aerosol (a; ), transmittance due
to ozone absorption (o; ), transmittance due to mixed gases absorption (g; ) and transmittance due to water vapor absorption(wa;) are given by:

r;
a;
o;
g;

wa;

=
=
=
=
=

450
4.193
3.319
2.868
2.527
1.129
0.33730.16820.12750.091110.086520.08767
0.10240.1368

550
4.177
3.329
2.878
2.536
1.138
0.34330.17300.13200.09591
0.09179
0.09352
0.10950.1447

650
4.147
3.335
2.883
2.542
1.142
0.34670.17570.1346
0.09871
0.09488
0.09697
0.11370.1495

850
4.072
3.339
2.888
2.547
1.147
0.35020.17850.13730.10167
0.098160.100650.11820.1566

Table 1: Scattering term  () for Mie scattering.

m = cos  + 0:15  (93:1885 ,  180 ),1:253 :
s

400
4.192
3.311
2.860
2.518
1.122
0.33240.16440.12390.08734
0.082420.08313
0.097010.1307

1
4
7
10
30
60
80
90
110
120
130
150
180

T )( , 2s ).
where  = ( 49 , 120
Zenith chromaticity (xz ; yz ):,4:08me,0:008735
;
, m
,

e
;
,
k
lm
o;
e
;
,
1:41kg; m=(1+118:93kg;m)0:45
e;,
0:2385kwa; wm=(1+20:07kwa;wm)0:45
e
;



2



2

xz = Tyz = T

where is Angstrom’s turbidity coefficient, is the wavelength exponent, ko; is the attenuation coefficient for ozone absorption, l isthe amount of ozone in m at NTP, kg; is the attenuation coefficientof mixed gases absorption, kwa; is the attenuation coefficient ofwater vapor absorption, w is the precipitable water vapor in m and is the wavelength in m. The coefficient varies with turbidity
T and is approximately given by 0:04608T , 0:04586. As originally suggested by Angstrom, we use = 1:3. A value of 0.0035m
for l and 0.02m for w is commonly used.
The spectrums ko; , kg; and kwa; are found in Table 2.

,0:0037 0:0021
0:0638,0:0320,0:2120 0:0605

0 000 6 ss2 7
7
0 0039 5 6
4 s5
0 25891:,::,0:0061 0:0032
0:0897,0:0415,0:2676 0:0667

0 000 6 ss2 7
7
0 0052 5 6
4 s5
0 26691TT

0 0028
1 4 0 0421
0 1535







2 33

:
,:
:

2

0 0017
1 4 0 0290
0 1169



2

:
:
:

:
:
:

3

3





2 33

A.3 Scattering Coefficients
In scattering theory, the angular scattering coefficient and thetotal scattering coefficient determine how the light is scatteredby particles. For our work Rayleigh scattering is used for gas
molecules and Mie scattering theory is used for haze particles.
Here we give the scattering coefficients for gas molecules and haze.
Notice that the total scattering coefficient is the integral Rof angular
scattering coefficient in all directions, for example =()dw .
For an elaborate discussion on scattering, see [21, 28].

A.2 Skylight Distribution Coefficients and Zenith
Values

The angular and total scattering coefficients for Rayleigh scattering for molecules are:

The distribution coefficients vary with turbidity and the zenith values are functions of turbidity and sun position.
The coefficients for the Y; x and y distribution functions are:

AY 3 2 0:1787 ,1:46303
0:42757  
6BY 7 6,0:3554
7 T
7 6
6
5:32517
6 CY 7 = 6,0:0227
4D 5 4 0:1206 ,2:57715 1
Y,0:0670 0:3703
EY
2 3 2
Ax
,0:0193 ,0:25923
6 Bx 7 6,0:0665
0:0008 7  
6 7 6
7 T
0:2125 7
6 Cx 7 = 6,0:0004
4D 5 4,0:0641 ,0:89895 1
x
Ex,0:0033 0:0452
2 3 2
Ay,0:0167 ,0:26083
6By 7 6,0:0950
0:0092 7  
6 7 6
7 T
0:2102 7
6 Cy 7 = 6,0:0079
4D 5 4,0:0441 ,1:65375 1
y,0:0109 0:0529
Ey
2

m ( )=m=

2(n2 , 1)2 ( 6 + 3pn )(1 + cos2 );
2N4
6 , 7pn3228 (n , 1) 6 + 3pn
(
);
3N4
6 , 7pnwhere n is refractive index of air and is 1.0003 in the visible spectrum, N is number of molecules per unit volume and is 2:545x1025
for air at 228.15K and 1013mb, pn is the depolarization factor and0.035 is considered standard for air.
The angular and total scattering coefficients for Mie scatteringfor haze are:

p ()

=

p

=

c 2 )v,2 12 ();2 v,2
0:434c ( ) K;:

0 434 (

where c is the concentration factor that varies with turbidity T andis (0:6544T , 0:6510)x10,16 and v is Junge’s exponent with avalue of 4 for the sky model. A table for  (; ) for v = 4 (Source:[3]) is given in Table 1, and the spectrum for K is given in Table 2.

Absolute value of zenith luminance in K cd m,2 :Yz = (4:0453T , 4:9710) tan  , 0:2155T + 2:4192;
98


nm)

A.4 Aerial Perspective Formulas(

The expression for aerial perspective is L(s) = L0  + Lin , where
L0 is the spectral radiance of the distant object. The extinction factor  is given by equation 7. The light scattered into the ray is handled differently depending upon the viewing angle  and distance s.
For j s cos j  1 we use equations 10 and 8. Otherwise we needto integrate the expression from equation 11. First the integration:

I1

380
390
400
410
420
430
440
450
460
470
480
490
500
510
520
530
540
550
560
570
580
590
600
610
620
630
640
650
660
670
680
690
700
710
720
730
740
750
760
770
780

Z u (s)1e,K1 (H1 ,v) g(v)dv
 u1 (0)1g(v) , g0 (v) + g00 (v) , g000 (v) )]u1 (s), cos[e,K1 (H1 ,v ) (

K1
K12
K13
K14 H1
1
1, cos  ((e,K1 (H1 ,u1 (s)) ) 
1
0
00
000
g
(u1 (s))(, g (uK12(s)) + g (Ku13(s)) , g (Ku14 (s)) ) ,K1
11
1
0 (H1 ) g00 (H1 ) g000 (H1 )
g
(H1 )
g
(
, K 2 + K 3 , K 4 ))
K1
1
1
1

,

=
=
=

1

1 cos

The values of a, b, c, and d for the function g (v ) = av 3 + bv 2 +cv + d to approximate f (x) = e,K2 (H2 ,u2 (x)) where v = u1(x)are determined by the solution to the following system of linear
equations:2

H13
H12 2 H1
6 u1 (s)3
u
(s)
u1(s)
1
6
4 3H122 H1
1
2
3u1 (s)2u1 (s)
1
:

32 3

a
b
c
d

2

3

1
1
17
7 6 ( )776=
4 5 4 (0)5
05
( )
0

The values for the exponential decay constant,1 and molecules = 0:1136 km,1 .

fs
f 00
f s

are:

0 8333 km

haze =

A.5 Converting Tristimulous Values to Spectral
Radiance

=

M2

=

S0

S1

S2

0.650
0.653
0.656
0.658
0.661
0.662
0.663
0.666
0.667
0.669
0.670
0.671
0.672
0.673
0.674
0.676
0.677
0.678
0.679
0.679
0.680
0.681
0.682
0.682
0.683
0.684
0.684
0.685
0.685
0.685
0.686
0.686
0.687
0.687
0.688
0.688
0.689
0.689
0.689
0.689
0.689

63.4
65.8
94.8
104.8
105.9
96.8
113.9
125.6
125.5
121.3
121.3
113.5
113.1
110.8
106.5
108.8
105.3
104.4
100
96
95.1
89.1
90.5
90.3
88.4
84
85.1
81.9
82.6
84.9
81.3
71.9
74.3
76.4
63.3
71.7
77
65.2
47.7
68.6
65

38.5
35
43.4
46.3
43.9
37.1
36.7
35.9
32.6
27.9
24.3
20.1
16.2
13.2
8.6
6.1
4.2
1.9
0
-1.6
-3.5
-3.5
-5.8
-7.2
-8.6
-9.5
-10.9
-10.7
-12
-14
-13.6
-12
-13.3
-12.9
-10.6
-11.6
-12.2
-10.2
-7.8
-11.2
-10.4

3
1.2
-1.1
-0.5
-0.7
-1.2
-2.6
-2.9
-2.8
-2.6
-2.6
-1.8
-1.5
-1.3
-1.2
-1
-0.5
-0.3
0
0.2
0.5
2.1
3.2
4.1
4.7
5.1
6.7
7.3
8.6
9.8
10.2
8.3
9.6
8.5
7
7.6
8
6.7
5.2
7.4
6.8

Sun
rad.
16559
16233.7
21127.5
25888.2
25829.1
24232.3
26760.5
29658.330545.430057.5
30663.728830.4
28712.1
27825
27100.627233.6
26361.3
25503.825060.2
25311.625355.925134.2
24631.5
24173.2
23685.323212.1
22827.722339.821970.2
21526.7
21097.920728.320240.4
19870.8
19427.2
19072.4
18628.9
18259.2
-

ko

kwa

kg

0.3
0.6
0.9
1.4
2.1
3
4
4.8
6.3
7.5
8.5
10.3
12
12
11.5
12.5
12
10.5
9
7.9
6.7
5.7
4.8
3.6
2.8
2.3
1.8
1.4
1.1
1
0.9
0.7
0.4
-

1.6
2.4
1.25
100
87
6.1
0.1
1e-03
1e-03
0.06

3.0
0.21
-

Table 2: Spectral quantities used in the model (SI units).

From Wyszecki and Stiles [32], the relative spectral radiant power
SD () of a D-illuminant is given by a linear combination ofmean spectral radiant power S0 () and first two eigen vector functions S1 () and S2 () used in calculating daylight illuminants.
SD () = S0() + M1 S1() + M2 S2(). Scalar multiples M1
and M2 are functions of chromaticity values x and y and are given
by

M1

K

Solar position (s , s ) can be computed from the solar declination angle, latitude and longitude.

s,1:3515 , 1:7703x + 5:9114y ;0:0241 + 0:2562x , 0:7341y
0:0300 , 31:4424x + 30:0717y:
0:0241 + 0:2562x , 0:7341y

s

=
=

 , arcsin(sin l sin  , cos l cos  cos t );2, cos  sin tarctan();
cos l sin  , sin l cos  cos t

12

12

12

where s is solar angle from zenith in radians, s is solar azimuth
in radians, l is site latitude in radians,  is solar declination in radians and t is solar time in decimal hours. Solar angles from zenith
are between 0 and =2 and angles above =2 indicate sun below
horizon. Positive solar azimuthal angles represent direction west ofsouth.

A.6 Sun Position
Sun position is given by angle from zenith (s ) and azimuth angle(s ) and they depend on the time of the day, date, latitude and longitude (see Figure 4). Solar time can be calculated from the standard
time by using the formula, 80) ) , 0:129 sin( 2(J , 8) ) + 12(SM , L) ;t = ts + 0:170 sin( 4(J373355


A.7 Spectra
There are several spectral quantities used in the model: K for v = 4
used in the calculation of Mie scattering coefficient; S0 ; S1 ; S2
spectrums [32]; the sun’s spectral radiance in Wm,2 nm,1 sr,1 .
The latter was calculated from the spectral distribution of solar radiation incident at top of the atmosphere as adopted by NASA as astandard for use in engineering design [5]. These quantities can befound in Table 2. The spectral curves ko in m,1 , kwa in m,1 andkg used in the sunlight computation are also listed (Source: [13]).

where t is solar time in decimal hours, ts is standard time in decimal hours, SM is standard meridian for the time zone in radians, L
is site longitude in radians and J is Julian date (the day of the year
as an integer in the range 1 to 365).
The solar declination ( ) in radians is approximated by, 81) ):
 = 0:4093 sin( 2(J36899

References[20] M AX , N. L. Atmospheric illumination and shadows. In Computer Graphics (SIGGRAPH ’86 Proceedings) (Aug. 1986),
vol. 20, pp. 117–24.

[1] B LINN , J. F. Light reflection functions for simulation ofclouds and dusty surfaces. vol. 16, pp. 21–29.

[21] M C C ARTNEY, E. J. Optics of the Atmosphere. Wiley publication, 1976.

[2] B RUNGER , A. P., AND H OOPER , F. C. Anisotropic sky radiance model based on narrow field of view measurements ofshortwave radiance. Solar Energy (1993).

[22] M INNAERT, M. Light and Color in the Open Air. Dover,1954.

[3] B ULLRICH , K. Scattered radiation in the atmosphere. In
Advances in Geophysics, vol. 10. 1964.

[23] NAKAMAE , E., K ANEDA , K., O KAMOTO , T., AND
N ISHITA , T. A lighting model aiming at drive simulators.
In Computer Graphics (SIGGRAPH ’90 Proceedings) (Aug.
1990), F. Baskett, Ed., vol. 24, pp. 395–404.

[4] CIE-110-1994. Spatial distribution of daylight - luminancedistributions of various reference skies. Tech. rep., International Commission on Illumination, 1994.

[24] N IMEROFF , J., D ORSEY, J., AND RUSHMEIER , H. Implementation and analysis of an image-based global illumination
framework for animated environments. IEEE Transactions on
Visualization and Computer Graphics 2, 4 (Dec. 1996). ISSN
1077-2626.

[5] C OULSON , K. L. Solar and Terrestrial Radiation. Academic
Press, 1975.
[6] DA V INCI, L. The Notebooks of Leonardo da Vinci, vol. 1.
Dover, 1970.

[25] N ISHITA , T., D OBASHI, Y., K ANEDA , K., AND YA MASHITA , H. Display method of the sky color taking into
account multiple scattering. In Pacific Graphics ’96 (1996),
pp. 117–132.

[7] D OBASHI, Y., N ISHITA , T., K ANEDA , K., AND YA MASHITA , H. Fast display method of sky color using basis
functions. In Pacific Graphics ’95 (Aug. 1995).
[8] E BERT, D., M USGRAVE , K., P EACHEY, D., P ERLIN , K.,
AND W ORLEY. Texturing and Modeling: A Procedural Approach, second ed. Academic Press, 1998.

[26] N ISHITA , T., S IRAI, T., TADAMURA , K., AND NAKAMAE ,E. Display of the earth taking into account atmospheric scattering. In Computer Graphics (SIGGRAPH ’93 Proceedings)(Aug. 1993), J. T. Kajiya, Ed., vol. 27, pp. 175–182.

[9] E LTERMAN , L. Aerosol measurements in the troposphere andstratosphere. Applied Optics 5, 11 (November 1966), 1769–
1776.

[27] R. P EREZ , R. S EALS, J. M., AND I NEICHEN , P. An allweather model for sky luminance distribution. Solar Energy(1993).

[10] G OLDSTEIN , E. B. Sensation and Perception. Wadsworth,1980.

[28] R AYLEIGH , L. On the scattering of light by small particles.
Philosophical Magazine 41 (1871), 447–454.

[11] G RACE , A. Optimization Toolbox for use with MATLAB:
User’s Guide. The Math Works Inc., 1992.

[29] TAKAGI , A., TAKAOKA , H., O SHIMA , T., AND O GATA , Y.
Accurate rendering technique based on colorimetric conception. In Computer Graphics (SIGGRAPH ’90 Proceedings)(Aug. 1990), F. Baskett, Ed., vol. 24, pp. 263–272.

[12] I NEICHEN , P., M OLINEAUX , B., AND P EREZ , R. Sky luminance data validation: comparison of seven models with four
data banks. Solar Energy 52, 4 (1994), 337–346.

[30] WARD , G. J. The RADIANCE lighting simulation and rendering system. In Proceedings of SIGGRAPH ’94 (Orlando,
Florida, July 24–29, 1994) (July 1994), A. Glassner, Ed.,
Computer Graphics Proceedings, Annual Conference Series,
ACM SIGGRAPH, ACM Press, pp. 459–472. ISBN 0-89791667-0.

[13] I QBAL , M. An Introduction to Solar Radiation. Academic
Press, 1983.
[14] K AJIYA , J. T. The rendering equation. In Computer Graphics(SIGGRAPH ’86 Proceedings) (Aug. 1986), D. C. Evans and
R. J. Athay, Eds., vol. 20, pp. 143–150.

[31] WARD -L ARSON , G. Personal Communication, 1998.

[15] K ANEDA , K., O KAMOTO , T., NAKAME , E., AND N ISHITA ,T. Photorealistic image synthesis for outdoor scenery undervarious atmospheric conditions. The Visual Computer 7, 5 and6 (1991), 247–258.

[32] W YSZECKI , G., AND W.S.S TILES. Color Science. WileyInterscience publication, 1982.
[33] Y U , Y., AND M ALIK , J. Recovering photometric properties of architectural scenes from photographs. In SIGGRAPH
98 Conference Proceedings (July 1998), M. Cohen, Ed., Annual Conference Series, ACM SIGGRAPH, Addison Wesley,pp. 207–218. ISBN 0-89791-999-8.

[16] K ARAYEL , M., NAVVAB , M., N E ’ EMAN , E., ANDS ELKOWITZ , S. Zenith luminance and sky luminance distributions for daylighting calculations. Energy and Buildings
6, 3 (1984), 283–291.
[17] K AUFMAN , J. E., Ed. The Illumination Engineering SocietyLighting Handbook, Reference Volume. Waverly Press, Baltimore, MD, 1984.
[18] K LASSEN , R. V. Modeling the effect of the atmosphere on
light. ACM Transactions on Graphics 6, 3 (1987), 215–237.
[19] LYNCH , D. K., AND L IVINGSTON , W. Color and Light inNature. Cambridge University Press, 1995.

100