Stratospheric Ozone and Human Health Project
 |
RADIATIVE TRANSFER MODEL OVERVIEW: DATASET DESCRIPTION AND METHODOLOGY |
Darryl H. Charache, Ph.D
Consortium for International Earth Science Information Network
Socioeconomic Data and Applications Center
University Center, MI, USA
TABLE OF CONTENTS
- Model Overview
- Incorporation of Ozone into the Model
- Ozone Absorption Cross Sections
- Number Density Profiles: Nimbus-7
- Total Column Ozone: TOMS
- Incorporation of Aerosol into the Model
- Boundary Layer Aerosol Extinction: SAMSON
- Particle Size Distributions and Phase Functions
- Incorporation of Cloud into the Model
- International Satellite Cloud Climatology Project (ISCCP)
- Solar and Meteorological Surface Observation Network
- Cloud Particle Size Distributions and Phase Functions
- Incorporation of Topography into the Model
- Temperature and Pressure Profiles
- Calculation of Biologically Effective Dose
- Action Spectra
- Daily and Seasonal Dose Estimation
- Model Validation and Sensitivity Analysis
- Clear Sky Validation: Modeled vs. Measured Irradiance
- References
The radiative transfer model developed for this study simulates the effects of
absorption by ozone and multiple scattering by aerosol, molecules, and clouds
in the atmosphere. Calculations of the direct and diffuse components of
radiation were carried out with a wavelength resolution of two nanometers (nm)
between 297 and 399 nm, using the adding and doubling method (Hansen and Travis, 1974; Lacis and Hansen, 1974). The vertical resolution of the model for this study is
five hundred meters for altitudes less than ten kilometers; above this height
one kilometer layers were employed. Transmission and reflection coefficients
are computed at each layer using scattering and absorption extinction profiles
for clouds, aerosols, molecules, and ozone. The model also has the capability to incorporate multiple cloud layers of varying horizontal and vertical dimensions, a feature discussed in detail in Charache et al. (1994). For the purposes of this study, we will focus on single layer computations. Extraterrestrial
solar radiation at the top of the atmosphere is taken from the solar flux
atlas of Kurucz et al. (1984). The ground albedo used in all simulations was
set equal to 0.1. A description of the data sets incorporated into the RTM is
given in the following sections.
Extinction of UV radiation by ozone is found by the product of the
cross-sectional area (sigma) and the number density for each
atmospheric layer n(z),
(1)
Cross sectional values as a function of wavelength and temperature were
obtained from Paur and Bass (1984) for the wavelength range 297-339 nm, and
Cacciani et al. (1989) from 340-355 nm. Cross sections from both sources were
given in 1 Angstrom increments; these were averaged over a 10 Angstrom
bandwidth centered on each nanometer to find cross sections at 1 nm increments.
The temperature dependence of the cross sections were estimated by a linear
interpolation between given values of 228K and 295K. Values for 228K and 295K
are shown in Figure 1.
Seasonally averaged number density profiles at increments of 1 km in the
stratosphere and 0.5 km in the troposphere were interpolated from Umkehr level
measurements taken by the Nimbus-7 Solar Backscatter Ultraviolet (SBUV)
instrument (McPeters et al., 1984), and typical surface-level mixing ratios.
Umkehr data are given in monthly averaged 10deg. latitude bands for 1979. In
addition to these fixed model profiles, pollution episodes in the boundary
layer and troposphere can be simulated by defining ozone mixing ratios at any
level that will overwrite pre-defined profile values at these altitudes.
Figure 1 -- Wavelength and temperature dependence of ozone absorption cross section.
Seasonal profiles were then derived by normalizing the number density profiles
to the total Dobson unit values measured over the region by the Total Ozone
Mapping Spectrometer (TOMS) instrument. The TOMS version 6 dataset provides,
with few exceptions, continuous daily measurements of total ozone column depth
on a 1x1.25 degree grid from December 1978 through January 1992. An example of
this thirteen year variation of seasonally averaged ozone column depths is
displayed in Figure 2 for the southeastern Michigan region. Total column depth
of ozone shows peak values during the spring season (see Figure 3). Column
depth decreases during the summer and fall months, and begins to rise again in
winter. This seasonal cycle prevails throughout the thirteen year period.
Figure 2 -- Variation of total column ozone over SE Michigan, USA
Figure 3 -- Individual season variation of ozone over SE Michigan
Aerosol particles are an important component in determining the amount of UV
radiation reaching the surface. They arise from a variety of sources, both
natural and anthropogenic, and can vary significantly in scattering and
absorption properties depending on their composition and origin. In this study,
aerosol have been delineated into two categories based on their scattering
characteristics: boundary layer and troposphere-stratosphere. In all cases,
aerosol particles are assumed to be spherical so that Mie theory can be applied
in calculating scattering cross sections and phase functions. Properties of
each aerosol layer is described in the following sections.
Boundary layer aerosol extinction coefficients have been found to be inversely
proportional to the visibility range, and can be estimated by the empirical
Koschmieder equation
(2)
giving extinction coefficients in units of km-1. The boundary layer
in this model extends from the surface up to 2 km, with an aerosol scale height
of 3 km. Hourly visibility measurements for the selected locations were taken
from the Solar and Meteorological Station Observation Network (SAMSON) CD-ROM
database, available from the National Climatic Data Center (NCDC). Hourly
station observations are given at most stations for the years 1961-1990. Cloud
cover observations are also available from SAMSON, which are discussed in a
later section. Extinction profiles for upper tropospheric and stratospheric
aerosol as a function of wavelength were taken from the LOWTRAN 5 atmospheric
model (Kneizys et al., 1980). Several aerosol profile scenarios are available
depending on existing conditions: background stratospheric loading,
moderate, high, and extreme volcanic loading conditions. All data are given in
units of km-1.
Particle size distributions for the aerosol layers included in the model were
fitted to a modified gamma distribution as described by Deirmendjian (1969),
(3)
where r is the particle radius, n(r) is the number density of particles per
radius interval, and a, b, [[alpha]], and [[gamma]] are parameters used to fit
the observed distribution. Once the size distributions of the aerosols were
determined, Mie scattering theory was employed to calculate scattering phase
functions for the boundary layer and upper troposphere/stratosphere, using
indices of refraction given in Jursa (1985). Although model distributions
differ slightly for the troposphere and stratosphere, a sensitivity study of
surface radiation showed no dependence on the use of separate phase functions
for upper tropospheric and background stratospheric aerosol. Table 1 lists the
size parameters and index of refraction for each aerosol layer. For a
polydispersion, the phase function is taken as the average of all phase
functions for each particle radius weighted by the scattering cross section and
number density at that radius.
Layer a b [alpha] [gamma] Refractive Index
Boundary Layer 6.0e+7 18 1.5 1 1.572 - 1.0e-2i
Tropo/Strato 5.0e+7 19 1.55 1 1.53 - 8.00e-3i
Table 1 Gamma distribution size parameters and refractive index for boundary layer and upper troposphere/stratosphere aerosol layers.
ISCCP is a NASA-funded program which provides both daily and monthly-averaged
cloud and atmospheric parameters on a global scale. These data are compiled
from a suite of polar-orbiting and geostationary platforms and aggregated into
2.5x2.5 degree latitude/longitude gridded data product. The current database
runs from 1983-1990, but recent updates have been made to extend to later
years. These data are available on CD-ROM from the NASA Langley Research
Center's (NASA-LaRC) Distributed Active Archive Center (DAAC) free of charge.
Updates are also available from LaRC. Cloud data used for this study are
monthly averaged cloud optical depth values.
As mentioned in the previous section on aerosols, SAMSON contains hourly
station observations for 237 locations in the United States from 1961-1990.
These include total sky cover, opaque cloud cover, and ceiling height. Opaque
sky cover refers to non-transparent cloud layers, while total cloud cover may
include semi-transparent cirrus layers. Values for these parameters are given
in tenths, while ceiling height is given in meters above ground. These data
are used to infer the type of cloud present (cumuliform, stratiform, or
cirroform) and the amount of sky coverage in order to assign cloud optical depth and the clear sky/cloud cover fractional value used to weight irradiance values in partly cloudy conditions.
The Mie scattering code computes scattering cross sections and the resulting
phase function at one degree scattering angle increments from
0o-180o. Since the multiple scattering routine uses
Gauss-Legendre quadrature technique for integrating in spherical coordinates,
evaluation of the phase function at scattering angles defined by the zeros of
the Legendre polynomials is carried out by quadratic interpolation between the
values given by the Mie routine. This procedure gives good results when the
phase function is smooth, as in the case of a tropospheric or stratospheric
aerosol layer. It does not work well, however, when the phase function changes
rapidly with scattering angle, as in the case with cloud.
The sharply forward-scattered peak resulting from the large size parameter
involved with cloud droplets presents several obstacles to incorporation into
the multiple scattering routine. With such a narrow peak extending over 4 to 5
orders of magnitude, the phase function approaches a delta function, making
integration of this function difficult in the radiative transfer code. Even
with high order quadrature (>80), the peak of the Mie phase function was not
resolved sufficiently. To remedy this problem, the Henyey-Greenstein phase
function, an analytic expression based on the asymmetry parameter, was employed
(Henyey and Greenstein, 1941). The asymmetry parameter, denoted by
<cos[[zeta]]>, defined mathematically as the first moment of the Mie
phase function integrated over all values of the cosine of the scattering
angle, is given by
(6)
where u=cos[[zeta]] (e.g. Irvine, 1965). Once the asymmetry
parameter is known, the Henyey-Greenstein phase function is found from the
relation
(7)
where g=<cos[[zeta]]>. Using a Henyey-Greenstein expression
allows phase function values to be evaluated explicitly at any scattering angle
interval desired, as opposed to the angles defined by the zeros of the Legendre
polynomials using n-th order quadrature. This results in better
resolution of the forward scattering lobe.
From the phase functions calculated using the mode radius, linear weighted
mean, and surface weighted mean of the cloud particle size distribution,
asymmetry parameters of 0.874, 0.898, and 0.903, respectively, were found when
numerically integrating the asymmetry parameter expression. These values are in
good agreement for the given size parameter with previous studies and reviews
that employed Henyey-Greenstein phase functions (e.g. Hansen, 1968, Frederick
and Lubin, 1988, Madronich, 1993).
Although the Henyey-Greenstein expression can reproduce the forward scattering
lobe using asymmetry parameters of 0.95 or greater, it underestimates the
amount of scattering at intermediate angles in comparison to the Mie phase
function, and can give unrealistic results. Despite the fact that the forward
peak is not precisely reproduced using asymmetry parameters on the order of
0.8-0.9, results show good agreement with modeling studies of Frederick and
Lubin (1988) and Frederick and Snell (1990). Figure 5 shows a comparison of the
Mie phase function for a water droplet equal to the surface weighted mean
radius (4.0um) at a wavelength of 300 nm (dotted line) with the
Henyey-Greenstein phase function used in the model simulations having an
asymmtery parameter of 0.9 (solid line).
Figure 5 -- Mie phase function for equivalent droplet radius of 4 microns, and Henyey-Greenstein phase function with asymmetry parameter of 0.9.
With the Henyey-Greenstein function, Gaussian quadrature integration sampling
size in the zenith coordinate was reduced to 30 points in comparison to a
sampling of greater than 80 required for the Mie phase function. For the
azimuth coordinate, 12 points were sufficient in the phase function integration
scheme to resolve the necessary variation in the scattering pattern, in
comparison to the 50 point integration required to resolve azimuthal dependence
with the Mie function.
Effects of topography on the amount of radiation reaching the surface can be
significant; higher altitude regions will receive more direct beam radiation
than will lower altitude regions since the number of scattering events decrease
with the density of the atmosphere, resulting in less radiation being scattered
back to space (e.g. Blumthaler and Ambach, 1990). The vertical resolution for
scattering and absorption properties for the various atmospheric constituents
is 1 km for altitudes greater than 10 km, and has been modified to 0.5 km below
this height in order to consider effects of topography on dose amounts.
Altitude variations are crudely taken into account with this 500 m vertical
resolution using RAND's global elevation and depth dataset obtained from the
National Center for Atmospheric Research (NCAR); global resolution of this
dataset is on a 1ox1o posting, while the
horizontally-averaged vertical resolution of 10m is rounded to the nearest
one-half kilometer increment to match the model resolution.
Vertical profiles of temperature and pressure as a function of season and
latitude was taken from the LOWTRAN 5 model (Kneizys et al., 1980). LOWTRAN 5
defines six model atmospheres corresponding to the 1962 U.S. Standard
atmosphere and five supplementary models based on season and latitude:
Tropical, Midlatitude Summer, Midlatitude Winter, Subarctic Summer, and
Subarctic Winter. Vertical resolution for each profile is 1 km. The appropriate
data set was incorporated into the model based on the latitude of each
1ox1o gridpoint. Both temperature and pressure variables
are used in finding the Rayleigh phase function, while temperature as a
function of altitude was used in calculating the ozone absorption cross section.
Once the spectral irradiance at the surface is found, it must be
weighted with an action spectrum to obtain exposure values in the UV wavelength
range. Action spectra are a function of wavelength, since certain wavelengths
are more effective at producing a biological response than others; every
biological effect having a dose-response relation to radiation (e.g., DNA
mutation, erythema, plant damage) is characterized by its own particular action
spectrum. The erythemal action spectrum given by McKinlay and Diffey (1987) is
shown for comparison along with the non-melanoma skin cancer (NMSC) spectrum of
de Gruijl and van der Leun (1993) in Figure 6. In most applications, a
normalized action spectrum is used in determining the potential biologically
effective dose of radiation for a given effect. A normalized action spectrum is
determined relative to the biological response produced at a particular
wavelength, resulting in action spectrum values ranging between zero (no
biological effect) and one (biological effect at the reference wavelength).
Erythemal effectiveness is normalized to the response at 290 nm, and NMSC is
normalized relative to about 296 nm. Both curves show similar relative
responses up to approximately 325 nm, beyond this point erythema shows a
greater response to the longer UV-A wavelengths. Note, however, that a direct
comparison of received dose amounts for any induced effect cannot be made on
the basis of the action spectrum alone; this only gives an indication as to the
relative effectiveness of particular wavelengths, not the actual dose amount
required to induce a certain biological response. Two UV-induced effects may
have similar looking action spectra indicating a similar wavelength dependence,
but one may require a much larger dose amount to produce the particular
response. In the case of non-melanoma and erythema, for example, erythema has a
reasonably well-quantified dose-response threshold, while skin cancer does not;
there is no quantifiable amount of UV that is known to induce the formation of
cancerous cells.
The potential biologically effective dose of UV at the surface for a given
instant in time is found by a convolution of the UV spectrum (UV) and the
action spectrum (A):
(9)
Figure 6 -- Non-melanoma action spectrum of de Gruijl and van der Leun
(1993), and the standardized erythema action spectrum of McKinlay and Diffey
(1987).
Potential dose is the total amount of biologically effective radiation received
on a horizontal surface assuming no confounding environmental or behavioral
factors modifying this quantity; it is the total amount available to begin
induction of a particular biological effect. It is only a potential quantity,
however, since factors like clothing habits and geometrical orientation to the
ground and sun can modify the actual amount of effective radiation reaching the
cells. Incorporation of these factors would require a much more detailed study
on a local scale, and are not considered in the scope of this study.
The radiation code calculates instantaneous spectral irradiance in units of
W/cm2/nm for a given solar zenith angle. Calculation of zenith angle
is dependent upon latitude, longitude, day of the year, and time of day, and is
found through the use of spherical trigonometry applied to the celestial sphere
coordinate system. The first step in evaluating the solar zenith angle is
calculating the declination angle (d) and the equation of time (EQT). These two
values are approximated by a Fourier series expansion based on a parameter
derived from the day of the year (dn), denoted by the Greek letter theta and defined
as
(10)
where dn runs from 0 (1. January) to 364 (31. December). The fourier series
expansion in [[theta]] for d and EQT is given by Spencer (1971):
(11)
(12)
Values returned are in radians. The equation of time, once converted to time
units, is then added to the local clock time and the longitudinal correction
factor to come up with the local apparent time. The longitude correction is
based on the rotation rate of the Earth; it takes approximately 4 minutes for
the sun to appear to move through one degree of longitude. For example, in the
eastern time zone of the United States, the sun is highest in the sky at
75deg.W. In southeast Michigan, located at approximately 84deg.W, it takes
about 36 minutes past local noon for the sun to reach its highest point in the
sky. This factor, included with the clock time and equation of time, gives the
local apparent time used in astronomical calculations:
Local Apparent Time = Clock Time + Longitude Correction + EQT (13)
The local apparent time is then used to find the hour angle of the sun in the
sky. The hour angle is defined as the angle measured clockwise from the south
meridian in the sky. This value changes continuously due to the rotation of the
Earth and is found first by subtracting local apparent time from local noon (12
hours). This number, which can be either positive or negative, is then
multiplied by the number of degrees of longitude that the Earth rotates in one
hour (15deg./hour) to give the hour angle in units of degrees.
Once the hour angle and declination angle are known, the zenith angle at
latitude [[phi]] can be derived using spherical trigonometry relations; the
final result is given by the equation
(14)
or equivalently
(15)
This procedure is included in the multiple scattering code to compute the sun's
position at any location on Earth; once the zenith angle is known, the spectral
irradiance can be found at any time of day.
The process involved in estimating daily UV totals begins with tracking the
progression of the sun's position over the course of the day. The scattering
routine first calculates spectral irradiance at 5deg. solar zenith angle
intervals from 0deg. to 80deg., and integrates over the UV wavelength region.
Then the zenith angle is calculated at 1 hour intervals over the course of the
given day, and the total UV flux for that instant in time is estimated by
interpolating between values calculated at the 5deg. spacing. Using this method
allows determination of flux values at any time and place without prior
knowledge of the zenith angle at a particular time, since exact zenith angle
values are computed internal to the routine. This saves both CPU, and
especially interactive operating time, as a result. Errors incurred by using
this technique as compared to explicit evaluation of fluxes at the zenith angle
for each time step are less than one percent. Daily totals are then found by
integrating over these 1 hour time steps:
(16)
8.1 Clear Sky Validation: Modeled vs. Measured Irradiance
The accuracy of the model was tested by comparing model output to clear sky
measurements taken with a Biospherical Instruments GUV-511 device. This
instrument was installed at the Department of Atmospheric, Oceanic, and Space
Sciences as part of a nationwide UV monitoring network funded by Bausch &
Lomb. The GUV-511 is a temperature-controlled unit providing spectral
measurements in units of uW cm-1 nm-1 at 305 nm, 320 nm,
340 nm, and 380 nm with FWHM bandwidths ranging from 8 to 10 nm. Erythemal dose
and Photosynthetically Active Radiation (PAR) measurements are also provided. A
Compaq ProLinea 4/25s PC collects the data using a data acquisition system
provided by Global Environmental Engineering Consultants, Inc. that records
one-minute averages in daily output files. Clear sky cases only will be
considered in this validation study due to lack of data on cloud optical depths
over the local area. An inherent difficulty in validating the cloud model is
the high spatial and temporal variability of cloud; partly cloudy skies are
continually changing over the course of minutes, and reflection from edges of
cloud along with diffuse contributions from clear and cloudy areas not directly
in the line of sight make this a difficult, if not impossible problem to
consider. Even in seemingly invariant overcast conditions, cloud parameters
remain fluid as evidenced in data collected from August 2, 1994 (Figure 7).
Overcast conditions persisted throughout the day, but considerable variation in
the stratus layer is evident in the UV irradiance measured at the surface.
Clear sky conditions were observed on September 11, 1994 in the southeast
Michigan region. Total column ozone in the 1ox1.25o grid
cell over the Ann Arbor coordinates of 42.15oN, 83.45oW
for this day was taken from the on-line database of TOMS Meteor-3 observations.
NASA made these data available in near real-time with the understanding that
these data are preliminary figures and should not be used for detailed research
purposes until validated. With this caveat, the daily column ozone value over
this location was 305 Dobson units (DU). This value was coded into the ozone
extinction program and the model ozone profile for mid-latitude autumn was
normalized to this value. Visibility was high on this day; a correspondingly
low surface layer aerosol extinction coefficient of 0.2 km-1 was
incorporated into the boundary layer aerosol model.
UV irradiance was calculated every half hour from 8AM to 7PM Eastern Daylight
Time (EDT) based on variation in solar zenith angle calculated by the procedure
outlined in section 7.2. Ozone and aerosol levlels were assumed constant over
the day for simplicity. A comparison of calculations and isntrument output at
305 nm, 320 nm, and 340 nm is shown in Figure 3.8. Model estimates showed
variation of less than 8% at 305 nm and 320 nm, and less than 4% at 340 nm
relative to measurements. With uncertainties in boundary layer aerosol and
ozone in addition to the non-validated total column ozone value, these
estimates give confidence that this model provides an accurate representation
of the UV budget at the surface of the Earth.
Figure 7 -- UV irradiance under overcast conditions on 8/2/94.
Figure 8 -- Model comparison to clear sky measurements on 9/11/94.
Blumthaler, M. and W. Ambach, Indication of increasing solar ultraviolet-B
radiation flux in alpine regions. Science, Vol. 248, pp. 206-208,
1990.
Cacciani, M., A. Di Sarra, G. Fiocco, and A. Amoruso, Absolute determination of
the cross sections of ozone in the wavelength region 339-355 nm at temperatures
220-293 K. Journal of Geophysical Research, Vol. 94, pp. 8485-8490,
1989.
Charache, D. H., V. J. Abreu, W. R. Kuhn, and W. R. Skinner, Incorporation of
multiple cloud layers for ultraviolet radiation modeling studies. Journal of
Geophysical Research, Vol. 99, pp. 23,031-23,040, 1994.
Deirmendjian, D., Electromagnetic Scattering on Spherical
Polydispersions, American Elsevier, New York, 1969.
Frederick, J. E., and D. Lubin, The budget of biologically active ultraviolet
radiation in the Earth-atmosphere system. Journal of Geophysical
Research, Vol. 93, pp. 3825-3832, 1988.
Frederick, J. E., and H. E. Snell, Tropospheric influence on solar ultraviolet
radiation: the role of clouds. Journal of Climate, Vol. 3, pp. 373-381,
1990.
Hansen, J. E., Exact and approximate solutions for multiple scattering by
cloudy and hazy planetary atmospheres. Journal of the Atmospheric
Sciences, Vol. 26, pp. 478-487, 1968.
Hansen, J. E., and L. D. Travis, Light scattering in planetary atmospheres.
Space Science Reviews, Vol. 16, pp. 527-610, 1974.
Henyey, L. C., and J. L. Greenstein, Diffuse radiation in the galaxy,
Astrophysical Journal, Vol. 93, pp. 70-83, 1941.
Irvine, W. M., Light scattering by spherical particles: radiation pressure,
asymmetry factor, and extinction cross section, Journal of the Optical
Society of America, Vol. 55, pp. 16-21, 1965.
Jursa, A. S., Handbook of Geophysics and the Space Environment, Air Force
Physics Laboratory, Air Force Systems Command, USAF, 1985.
Kneizys, F. X., E. P. Shettle, W. O. Gallery, J. H. Chetwynd, Jr., L. W. Abreu,
J. E. A. Selby, R. W. Fenn, and R. A. McClatchey, Atmospheric
transmittance/radiance: Computer Code LOWTRAN 5, no. 687, AFGL-TR-80-0067, Air
Force Geophysical Laboratory, Hanscom AFB, Massachusetts, 1980.
Kurucz, R. L., I. Furenlid, J. Brault, and L. Testerman, Solar Flux Atlas from
296 to 1300 nm, National Solar Observatory Atlas No. 1, June 1984.
Lacis A. A., and J. E. Hansen, A parameterization for the absorption of solar
radiation in the Earth's atmosphere, Journal of the Atmospheric
Sciences, Vol. 31, pp. 118-133, 1974.
Madronich, S., Chapter 2, UV radiation in the natural and perturbed atmosphere,
in UV-B Radiation and Ozone Depletion: Effects on Humans, Animals, Plants,
Microorganisms, and Materials, M. Tevini, ed., Lewis Publishers, Boca Raton,
Florida, 1993.
McKinlay, A. F., and B. L. Diffey, A Reference action spectrum for ultra-violet
induced erythema in human skin, in Human Exposure to Ultraviolet Radiation,
Risks and Regulations, W. F. Passchier and B. F. M. Bosnjakovic, ed.,
International Congress Series 744, pp. 83-87, 1987.
McPeters, R. D., D. F. Heath, and P. K. Bhartia, Average ozone profiles for
1979 from the Nimbus-7 SBUV instrument, Journal of Geophysical Research,
Vol. 89, pp. 5199-5214, 1984.
Paur, R. J., and A. M. Bass, The ultraviolet cross-sections of ozone: II.
Results and temperature dependence, in Atmospheric Ozone, edited by C.S.
Zerefos and A. Ghazi, pp. 611-616, D. Reidel, Dordrecht, 1984.
Spencer, J. W., Fourier series representation of the position of the sun.
Search, Vol. 2, p. 272, 1971.
