For the user interested in the details of the methods used, here are the actual algorithms used in selecting the ``nearest'' sky, and in the sky-brightness model:
In the ``nearest'' method, a distance estimator is computed that includes separations both in time and on the sky. The estimator is
where the t's are times in decimal days,
the AM's are air masses, and the AZ's are
azimuths (in radians) of the two observations being compared.
(The azimuth difference is always taken to be less than 
radians.)
S crudely takes account of the greater variation in sky brightness with
position near the horizon.
At moderate air masses, it makes a separation of a minute in time about
equivalent to a degree on the sky.
This estimator is computed for the two sky samples closest in time to each star observation, one before the star and the other after it. Only observations made with the same filter (and, if diaphragm size is available, with the same diaphragm) are used. The sky observation that gives the smaller value of S is the one used in the ``nearest'' sample method. Obviously, if the star is the first or last of the night, there may be no sky sample on one side; then the one on the other side (in time) is used.
The angular part of S is necessary to prevent problems when groups of variable and comparison stars are observed. It can happen that the sky is observed only after two stars in a group have been observed, if only a single sky position is used for the whole group. If sky was measured after the last star before the group, that previous sky may be closer to the first star of the group, in time alone, than the appropriate sky within the group. Then a purely time-based criterion would assign the (distant) previous star's sky to the first star of the group, instead of the correct (following) sky. Similar effects can occur at the end of a group, of course.
While this problem can be prevented by careful observers, not all observers are sufficiently careful to avoid it. The crude separation used here is adequate to resolve the problem without going into lengthy calculations. Note that both the airmass dependence of sky brightness and the possible presence of local sky-brightness sources around the horizon (e.g., nearby cities) make the horizon system preferable to equatorial coordinates for this purpose.
This brings us to the sky model. The general approach is to represent the sky brightness as the sum of two terms, a general one due to scattered starlight, zodiacal light, and airglow; and an additional moonlight term that applies only when the Moon is above the horizon.
The airglow and scattered starlight are proportional to the airmass, to a first approximation. Actually, extinction reduces the brightness of the sky light near the horizon. However, the full stellar extinction appears only in the airglow and zodiacal components, not in the scattered light. All three components are of comparable magnitude in the visible part of the spectrum.
While Garstang [6,7,8] has made models for the night-sky brightness, these require much detailed information that is not usually available to the photometric observer. Garstang's models also were intended to produce absolute sky brightnesses, while data at this preliminary stage of reduction do not have absolute calibrations. Finally, they do not include moonlight. Therefore, a simpler, parametric model is used.
Some guidance regarding the scattered starlight can be obtained from the multiple-scattering results given by van de Hulst [24]. In photometric conditions, the aerosol scattering is at most a few per cent, and the total optical depth of the atmosphere is less than unity. In the visible, the extinction is dominated by Rayleigh scattering, which is not far from isotropic, and nearly conservative. Therefore, we are interested in cases with moderate to small optical depth, and conservative, nearly isotropic scattering. Because the light sources (airglow, stars, and zodiacal light) are widely distributed over the sky, we expect small variations with azimuth, and can use the values in van de Hulst's Table 12 to see that azumuthally-averaged scattering has the following properties:
1. For air masses such that the total optical depth along the line of sight is less than 1, the brightness is very nearly proportional to air mass, regardless of the altitude of the illuminating source.
2.
For vertical optical depths 
less than about 1.5, the sky brightness reaches a
maximum at an air mass on the order of 
, and then declines to a fixed
limiting value as 
(remember that for the plane-parallel
model, the airmass does go to infinity at the horizon).
The decrease in the scattered light at the horizon is also to be expected in the direct airglow and zodiacal components attenuated by extinction; so the same general behavior is expected for all components. The simplest function that has these properties is
where M is the airmass; the limiting brightness at the horizon is just 
.
Actually, a substantially better fit can be obtained to the values in van de
Hulst's Table 12 by including a linear term in the denominator; so
the approximation
is used to represent the airglow and scattered light.
There is a problem in fitting such a function to sky data that cover a limited range of airmass. Except for optical depths approaching unity (i.e., near-UV bands), the maximum in the function lies well beyond the range of airmasses usually covered by photometric observations. That means that the available data usually do not sample the large values of M at which the squared terms become important. Thus, one can usually choose rather arbitrary values of these terms, and just fit the well-determined linear terms. It turns out that choosing b = 0 is often satisfactory. If the data bend over enough, d can be determined; otherwise, it defaults to 0.01 times the square of the largest airmass in the data.
An example of data that extend to large enough airmass to determine all four
parameters is the dark-sky brightness table published by Walker [25],
which were used by Garstang to check his models.
These data extend to 
.
The model above fits them about as well as do Garstang's models; typical errors
are a few per cent.
Various subsets, with the largest-Z data omitted, give similarly good fits.
This indicates that the model is adequate for our purposes here.
In principle, one could add separate terms for zodiacal and diffuse galactic light that depend on the appropriate latitudes; but this seems an excessive complication, as these components vary with wavelength and longitude as well. We have also neglected the ground albedo. Unless the ground is covered with snow, this is a minor effect except near the horizon. Furthermore, the airglow can vary by a factor of 2 during the night; so we cannot expect a very tight fit with any simple formula.
The moonlight term is more complicated. In principle, it consists of single scattering, which in turn depends on the size and height distributions of the aerosols, as well as Rayleigh scattering from the molecular atmosphere; and additional terms due to multiple scattering. The radiative-transfer problem is complicated further by the large polarization of the Rayleigh-scattered component, which can approach 100%. Rather than try to model all these effects in detail, we adopt a simple parametric form that offers fair agreement with observation, but does not have too many free parameters to handle effectively.
First of all, van de Hulst [24] points out that the brightness of the solar aureole varies nearly inversely with elongation from the Sun. We assume the lunar aureole has the same property. And, for the small optical depths we usually encounter, and the moderate to small airmasses at which we observe, we can simply assume the brightness of the lunar aureole is nearly proportional to airmass.
Second, interchanging the directions of illumination and observation would
give geometries related by the reciprocity theorem if the ground were black.
For typical ground albedoes, we can still assume approximate reciprocity.
We can also assume 
, where M is the airmass, accurate to a few
per cent for actual photometric data, in calculating elongations from the Moon.
The adopted form is
where M is the airmass in the direction of observation, E is the angular elongation from the Moon, S is the sum of observed and lunar airmasses, and P is their product. (Note that the lower-case parameters here are different from the ones in the dark-sky model.)
The factor in parentheses mainly represents the single-scattering phase
function, and should be nearly constant in good photometric conditions.
It is plotted as a ``normalized'' sky brightness.
Its parameter a is a measure of the lunar aureole strength; if 
is large, you probably have non-photometric conditions.
The factor in brackets handles the reciprocity effects.
The 
term produces the correct asymptotic behavior for a
homogeneous atmosphere; however, it cannot represent ``Umkehr'' effects at
wavelengths where ozone absorbs strongly.
Both factors have the symmetry required by the reciprocity theorem. This condition may be violated by ground-albedo effects, and by photometers that have a large instrumental polarization.
Unfortunately, when the telescope is pointed so that the Moon can shine on the objective, the apparent aureole is nearly always dominated by light scattered from the telescope optics, not from the atmosphere. Even if the mirror has been freshly aluminized, the scattered light may not be negligible, because of surface scattering. This scattering has different angular and wavelength dependences from those of the true sky brightness, and so is not represented by the sky model. This means we should avoid observing so near the Moon that it can shine directly on the primary mirror.
However, we cannot avoid having the star shine on the mirror; so we must
include a term in the ``sky'' model that is
proportional to the brightness of the measured star, to allow for the wings of
the star image that we include in our ``sky'' measurements.
Even if this fraction is so small (say, below 
) as to have no effect
on the sky-subtracted stellar data, it can be important in the sky model,
if fairly bright stars are observed.
The fraction of starlight measured depends on the distance from the star to the
sky area chosen; as mentioned before,
one must keep this distance fixed for all stars when observing.