# Self-Consistent Atmosphere Representation and Interaction in Photon Monte
Carlo Simulations

J. R. Peterson, G. Sembroski, A. Dutta, C. Remacaldo  Department of Physics
and Astronomy, Purdue University, West Lafayette, IN 47907, USA John R.
Peterson<EMAIL_ADDRESS>
###### Abstract

We present a self-consistent representation of the atmosphere and implement
the interactions of light with the atmosphere using a photon Monte Carlo
approach. We compile global climate distributions based on historical data,
self-consistent vertical profiles of thermodynamic quantities, spatial models
of cloud variation and cover, and global distributions of four kinds of
aerosols. We then implement refraction, Rayleigh scattering, molecular
interactions, Tyndall-Mie scattering to all photons emitted from astronomical
sources and various background components using physics first principles. This
results in emergent image properties that include: differential astrometry and
elliptical point spread functions predicted completely to the horizon,
arcminute-scale spatial-dependent photometry variations at 20 mmag for short
exposures, excess background spatial variations at 0.2% due the atmosphere,
and a point spread function wing due to water droplets. We use a common
atmophere representation framework to self-consistently model all phenomena by
simulating individual photons. We reproduce the well-known correlations in
image characteristics: correlations in altitude with absolute photometry
(overall transmission) and relative photometry (spectrally-dependent
transmission), anti-correlations of altitude with differential astrometry
(non-ideal astrometric patterns) and background levels, and an anti-
correlation in absolute photometry with cloud depth. However, we also find
further subtle correlations including an anti-correlation of temperature with
background and differential astrometry, a correlation of temperature with
absolute and relative photometry, an anti-correlation of absolute photometry
with humidity, a correlation of humidity with Lunar background, a significant
correlation of PSF wing with cloud depth, an anti-correlation of background
with cloud depth, and a correlation of lunar background with cloud depth.

telescopes– instrumentation: detectors– atmosphere

##  1 Introduction

Simulations have become increasingly indispensible in modern astrophysics. In
particular, a wide range of simulation approaches have been pursued in X-ray
astronomy (Peterson, Jernigan, & Kahn 2004; Peterson, Marshall, & Andersson
2007; Andersson, Peterson, & Madejski 2007; Davis et al. 2012) and optical
astronomy (Ackermann et al. 2012; Lane, Glindemann & Dainty 1992, Ellerbroek
2002; Le Louarn 2002; Britton 2004; Jolissaint 2010; Bertin 2009; Dobke et al.
2010; Peterson et al. 2015, Rowe et al. 2015). In our previous work, we have
created an ab initio physics simulator (PhoSim) that follows individual
photons through the atmosphere telescope and camera (Peterson et al. 2015,
Peterson et al. 2019, Burke et al. 2019a, Peterson et al. 2019). The detailed
microphysics of photon and electron interactions are followed to predict the
final pixel (or removal) of individual photons. Recent work has extended the
first principles physics to include the deformation physics of optics
(Peterson et al. 2019) and the electrostatic physics of sensors (Peterson et
al. 2020).

A number of publications have demonstrated a diverse set of applications for
an ab initio simulator in astronomy using PhoSim (https://www.phosim.org).
PhoSim has been used to: 1) to test the design of Rubin/LSST (Angeli et al.
2016, Thomas et al. 2016 Xin et al. 2015) and JWST (Burke et al. 2019a) 2) to
plan for future observations (e.g. Merlin et al. 2022, Bretonniere 2022,
Sanchez et al. 2020, Thomas & Kahn 2018 Bard et al. 2016, Mandelbaum et al.
2014, Bard et al. 2013 Chang et al. 2013a, Chang et al. 2013b, Chang et al.
2012), 3) for advanced image processing algorithm development (e.g. Nie et al.
2021a, Nie et al. 2021b, Carlsten et al. 2018, Li et al. 2016 Meyers & Burchat
2015), 4) to understand physical effects (Walter 2015, Beamer et al. 2015),
and 5) for advanced AI/machine learning development by simulating training
sets (Burke et al. 2019b). In addition to published studies, hundreds of users
have used PhoSim for informal studies. Most commonly the applications include:
1) simulations to plan for future observations, 2) testing analysis algorithms
on a simulation with perfect information, and 3) developing the designs of
future telescopes. Several terabytes of images have been produced for a
variety of data challenges for multiple observatories. These applications have
also led to the community collectively validate the physics implementation in
PhoSim.

It is increasingly important to consider the full self-consistent
representation of the objects while simulating the photon (and electron)
propagation and interactions. This is important since emergent effects that
are imprinted on images, including the point-spread-function (PSF) size and
shape and the photometric and astronometric patterns, are directly affected in
correlated ways. This is difficult since it involves simulating the ab initio
physics within a common framework and does not allow individual effects to be
treated in isolation. Unfortunately, that leads to a much more complicated
implentation. Then, it is a difficult to build more detail without the
simulator both becoming opaque to the user that contains hidden assumptions
and degenerating into one of many possible models for a particular aspect of
the simulation. However, in previous work with PhoSim, we have left the
modeling aspect of the simulation (i.e. the instrument design details) exposed
as a tiered set of user inputs, and built the simulator entirely around the
physical response of the photons and subsequent electrons. In that way, the
implementation of the latter would be unambiguous and follow ab initio
physics.

We apply this concept in this work with the simulation of a self-consistent
representation of the atmosphere. Then, the photon interactions are applied
unambigously on that particular representation of the atmosphere. However, the
simulator is built so that representation of the atmosphere can also be
modified or simplified by the user. The observational features are not
affected by the representation of the atmosphere because they only depend on
the photon interactions, but it is critical that the photon interactions are
simulated in a self-consistent basis with the same atmosphere representation
in order to reproduce correlations of observables. Therefore, we first compile
a global representation of the atmosphere based on a number of measurements in
the literature in §2. The subsequent photon interactions with the self-
consistent representation are described in §3. The resulting observational
effects are described in §4, and the correlated patterns are outlined in §5.
Conclusions and future work is described in §6. We limit the scope to the non-
turbulent aspect of the atmosphere as the turbulent atmosphere can be handled
as in Peterson et al. 2015 with one exception.

This approach improves over standard methologies in three key ways. The first
is that we compile a complete global representation of the atmosphere. We use
theoretical and empirical descriptions of the atmosphere to describe
thermodynamic quantities for every location on the globe. As discussed above,
although the reference implementation may not be perfect and it depends on the
current fidelity of the atmosphere modelling in the literature, the response
of light to the represented atmosphere follows exact first principles physics.
Having a complete global description is increasingly important in the PhoSim
package since dozens of ground-based telescopes are implemented, and it lowers
the burden to the user to add additional telescopes. This leads to the ability
to compare image of different observatories across the globe in a straight-
forward manner.

The second improvement is that emergent phenomena are uniquely predicted with
no additional parameters. In particular, both the background emission pattern
and the point-spread-function wing due to the atmosphere are self-consistently
predicted along with layer transmission absorption functions and differential
chromatic refraction. In a photon by photon approach through segments of the
atmosphere, the emergent observables are predicted. Thus, observables such as
the differential refraction pattern, the background level, point-spread-
function wing, and transmission properties are uniquely predicted from the
choice of atmosphere representation.

Finally, the third improvement is by using a physics-based approach some of
the common approximations and simplification are not used in favor of an exact
approach. For instance, given that refraction affects the path length of light
through the atmosphere and that the atmospheric composition varies with
altitude, means that using a single airmass quantity is not appropriate.
Similarly, in some cases the transmission of light calculations assume that
light is lost when it is actually scattered and still present in the image.
Additionally, the curvature of the Earth affects transmission and refraction
patterns at large zenith angles. In most cases, these complications are
actually solved in a straight-forward manner with a photon Monte Carlo
approach, despite being difficult to implement with analytic approaches.

We expect that the detailed implementation in this work may be insufficient to
predict exact long-term atmosphere patterns for a given site. That application
is best studied with specific observations. However, there are many
applications to predict realistic astronomical images that could be obtained
at a given site and this work combines ab initio physics and careful
representations cover the globe. In particular, three common uses are: 1)
detailed simulation of a telescope while it is being designed, built, and
commissioned to understand scientific performance and explain complex
phenomena, 2) planning of future observations for a given telescope before the
data is obtained, 3) simulation of realistic images with perfect input
knowledge to test any complex image processing pipelines with associated
algorithms, and 4) simulation of realistic training sets for machine
learning/artificial intelligence applications. We also estimate in the
conclusion future and current validation.

##  2 Input Self-Consistent Atmosphere Representation

In this section, we develop the framework for the input to PhoSim of a self-
consistent atmosphere representation. We will need to describe the relative
quantity of different constituents at different locations both as a function
of altitude and across the globe. This includes molecules, water droplets, and
aerosols, including the size distributions of water droplets and aerosols.
Similarly, we will self-consistently represent thermodynamic quantities, such
as temperature and pressure. Various parts of this will be used for the photon
interactions used in §3. The most important feature is that we use the same
atmosphere representation for all of the photon interactions. We also allow
any aspect of this representation to be modified by the user. For example, in
§4 we consider modifying the ground temperature but then letting the
temperaure as a function of altitude be determined from the input data below.

The complete description of the atmosphere is a very complex problem where
there is no current theoretical solution. Therefore, we will rely on a
combination of empirical historical data augmented with theoretical profiles
where there is insufficient detail. The accuracy required is not simple to
define, since we will use the thermodynamic quantities and molecular and
particle abundances to predict a number of interactions including
transmission, refraction, and scattering in §3. However, we will use as
comprehensive representation as possible that will reasonably approximate the
range of atmospheric conditions.

###  2.1 Climate Historical Data

Figure 1:  Temperature maps across the globe for the month of January. The top
left shows the mean temperature, the top right is the standard deviation for
that month, the bottom left shows the mean diurnal range, and the bottom right
shows the standard deviation of the diurnal range. We compute the statistics
using the data from Kalnay et al. (1996) as described in the text. The scales
go from purple to red from 257 to 296 K, 0.50 to 5.7 K, 0.44 to 5.7 K, and
0.44 to 2.7 K for the plots, respectively. Figure 2:  Pressure maps across the
globe for the month of January. The top left shows the mean pressure, the top
right is the standard deviation for that month, the bottom left shows the mean
diurnal range, and the bottom right shows the standard deviation of the
diurnal range. We compute the statistics using the data from Kalnay et al.
(1996) as described in the text. The scales go from purple to red from 87000
to 105963 Pa, 270 to 1135 Pa, 269 to 622 Pa, and 130 to 440 K for the plots,
respectively. Figure 3:  Relative humidity maps across the globe for the month
of January. The top left shows the mean humidity, the top right is the
standard deviation for that month, the bottom left shows the mean diurnal
range, and the bottom right shows the standard deviation of the diurnal range.
We compute the statistics using the data from Kalnay et al. (1996) as
described in the text. The scales go from purple to red from 66 to 94 %, 3.0
to 15 %, 5.7 to 25 %, and 2.7 to 16 % for the plots, respectively. Figure 4:
Wind patterns across the globe at high altitude consisting of a pressure of 10
mb in January. The top left shows the mean longitudinal velocity component,
the top right shows the variance of the logitudinal velocity component, the
bottom left shows the mean latitudinal velocity component, and the bottom
right shows the variance of the latitudinal velocity component. We compute the
statistics using the data from Kalnay et al. (1996) as described in the text.
The scales go from purple to red from -3.9 to 4.0 m/s, 0.85 ot 2.5 m/s, -2.6
to 2.2 m/s, and 0.77 to 2.1 m/s for the plots, respectively. Figure 5:  Wind
patterns across the globe at high altitude consisting of a pressure of 1000 mb
in January. The top left shows the mean longitudinal velocity component, the
top right shows the variance of the logitudinal velocity component, the bottom
left shows the mean latitudinal velocity component, and the bottom right shows
the variance of the latitudinal velocity component. We compute the statistics
using the data from Kalnay et al. (1996) as described in the text. The scales
go from purple to red from -21 to 17 m/s, 2.5 to 15 m/s, -8.8 to 8.7 m/s, and
-0.86 to 8.3 m/s for the plots, respectively.

To determine the distribution of thermodynamic quantities at the surface, we
use the National Centers for Environmental Protection/National Center for
Atmospheric Research (NCEP/NCAR) database which provides extensive information
across the globe from 1948 to the present day (Kalnay et al. 1996). For the
surface temperature, surface pressure, and surface relative humidity, we use
the entire database on 2.5x2.5 degree grid that makes a series of measurements
4 times per day. For a given latitude/longitude position, we bin the
measurements into monthly groups including all years, so we can capture the
seasonal variation. Then, within each monthly dataset we calculate four
quantities: the mean value, the standard deviation, the mean diurnal variation
from the 4 daily points, and the standard deviation of the diurnal variation.
The diurnal variation mean and standard deviation are multiplied by a
correction factor of due to the underestimate of the variation. If the daily
variation is only sampled 4 times and is approximately sinusoidal then the
variation is underestimated by a factor of 1.11. We repeat this calculation
for the surface temperature, pressure, and relative humidity. This then allows
a representation of the seasonal and dirunal variation for every location on
the Earth at a grid spacing of 2.5 degrees. The results are shown in Figures
1, 2, and 3.

We also capture the wind vector patterns using another related data set from
Kalnay et al. (1996). We use the vector mean and standard deviation can be
captured on a monthly basis since 1948. However, here we can obtain different
wind vectors at a series of different pressure locations which allows us to
construct the velocity vector at different altitudes. There are 17 altitudes
corresponding to different pressure levels and therefore different altitudes.
We then compute the mean and standard deviation of each component of the
velocity vector for each month of the year on each location on the Earth at a
2.5 degree grid. The global wind patterns for two altitudes are shown in
Figure 4 and 5.

###  2.2 Vertical Profiles of Thermodynamic Quantities

Figure 6:  The temperature profile as a function of altitude. The colored
lines indicate the 5 different reference models (equator, mid-latitude
winter/summer, sub-arctic from Anderson et al. (1986)), and the profiles
diverge at high temperature due to different asymptotic exosphere temperatures
from Jaccia (1971). This interpolates and combines the results from Anderson
et al. (1986), Champion, Cole, & Kantor (1976), and Jaccia (1971). Figure 7:
The pressure profile as a function of altitude. The colored lines indicate the
5 different reference models (equator, mid-latitude winter/summer, sub-arctic
from Anderson et al. (1986), and the profiles diverge at high temperature due
to different asymptotic exosphere temperatures from Jaccia (1971). This
interpolates and combines the results from Anderson et al. (1986), Champion,
Cole, & Kantor (1976), and Jaccia (1971). Figure 8:  The relative abundance
profiles as a function of altitude. The most important molecular constituents
are labelled. The colored lines indicate the molecules that are most variable
and are therefore different in the 5 reference models (equator, mid-latitude
winter/summer, sub-arctic from Anderson et al. (1986)). The abundances diverge
at high temperatures due to different exosphere temperatures from Jaccia
(1971). This interpolates and combines the results from Anderson et al.
(1986), Champion, Cole, & Kantor (1976), and Jaccia (1971).

The vertical profiles of thermodynamic and abundance quantities are not
measured comprehensively for all sites at all times. Consequently, they are
produced from theoretical models that take into account the atmospheric
chemistry (Smith 1982, WMO 1982, Brasseur & Solomon 1984, WMO 1986). We use
the Air Force Geophyiscal Laboratory (AFGL) constituent profiles (Anderson et
al. 1986) which comprehensively described thermodynamic and abundance profiles
for different latitudes and seasons from ground level to 120 km. The AFGL
profiles have five generic profiles for the equator (tropical), a mid-latitude
(45 degree) for winter and summer seasons, a sub-arctic (60 degrees) for
winter and summer seasons, and a US standard. The first 5 have pressure and
temperature profiles and abundanced profiles for water, ozone, N2O, CO, and
CH4. The US standard profile includes the same quantities but adds many other
molecules (CO2, O2, NO, SO2, NO2, NH3, HNO3, OH, HF, HCl, HI, ClO, OCS, H2CO,
HOCl, N2, HCN, CH3, H2O2, C2H2, C2H6, PH3, Ar, O, Ne, Kr). These additional
molecules are not significantly time-variable unlike the other molecules
(water, ozone, N2O, CO, and CH4). For this reason, we can reasonably use the
US standard abundance profile for the non-variable constituents regardless of
the location or season. Note for the last 4 elements above (Ar, O, Ne, Kr) we
had to augment the data by assuming that Ar, Ne, and Kr assume their surface
level value of Champion, Cole, & Kantor (1976) and the difference from the
total is pure $O$, since the Noble gas species were left out of Anderson et
al. (1986) because of their negligible chemical interations.

Above 120 km, the chemical differences of the atmosphere are considerable.
Here we can use empirical-based models of Jaccia (1971), which vary as a
function of the exosphere asymptotic temperature. The exosphere is relevant in
our work, since it plays an important role in scattering sun light at
twilight. The profiles in Jaccia (1971) predict temperature and pressure
profiles as well as molecular abundances of N2, O2, $O$, $Ar$, $He$, and $H$.
In Jaccia (1971) the exosphere ranges from 500 to 1600 K. The exosphere
temperature can either be specified by the user or we can estimate it based on
the time of the observation. Since the exosphere temperature is highly
correlated with solar activity, we predict it using a sinusoidal function with
the 11 year solar cycle period. The mean value is set to 1050 K with
sinusoidal amplitude of 550K. We also linearly interpolate the previous values
between 100 km and 120 km to asymptotically match the values of temperature,
pressure, and abundances at 120km in order to avoid discontinuities. The
profiles combining the compilations are shown in Figures 6, 7, and 8.

For the temperature, pressure, and abundance profiles of the variable
components, we use the following procedure to obtain realistic vertical
profiles for any location for any time in the year. Given the month and the
latitude and longitude of the observation, we choose the appropriate mean,
standard deviation, diurnal range, and diurnal range variation for the
pressure, temperature, and relative humidity. Then given the time of day of
the observation, we then determine the surface pressure, temperature, and
relative humidity by first gaussian sampling the mean value and the diurnal
variation using the climate data in §2.1. Then we determine the temperature,
pressure, and relative humidity by the following equations.

For the temperature variation, we use the Göttsche & Olesen (2001) model. The
temperature prior to sunset is given by

$T=T_{0}+\left(T_{1}-T_{0}\right)\cos{\left[\frac{\pi}{\omega}\left(h-14.2\right)\right]}$

The temperature after sunset is given by

$T=T_{0}+\left(T_{1}-T_{0}\right)\cos{\left[\frac{\pi}{\omega}\left(T_{s}-14.2\right)\right]}e^{-\frac{h-T_{s}}{k}}$

where the time of Sunset is $T_{s}$ and the constant, $k$ is defined below. We
calculate the day length, $\omega$, as given by Brock (1981) as

$\omega=\frac{2}{15}\arccos{\left(-\tan{\phi}\tan{\delta}\right)}$

where $\delta$ is the Solar declination and is given by

$\delta=23.45\sin{\left[\frac{360}{365}(284+d)\right]}$

where $d$ is the current number of days in the year. Then the constant, $k$,
is defined by

$k=\frac{\omega}{\pi}\arctan{\frac{\pi}{\omega}(h-14.2)}$

Similarly, the pressure, $P$, is defined by using a sinusoidal variation with
a 12 hour period using

$P=P_{\mu}+P_{\sigma}\cos{\left(2\pi\frac{h-10.5}{12}\right)}$

and the relative humidity is defined using a full day sinusoidal variation

$h=h_{\mu}+h_{\sigma}\cos{\left(2\pi\frac{h-3.0}{24}\right)}$

After selecting the temperature, relative humidity, and pressure, the profiles
described above will be selected with some modifications. We use a linear
combination of the five profiles with weights according to:

$\frac{1}{(\frac{\Delta T}{10})^{2}+\Delta\phi^{2}}$

where $\Delta T$ is the difference in temperature from the reference profile
at ground level and $\Delta\phi$ is the difference in latitude from the
reference profile. This allows us to crudely sample a profile that would more
closely resemble the location and season of the observation. The appropriate
linear combination of the temperature, pressure, and abundance profiles are
calculated. We temperature profile by a factor of

$\frac{T_{0}}{T_{ref}}\left(1-\frac{a}{a_{0}}\right)+\frac{a}{a_{0}}$

where $T_{ref}$ is the reference temperature at ground level of the five
profiles, $T_{0}$ is the chosen temperature, $a$ is the altitude, and $a_{0}$
is a reference altitude set to 120 km. The allows us to match the pressure,
temperature, and water fraction at ground level perfectly and then
asymptotically match the reference profile at high altitude. This is a
reasonable approximation to the weather conditions modifying lower layers of
the atmosphere, but greater stability at higher altitudes. We use an anologous
factor for the pressure and water fraction as well (replacing the temperature
above with the appropriate variable). For the water fraction, we first
calculate equivalent vapor saturation pressure by Wagner & Pruss (1993) and
Wagner, Saul, & Pruss (1994) using the selected temperature. We then use the
chosen relative humidity to predict a ground level abundance fraction using
the selected pressure and the saturation pressure. Finally, we use the same
non-variable molecules abundance profile, but adjust the overall weight so the
same proportion of non-water molecular ratios are obtained. We then calculate
the mean molecular weight as a function of altitude using values from PubChem
(2019). The density profile can be constructed assuming the ideal gas law, and
using the temperature, pressure, and mean molecular weight profiles.

We also predict the potential temperature profile. This is given by

$\theta=T\left(\frac{P_{0}}{P}\right)^{0.286}$

where $T$ is the temperature, $P$ is the pressure, and $P_{0}$ is the sea
level pressure. Since, potential is conserved for a parcel of fluid moving
adiabatically , the derivative of potential temperature ($\frac{d\theta}{dz}$)
can be use to predict the relative fraction of turbulence. Specifically, we
find

$\mu_{C_{N}^{2}}=50.0nk\frac{d\theta}{dz}$

$\sigma_{C_{N}^{2}}=2\times 10^{-18}n\mu_{C_{N}^{2}}$

where $n$ is the number density. The result is to be able to predict the
temperature, pressure, humidity, and turbulence intensity as a function of
altitude for any time and location on the globe while anchoring the results to
the theoretical calculations in the literature. As discussed above, these
quantities can be input by the user, but the ability to predict reasonable
values for any time and location will be useful.

###  2.3 Cloud Representation

Figure 9:  Cloud opacity simulations and comparisons. The top left panel shows
various realizations of the structure function from PhoSim compared with the
single SDSS measurements of Ivezic et al. (2007). The average of the PhoSim
curves is similar. The structure function normalization will change with
various levels of opacity, but we match the average structure function. The
top middle panel shows the cumulative distribution comparison with Essence by
Miknaitis et al. (2007). The other panels show the variety of opacity patterns
on 2 km scales setting the areas where there is no absorption to blue and high
absorption in white.

There are a number of aspects of clouds to fully capture their effect on
astronomical images. Their spatial structure can be characterized by a series
of frozen screens at various altitudes that drift during the exposure
according to the local wind vector. We first need to characterize the optical
depth to extinction (scattering or absorption) at each point. The morphology
of clouds patterns and types is not well characterized in a mathematical
representation in the literature. Therefore, we first attempt to represent
this in a complex manner to capture the major details for physics simulation.
We present the optical depth spatial position with a spatial power spectrum.
We find that a spatial power spectrum following the scale, $l^{\frac{1}{3}}$,
in 1-dimension that reproduces spatial statistics described below. This is
flatter than the Kolmogorov spectrum that follows, $l^{\frac{5}{3}}$, commonly
considered for atmospheric turbulence. We then generate integrated screens
having spatial power spectra following this spectrum. The amplitude therefore
has a unitless variance. We then represent the optical depth pattern in
proportion to the screen amplitude, however, an important physical aspect of
clouds is that there are many lines of sight where the optical depth is
negligible (i.e. they are scattered) that possibly represents a threshold
thermodynamic quantities required for condensation. To represent this, we make
the magnitude of the optical depth, $m$, equal to

$m(x,y)=max(\mu+\sigma s(x,y),0)$

where $s(x,y)$ is the screen pattern, $\mu$ is the mean, $\sigma$ is the
standard deviation. The maximum plays the role of a ramp or rectifier function
that then will make a significant fraction of the pattern of being cloud free.
The mean and standard deviation can then be predicted from historical
extintion distributions. The mean, which affects the cloud cover in the above
formula, is also significantly affected by the relative humidity with
considerable scatter. Therefore, we choose the cloud cover, $c_{c}$ according
to

$c_{c}=\alpha e^{-\frac{1-r}{\beta}}+\gamma g$

where $r$ is the relative humidity at ground level, $\alpha$ is the cloud
cover at 100 percent humidity, $\beta$ sets the fall off with humidity,
$\gamma$ is the amplitude of scattering using a Gaussian random number, $g$
which expands on the representation of Walcek (1994). This then will capture
the cloud cover correlations with humidity and therefore diurnal and seasonal
cycles represented above. The other variable besides cloud cover will be the
cloud depth, $c_{d}$, which we define as the average cloud extinction. We
choose this variable as a uniform distribution up to a maximum value,
$c_{d}^{M}$. We derive the mean, $\mu$ and standard deviation, $\sigma$ from a
truncated distribution from the cloud cover, $c_{c}$ and the cloud depth,
$c_{d}$, as

$\sigma=\frac{c_{d}}{\sqrt{2}fc_{c}+\frac{1}{\sqrt{2\pi}e^{-\sqrt{\pi}f}}}$

$\mu=\sqrt{2}\sigma f$

where $f$ is the inverse error function of $2c_{c}-1$. Thus, we can generate
the screens, $\mu$, and $\sigma$ by from two variables, $c_{c}$ and $c_{d}$.
Those variables can either be specified for a given simulation or predicted
randomly that also incorporates the humidity distributions.

To set the constants in the above equation, we use extinction distribution in
the ESSENCE survey (Miknaitis et al. 2007, Wood-Vasey et al. 2007, Narayan et
al. 2016) where the photometry across a multi-year survey was carefully
tracked in order to estimate accurate Supernovae photometry. The extinction
distribution of this survey has a narrow peak and a long tail where the clouds
dominate as shown in Figure 9 (C. Stubbs, priv. comm.). The narrow peak is
consistent with the affect of aerosols having $0.12\pm 0.05$ magnitudes of
extinction that is similar to what is expected from the aerosols distributions
described below for the Blanco site. The tail of the distribution can then be
used to predict the mean, $\mu$ and the standard deviation, $\sigma$ to be
used in the screen calculation above. Using the humidity distribution for the
Blanco site, we predict the cloud cover fraction and random predict the cloud
depth up to a maximum values. We find $\alpha=0.65$, $\beta=0.3$,
$\gamma=0.1$, $c_{d}^{M}=0.5$, fits the cumulative distribution well as shown
in Figure 9. The first three parameters are not strongly constrained by a
single data set, but this is also consistent with studies of cloud cover in
Walcek (1994) as well as the maximum global cloud cover being near 65% for
humid sites.

We then intergrate the screens by using the wind distributions for the screen
height as described above for the wind distributions. The height of cloud in
the simulation is also a complex topic. However, a simple prescription is to
place the screens at the approximate height where the water would achieve 100%
humidity given the relative humidity at ground level. Therefore, this is given
approximately by $8.5\mbox{km}\log{h}$ where $h$ is the relative humidity at
ground level for a pressure scale height of 8.5 km. We also split the cloud
distributions at the two closest levels to give some verticality to the
distribution and allow for different velocity vectors. The cloud depth of an
individual layer, $c_{i}$ has to be $c_{i}=1-\sqrt{1-c_{t}}$ to preserve the
total cloud cover, $c_{T}$. The cloud depth can be divided between the two
layers equally. Despite, the complexity the entire cloud representation, two
parameters (the cloud cover and cloud depth) capture most of the observational
characteristics.

We then consider the validation of the spatial statistics of the extinction
variation. The photometric variation of bright stars on small spatial scales
using the SDSS observations was studied in detail by Ivezic et al. (2007).
Here we can match the spatial statistics which tests both the input power
spectrum as well as the use of a ramp function described above. Figure 9 shows
the structure function of the Ivezic et al. (2007) results as well as the
structure function of the screens. On average there is good agreement with the
average structure function. Future work may allow us to look at structure
functions in a variety of conditions and sites as well as the average
photometry statistics, but this representation is sufficient given the current
set of measurements. The interaction of light with the clouds is described in
§3.

###  2.4 Aerosol Representation

Figure 10:  Global fit to the AERONET data (Holben et al. 1998, AERONET 2022)
for various sites across the global. The color corresponds with the optical
depth where blue indicates negligible optical depth and red indicates optical
depth of 0.5. Figure 11:  Global fit to the AERONET data (Holben et al. 1998,
AERONET 2022) for various sites across the global. The color corresponds with
the spectral index where blue indicates an index of -2 and red indicates an
index of 2. Figure 12:  Global fit to the AERONET data (Holben et al. 1998,
AERONET 2022) for various sites across the global. The color corresponds with
the spectral curvature where blue indicates a curvature of -0.5 and red
indicates a curvature of 0.5.

The final component of the atmosphere is aerosols that arise from a number of
sources. The dominant cause of large particles are from dust, pollution,
smoke, and sea salt. There is a wide variation from site to site and there is
likely a mixture at all locations. In addition, there is a significant
seasonal variation as well. To address this comprehensively, we use the
extensive database of AERONET sites that covers a series of measurements
across the globe (Holben et al. 1998, AERONET 2022). We use 1271 sites where
monthly data is collected on aerosol optical depths as a function of
wavelength for up to the last 25 years. For each monthly measurement at a
given site, we fit the optical depth as a function of wavelength to a power
law with possible curvature as:

$\tau=e^{p_{0}+p_{1}\log{\frac{\lambda}{500.0}}+p_{2}\left(\log{\frac{\lambda}{500.0}}\right)^{2}}$

where $\tau$ is the optical depth, $\lambda$ is the wavelength in $nm$, and
$p_{0}$ , $p_{1}$, and $p_{2}$ are the normalization, index, and curvature of
the power law. Typically, at least a half dozen wavelength measurements are
obtained so it is well-constrained. The measurements are shown in Figure 10,
11, and 12.

To separate out the aerosol components the following procedure is adopted. We
first found regions on the globe where each of the four components (dust,
pollution, smoke, and sea salt) was the dominant aerosol. We use the Sahara
desert for dust, the Amazon and Congo rainforests for smoke, isolated islands
far from continents for sea salt, and the largest population density areas for
pollution. For each case, the mean value of the three parameters can be
obtained as well as the standard deviation. In addition, there is sufficient
data to crudely estimate the altitude dependence of each component by
comparing different sites at various altitudes. We use altitude for every site
using the NOAA GLOBE (Global Land One-Kilometer Base Elevation) data set on 1
km scale (Hastings & Dunbar 1999). We use the population density for every
point on the globe from the CIESIN data set (CIESIN 2018). The altitude
dependence has an exponential scale height of 2.2 km for seasalt, 3.6 km for
dust, 6.6 km for smoke, and 1.7 km for pollution.

Then, we need a method to estimate the relative contribution of various
components at a given site. There are a variety ways to do this with a
detailed model, but a relatively simple prescription is to predict: 1) the sea
salt component in proportion to the nearby land area, 2) the smoke component
in proportion to the temperature and rainfall, 3) the dust component
negatively in proportion to temperature and rainfall, and 4) the pollution
component in proportion to the population density. The rainfall is predicted
by computing the average daily water vapor pressure times the humidity divided
by the density of water and gravitational acceleration. For a given site, we
consider the site and four points 1000 km away in each cardinal direction to
form a weighted sum.

$S=\sum_{i}\left(1-L_{i}\right)w_{i}e^{-2.56+0.29g}$
$D=\sum_{i}L_{i}\max[\left(-r_{i}+0.2t_{i}\right),0]\frac{1}{0.12}e^{-0.8+0.46g}$
$S_{m}=\sum_{i}L_{i}w_{i}\max[\left(r_{i}-0.2t_{i}\right),0]\frac{1}{0.06}e^{-1.83+0.8g}$
$P=\sum_{i}w_{i}p_{i}\frac{1}{0.14}e^{-1.28+0.4g}$

where $S$ is the seasalt, $D$ is the dust, $S_{m}$ is the smoke, $P$ is the
pollution, $p_{i}$ is the scaled population density, $t_{i}$ is the scaled
temperature, $r_{i}$ is the scaled rainfall, $L_{i}$ is a function that is 1
when the location is above sealevel and 0 otherwise, and $w_{i}$ is the
weighting function.

Finally, we need to determine the aerosol particle sizes in order to perform
Mie scattering in §4. To do this, we take the spectral index and curvature
determined for each of the four components. Then we consider particles having
a log normal distribution with a mean and variance. For each component, using
Mie scattering in §4, we choose combinations of mean and variance that are
consistent with the range of spectral indices and curvature for the
extinction. We also consider water droplets for the clouds and assume that the
extinction is approximately flat. We use the index of refraction for the
dominant molecular for each component: water (clouds), sodium chloride (sea
salt), carbon (smoke), silica (dust), and ammoniun sulfate (pollution). In
general, there are a range of particle size distributions that are consistent
with the extinction and may represents a combination of the true possible
distributions as well as some degeneracy due to the finite information in the
extinction spectrum. We find the following parameters to choose the log-normal
distributions with a correlation between the according to:

$\mu=p_{1}+p_{2}g_{1}$ $\sigma=p_{3}+p_{4}g_{2}+(\mu-p_{1})p_{5}$

where $p_{1}$ and $p_{2}$ set the mean particle size in microns, $p_{3}$ and
$p_{4}$ set the variation, $p_{5}$ induces a correlation, and $g_{1}$ and
$g_{2}$ are gaussian distributed random numbers. The particle distribution is
then chosen according to $e^{\mu+g\log\sigma}$, where $g$ is a gaussian random
number chosen for each particle. The parameters are given in the Table 1.
Alternatively, the aerosol optical depth for the four components can be
specified by the user for a particular location and observation time, but the
distributions by the calculation above should produce reasonable values for
any location and observation time.

Table 1:  Relationship between important environmental variables and image characteristics  Type | $p_{1}$ | $p_{2}$ | $p_{3}$ | $p_{4}$ | $p_{5}$
---|---|---|---|---|---
Water | 0.7 | 2.6 | 0.93 | 0.40 | -0.08
Sea Salt | -4.2 | 1.3 | 1.13 | 0.11 | -0.20
Dust | -3.6 | 1.2 | 1.16 | 0.10 | -0.22
Pollution | -4.3 | 1.3 | 0.99 | 0.04 | -0.21
Smoke | -4.9 | 1.2 | 0.95 | 0.06 | -0.22

##  3 Photon Interactions

The simulation of photons proceeds by following the photon through layers of
the atmosphere in a forward raytracing scheme (Peterson et al. (2015)). This
is implemented using first principle physics and acts on the atmosphere
representation described in §2. In each layer, we consider 1) the refraction,
2) the interaction with the layer of air through either Rayleigh or molecular
scattering and absoprtion, and 3) the possible interaction with the aerosols
in the bulk layer and cloud droplets in the cloud screens by Tyndall-Mie
theory. Following, we discuss each of these photon interactions. When dealing
with background light, a modified procedure that is more computationally
efficient is used, that is discussed separately below.

###  3.1 Refraction

The index of refraction in air, $n$, is given by the simplified formula of
Owens (1967)

$(n-1)10^{8}=(2371.34+\frac{683939.7}{130-\sigma^{2}}+\frac{4547.3}{38.9-\sigma^{2}})D_{s}$
$+(6487.31+58.058\sigma^{2}-0.71150\sigma^{4}+0.08851\sigma^{6})D_{w}$

$D_{s}=\frac{P_{s}}{T}\left[1+P_{s}\left(57.910^{-8}-\frac{9.325010^{-4}}{T}+\frac{0.25844}{T^{2}}\right)\right]$

$D_{w}=\frac{P_{w}}{T}\left[1+P_{w}\left(1+3.7e-4P_{w}\right)\left(-2.3732110^{-3}+\frac{2.23366}{T}-\frac{710.792}{T^{2}}+\frac{7.7514110^{4}}{T^{3}}\right)\right]$

where $T$ is the temperature in Kelvin, $P_{s}$ is dry air pressure, $P_{w}$
is the water pressure, and $\sigma$ is the inverse wavelength. Light will be
refracted throughout its trajectory as it passes through the atmosphere. It is
important to note that if the Earth were flat, then only the index of
refraction near the surface would matter. This is because plane-parallel
applications of Snell’s law essentially would cancel and only the final
refraction step would be relevant (Filippenko 1982). Since, the Earth is
curved this simplification is not valid, and the complete layer by layer
approach is necessary taking into account the full geometry (Stone (1996)). We
extend this approach in three important ways by using a numerical Monte Carlo
integration approach. One is that we have a range of temperature, pressure,
and water vapor profiles which are naturally that are handled naturally by
these methods. This allows us to assess the variation of atmospheric
dispersion under different conditions. Second, we can more accurately predict
dispersion particularly at large zenith angles, since other analytic
approaches involve expansions in tangent of the zenith angle which are known
to be catastrophically inaccurate at large zenith angles. Finally, we can also
perform dispersion for a individual photon at a time, which allows us to
accurately predict the range of emergent effects for non-monochromatic and
spectrally-varying spatialy-extended sources as we highlight in §4.1.

The atmospheric dispsersion is handled by first considering the direction to
the zenith position in the local coordinate system. We use the fact that the
refraction invariant, $r_{0}$, remains constant throughout the refraction
process (Fabritius 1878, Young 2004).

$r_{0}=(h+R_{E})n\sin{z}$

where $R_{E}$ is the radius of the Earth, $h$ is the height above sea level,
and $z$ is the zenith angle. The geometry can be written in differential form
so that the change in angle, $\Delta\phi$ at any layer due to refraction is
given by

$\Delta\phi=-\tan{z}\frac{\Delta n}{n}$

where $\Delta n$ is the change in index of refraction, $n$ is the local index
of refraction, and $z$ is the zenith angle (Young 2004). Before the photon
enters the atmosphere at height, $h_{0}$, the refraction invariant is simply,
$r_{0}=(R_{e}+h_{0})\sin(z_{0})$, and $\tan{z}$ can be substituted as

$\Delta\phi=-\frac{(R_{e}+h)\sin{z_{0}}}{\sqrt{{\left[n(R_{e}+h)\right]}^{2}+{\left[(R_{e}+h_{0})\sin(z_{0})\right]}^{2}}}$

Then the shift is performed in the direction of the zenith at each layer.
During the process both the index of refraction and the current height is
altered, so that the approximately curved trajectory is created through
numerical integration. We also allow for the direction of the angle to the
zenith to be modified during the exposure due to the field rotation, so
depending on the arrival time of the photon the refraction may be different.
In many applications, we also allow an option to simultaneously subtract off a
nominal refraction for a photon at the center of the field of view at a
nominal reference wavelength for the filter, so that the overall nominal
dispersion is cancelled as is customary for telescope operations. We describe
the effect of the dispersion in §4.1.

###  3.2 Rayleigh Scattering

The Rayleigh scattering cross-section can be written as

$\sigma=\frac{24\pi^{3}\left(n^{2}-1\right)^{2}}{\lambda^{4}N^{2}\left(n^{2}+2\right)^{2}}F$

where $n$ is the index of refraction (a function of wavelength and
thermodynamic variables), $N$ is the number density, $\lambda$ is the
wavelength, and $F$ is the polarization factor (e.g. van de Hulst (1957)). We
use the index of refraction described in §3.1 and the density given the
temperature, density, and water fraction at a given height when we evaluate
this cross-section. The calculation of Rayleigh scattering then resembles the
Tomasi (2005) calculation but using Owens (1967) for the index of refraction.

###  3.3 Molecular Interactions

We simulate molecular interactions by using the HITRAN data base (Rothman et
al. 2009, Rothman et al. 2005, Rothman et al. 2003, Rothman et al. 1998,
Rothman et al. 1992, Rothman et al. 1987, HITRAN 2014). We calculate the
cross-sections for molecular transmissions as described in Peterson et al.
(2015). In our new comprehensive representation, however, the molecular
interactions lead to absorption and re-emission rather than simply removing
photons (greenhouse effect). This is important for the simulation of the
scattered background light.

###  3.4 Tyndall-Mie Scattering

The aerosols and cloud droplets are large enough to require a full
electromagnetic treatment. Assuming the particles are quasi-spherical, the
interactions with photons can be described by Tyndall-Mie theory. This
performs a self-consistent electromagnetic calculation, and predicts the
absorption and scattering probability as a function of angle. The result
mainly depends on the dimensionless ratio, $x=\frac{2\pi r}{\lambda}$, where
$r$ is the size of the particle and $\lambda$ is the photon wavelength. It
also depends on the real and imaginary part of the index of refraction as a
function of wavelength. For water, we use the optical constants of Hale &
Querry (1973) for water, Sodium Chloride constants of Li (1976) for seasalt,
Silicon constants of Aspnes & Studna (1983) for dust, Ammonium Sulfate
constants of Cotterell et al. (2017) for pollution, and Carbon constants of
Phillip & Taft (1964) for smoke.

Numerically, we perform the Mie calculation for two-dimensional grid of values
of $x$ and $\lambda$. We use the Bohren-Huffman algorithm originally described
in Bohren & Huffmann (1983). This solves the electromagnetic boundary
conditions iteratively, and results in a scattering angular probability
distribution and absoprtion probability. We then save these probability
distributions and the total effective cross-section of both scattering and
absorption.

In order to properly simulate the photon interactions with these particles, we
need to construct a three-dimensional slab and perform photon Monte Carlo
techniques to handle multiple scattering properly. To do this, we simulate a
photon incident on the slab and consider its three-dimensional trajectory. We
then choose a particle from the log-normal distribution for the aerosol type
or cloud droplet distributions. We calculate the mean free path given that
particle size as

$\frac{f^{2}}{\tau\sigma^{2}}\frac{1}{2}$

where $\tau$ is the optical depth, $\sigma$ is the particle size, the factor
$\frac{1}{2}$ ensures that we allow for a cross-section four times more than
the geometric area, and $f$ is a numerical constant that we fix depending on
the log-normal distribution to achieve the correct optical depth. $f$ can be
determined by making

$1-e^{-\frac{\tau\sigma^{2}}{f^{2}}}$

match $1-\exp{-\tau}$ when the full log-normal distribution is considered for
$\sigma$. Then, the photon is shifted by a distance equal to the mean free
path following an exponential distribution.

Applying the photons angular trajectory, a new position in the slab is
calculated. Then, photons are removed if they are at a new position above the
slab and are therefore reflected. We roll a uniform random number to allow for
the possibility that the photons are neither scattered or absorbed since we
are assuming that the particle size are a factor of 4 above there geometric
area, and for most cases the effective cross-section is near a factor of 2
more than the geometric area. If they are still in the slab, then there is
some probability of photon absorption according to the calculated absorption
cross-sections. If they are not absorbed then, they will be scattered. We then
use the cumulative probability distributions given from the Mie probability
tables calculated earlier. These are formulated in terms of the cosine of the
polar angle, so the azimuthal angle is then chosen. Then, the three-
dimensional position is updated based on the new trajectory. Finally, the
process is repeated until the photon either emerges from the slab (either
forward or backwards) or is absorbed. Then, we will have a new perturbed
trajectory for photons that have successfully scattered off of an aerosol
particle or cloud droplet that is added to the existing trajectory before
interactiion with the layer. This process correctly handles multiple scatters
and matches the local optical depth which can vary considerably in the case of
clouds.

###  3.5 Background Methodology

Figure 13:  The calculation of the reflection fraction as a function of
optical depth from a series of Monte Carlo simulations. This is a Monte Carlo
simulation of a generic slab of material to predict the amount of light that
would reach the ground as a function of optical depth. At high optical depth
all light is reflected, whereas at low optical depth half is reflected. The
empirical fitting function is described in the text.

The simulation of the background is identical to the simulation of normal
sources with some important exceptions to improve efficiency. We start with
background light from the Sun, reflected sunlight off the Moon, airglow with
spatial variation, zodiacal light, and artificial light. We then modify the
simulation of light from that of normal sources to avoid simulating all
background at all angles with possible scattering before deciding which
photons would make it into a telescope aperture, which would be
computationally prohibitive. Essentially, we attempt to predict light from
which angles would make it into the the telescope aperture before any
significant computation. This would be trivial, if there was no scattering,
but scattering means the processes in the atmosphere have to be handled with
care. Essentially, this involves backward raytracing from the pupil plane.
This method is well known from computer graphics raytracing applications where
most of the interest is in raytracing prior to the pupil plane (Goldstein &
Nagel 1971). With the specific scattering in the atmosphere, we can perform
this efficiently with the following methods.

The first exception is for direct backgrounds at the top of the atmosphere
that includes zodiacal light and airglow. These background are approximately
flat, but have some small-scale spatial variation that is built into the
simulator by having a spatially-dependent emission pattern. Some fraction of
this light can be forward Mie scattered by the methodology described above and
would therefore be not detected in the field of view. However, a nearly
identical fraction of the light initially at angles outside the field of view
would scatter back in the field of view, and the net effect would be simply to
slightly blur the original spatial pattern of the emission. Therefore, to
simulate this we simply ignore forward Mie scattering for backgrounds of this
type to a very good approximation.

A second related consideration is to properly handle Rayleigh and molecular
scattering for these quasi-flat backgrounds. These processes scatter light to
much larger angles, and therefore remove a significant fraction of light when
the optical depth is high for a given wavelength. This is a straight-forward
numerical calculation and is simply dependent on the optical depth. From Monte
Carlo simulations, we find the fraction of light that is forward scattered, if
there is at least one scatter, as a function of optical depth is given by:

$0.5+\frac{5.28\tau+0.62\tau^{2}}{1+\frac{5.28\tau+0.62\tau^{2}}{0.5}}$

where $\tau$ is the optical depth. We performed Monte Carlo simulations for
both an isotropic scatterer (i.e. molecular transitions) as well as Rayleigh
scattering, and both were approximately the same. Therefore, at low optical
depth approximately half of the light that would continue to be scattered to
the ground if there is at least one scatter, whereas at high optical depth no
light would reach the ground. Therefore, we use this function to reduce the
opacity of the background light according to the optical depth, and therefore
properly predicts the total light fraction. The results are plotted in Figure
13.

A third consideration is to handle the scattered Moonlight and Sunlight. In
these cases, the scattered light itself is the main component of the
background. To simulate this light, consider an obervation at a point on the
sky at some particular angle. That angle will define a line segment from the
telescope aperture to the top of the atmosphere. Then, consider that a light
ray from the Sun or Moon will intersect that segment uniquely. Then, there
will be a particular angle that forms by the intersection of those line
segments. Therefore, in order for light to scatter it would need to scatter by
this particular angle. Therefore, to simulate this process we divide the
observing segment into a set of smaller segments corresponding to layers of
the atmosphere. At each segment we consider the probability that light could
scatter from a given process (Mie, Rayleigh, molecular). Then, if there is a
point where we obtain a successful scatter, we consider removing the photon by
attenuation prior to the scatter using a Monte Carlo approach. Then, we
continue the photon through the observing segment until it reaches the
atmosphere. It also could be scattered along this final path, so to properly
handle this we simply require that the photon only scattered once. This
complex propagation scheme where we essentially are stepping though two
partial line segments takes the same amount of computational time as a normal
photon. At the same time, we can treat all kinds of scattering properly and
deal with a complex altitude dependence of scatterers including clouds and
aerosols.

A final consideration is the simulation of artificial light (light pollution).
To simulate this, we simply simulate the light as if it came from the top of
the atmosphere. Then, we remove photons if there was not a scatter during the
propagation. This is equivalent to the correct physical simulation of light
starting at the bottom of the atmosphere and only considering photons where
they scattered backwards. In all of these four cases, we are basically
avoiding simulating any photons that do not have a significant chance of reach
the pupil and then performing the appropriate physical scattering along the
paths.

##  4 Observational Effects

There are a number of emergent effects that result from the photon
interactions that are imprinted on the images. Following, we discuss the
effect on astrometry, point spread ellipticity, photometry, point spread size,
and background.

###  4.1 Astrometry and PSF Ellipticity

Figure 14:  Differential chromatic refraction measured by the astrometric
positions of two monochromatic sources of 420 and 520 nm as a function of
zenith angle. We include points at various offsets from the field center by
1.4 degrees from a generic telescope to demonstrate that the refraction works
locally within a field. The points are compared to three analytic calculations
with good agreement at high altitudes (zenith angles less than 80 degrees).
Near the horizon the two plane parallel approximations fail since they use
expansions in the tangent to the zenith angle. The curved Earth two term model
of Stone (1996) has better agreement at high zenith angle, but still diverges
at high angles. PhoSim can work to arbitrarily low angles since it does not
involve any approximations, and can simulation telescopes at various
altitudes.

There are several aspects of refraction that affect the astrometry of the
field and ellipticity of effective point spread function. First, the overall
field is shifted by a significant amount as compared to the apparent position
if the atmosphere was not present. Second, there is a relative astrometric
shift of some objects with respect to other objects depending on the relative
position with respect to the zenith and the average wavelength of the detected
photon (i.e. differential chromatic refraction). Third, the detected photons
with shorter wavelength will be shifted more to the zenith direction than the
longer wavelength photons resulting in an elliptical shape. A fourth effect is
present when the astronomical source has a spatial distribution that varies
significantly based on wavelength (i.e. has a color gradient). This effect
then results in a morphological chromatic differential refraction where the
final distorted image is not simply described by the distortion by an
elliptical point spread function as the refraction affects the morphological
pattern in a wavelength-dependent manner.

The effect of differential chromatic refraction is shown in Figure 14. Here we
compare the predicted chromatic refraction with analytic models in the
literature. All four effects are essentially emergent effects due to the
single physical effect of the refraction of a photon depending on the
wavelength as it propagates through the atmosphere. Therefore, in a photon
Monte Carlo approach it is straight-forward to predict the complex emergent
behavior. In order to predict the lower altitude effect correctly, it is
critical to pursue a layer by layer approach including the correct geometry of
the curvature of the Earth otherwise the layer by layer refraction would
cancel and it would predict infinite refraction at the horizon.

###  4.2 Photometry

Figure 15:  Examples of absorption profiles as a function of wavelength in nm
at zenith. The variation is significant the variation of different components
(e.g water vapor, dust, clouds) including the height profile. The different
transmission functions are from completely different conditions (i.e.
different day). Figure 16:  Examples of absorption profiles as a function of
wavelength in nm at a zenith angle of 30 degrees. The variation is significant
the variation of different components (e.g water vapor, dust, clouds)
including the height profile. The different transmission functions are from
completely different conditions (i.e. different day).

Light is attentuated by essentially all of the physical photon interactions
described above. Rayleigh scattering decreases short wavelength light, atomic
and molecular absorption impacts specific frequencies, and Mie scattering and
absorption decreases transmission depending on wavelength, particle size, and
particle type. The calculation correctly calculates a local path length (i.e.
a local airmass), including the effect of the curvature of the Earth. Mie
scattering has a stronger wavelength dependence with smaller particles, so the
pollution, smoke, and seasalt distributions attenuate blue light
significantly. Dust has a mild wavelength dependence, whereas the wavelength
dependence of water droplets is minimal.

The atmosphere transmission is highly time variable, but this is predictable
from the self-consistent atmosphere representation pursued in this work.
Figure 15 and 16 show a range of transmission spectra in the optical band for
a zenith angle of 0 degrees and 30 degrees. Significant time variation is
present due to the variation of thermodynamic quantities (temperature,
density), atmosphere molecular constituents (especially ozone, water),
variation in cloud cover and depth, and variation in aerosols. The cloud cover
variation is clearly visible in a significant overall variation in the ceiling
of transmission functions. Interestingly, the aerosol variation is rather
significant in the short wavelength part of the spectrum and points to a
larger calibration challenge to u-band observations. Standard techniques of
adjusting zeropoints to match the photometry of bright calibration stars can
clearly correct for some of this variation. However, we can estimate the
possible photometric variation that is not easily correctable. If we simulate
two stars separated by 3 arcminutes the average relative photometric variation
is 20 millimags/$\sqrt{1+\frac{t_{exp}}{2s}}$. This would mean observations
with short exposures (between 10 seconds and a minute) could be limited to
photometric accuracy typically to $0.3$ to $0.8\%$ with large-scale empirical
calibration methods. Only more elaborate physically-based interpolations that
essential map cloud patterns and molecular spatial gradients could achieve
higher photometric accuracy.

###  4.3 Point-Spread-Function Wing

Figure 17:  Average atmosphere point spread function for 200 different
atmospheres mostly due to atmospheric turbulence. The red curve shows the
effect of Mie scattering of clouds whereas the blue curve removes the
scattering. The PSF wing starts to become significant beyond 1 arcsecond if
other sources of the PSF wing (such as mirror microroughness) are negligible.

The Mie scattering from both clouds and aerosols results in considerable
forward (small-angle) scattering. This is shown in Figure 17. If this
scattered light is at sufficiently small angles (i.e arcminutes), it will
appear as a wing on the point spread function. The scattering of particles
happens at all angles according to the Mie formalism, but some fraction of the
light will scatter at smaller angles. This will depend primarily on the ratio
of particle size to photon wavelength. The forward scattering happens
disproportionately for larger particles for a given wavelength in the
atmosphere, since the smaller particles asymptotically approaches the Rayleigh
limit (which has a modest $1+\cos{\theta}^{2}$ dependence) for a given
wavelength. Similarly, the shorter wavelength light will be scattered more for
a fixed particle size.

In our particle distribution for clouds and aerosols described above, the
largest particles will be the water droplets. Therefore, we the PSF wing is
dominated by scattered light from clouds. Similarly, the second largest
contributor to the PSF wing is the dust aerosols that have a distributions
that extend to larger sizes. The seasalt, pollution, and smoke scattering
results in a minimal distribution to the wing only from the largest particles.
In addition, for the same reason the scattering to small angles occurs
disproportionately for shorter wavelength light since the Mie forward
scattering occurs for larger ratios of the particle size over the wavelength.
For this reason, we predict a larger PSF wing in the UV part of the spectrum.
The wing from this scattering is the dominant contributor at larger angles
(arcminute). At smaller angles (several arcseconds), the PSF wing can have a
significant contribution from turbulence, diffraction, or mirror
microroughness.

Because the predicted PSF wing results primarily from water droplets in
clouds, there will be a correlation with the overall extinction from clouds,
since the larger extinction results in more scattering. This is intuitively
familiar as a fog-lamp effect where a car’s headlight in fog is both
attenuated and blurred (scattered) depending on the optical depth. 111Perhaps
the best artistic visualization of this is the exaggerated haloes around V.
Van Gogh’s stars. Thus, we also expect considerable variation in the PSF wing,
since the level of extinction can vary significantly in time as well as across
the field. There will then be an anti-correlation between an object’s
intensity and the PSF size (see §4).

###  4.4 Background

Figure 18:  The simulated background level as a function of some of the more
important variables: altitude of the observation (top left), moon altitude
(top right), distance to the Moon (bottom left) and the lunar phase (bottom
right). The different colors represent the u, g, r, i, and z intensities
(blue, light blue, green orange, red, respectively). The significant
complexity is evident as the intensity varies by four orders of magnitude. The
simulation is for a 1 m f/4 telescope located at a latitude and longitude of
(0,0) with 10 microns pixels and a exposure time of 15 seconds. Results can be
scaled for other telescopes and exposure times by multiplying by
$\frac{t_{exp}}{15~{}s}{\left(\frac{p_{size}}{10~{}\mu
m}\right)}^{2}{\left(\frac{4}{f_{\\#}}\right)}^{2}$. Figure 19:  Simulation of
the Moon with two observable effects. The cloud opacity is visible across the
Moon, and the Mie scattering halo is present near the edge of the Moon.

The overall background level is an emergent quantity that is predicted from a
combination of physical processes described above. Background, in particular,
is related to two major effects: 1) the deviation from an ideal astrometry
mapping patterns which distorts the incident background pattern and makes the
pixels cover a varying solid angle, and 2) the wavelength-dependent photometry
that has some spatial dependence in the direction of the horizon due to
atmospheric transmission. We also performed a simulation of the non-Lunar
background on 10 by 10 arcminute scales, and found a variation in excess of
the usual Poisson variation of 0.2%.

The overall background level is a function of several variables. The most
important are the altitude of the observation, the Lunar phase, the altitude
of the Moon, and the distance between the observation and the Moon. The cloud
patterns, aerosols, and the variation of atmospheric absorption lead to
additional variation. To assess this, we simulate 1000 random simulations
where: 1) the observation’s position is chosen randomly, 2) the Moon’s
position is chosen randomly (including below the horizon), 3) the Lunar phase
is chosen randomly, 4) the filter is chosen randomly, and 5) all other
variables (clouds, aerosols, atmospheric consituents) are chosen according to
the distributions described above. In Figure 18, we demonstrate the variation
in background by nearly 4 orders of magnitude. Predictably, the background
increases dramatically as the Moon is above the horizon and closer to the full
Moon phase. The background also increases dramatically at smaller lunar
distances due to increase of Mie scattering at smaller angles in comparison to
Rayleigh scattering. Interestingly, the non-Lunar background also shows a net
increase at low altitudes due to the relative increase in emission volume from
airglow. This is only partly modified by atmospheric attenuation which is more
visible in the bluer bands.

A demonstration of the complexity of the background is visible in Figure 19.
Here we simulate the Moon directly, and both the Mie scattering halo is seen
as well as the attenuation of direct light due to the spatial pattern of the
clouds.

##  5 Correlated Patterns

Figure 20:  Changes in differential astrometry, absolute photometry, relative
photometry, PSF wing fraction, background level, near Moonlight background
level, and far Moonlight background level as a function of the ground
temperature in Celcius. Figure 21:  Changes in differential astrometry,
absolute photometry, relative photometry, PSF wing fraction, background level,
near Moonlight background level, and far Moonlight background level as a
function of the relative humidity. Figure 22:  Changes in differential
astrometry, absolute photometry, relative photometry, PSF wing fraction,
background level, near Moonlight background level, and far Moonlight
background level as a function of the ground pressure in mmHg. Figure 23:
Changes in differential astrometry, absolute photometry, relative photometry,
PSF wing fraction, background level, near Moonlight background level, and far
Moonlight background level as a function of observation altitude. Figure 24:
Changes in differential astrometry, absolute photometry, relative photometry,
PSF wing fraction, background level, near Moonlight background level, and far
Moonlight background level as a function of cloud depth in magnitudes.

The net effect of the comprehensive implementation is that certain external
variables (temperature, pressure, etc.) produce correlated behaviours in the
emergent observable properties. Some of these correlations are well known, but
others are more subtle. To demonstrate this, we perform 150 sets of
simulations with the 5 most important external variables modified to 5
different values and all other variables randomly generated. The five external
variables are: pressure, temperature, humidity, altitude (zenith angle), and
cloud depth. For each of these variables, we varied the quantities between
extreme values that would likely be encountered at some point in an observing
run. Thus, temperature was varied between -10 and 40 degrees Celcius, humidity
was varied between 0 to 100%, pressure was varied between 740 and 780 mmHg,
altitude was varied between 30 and 90 degrees, and cloud depth was varied
between 0 and 0.2. In many real world cases, the variation may not be as large
(e.g. more cloudy nights might be avoided or switching to a new target at
different zenith angles). However, as we show below the effect on observables
is often close to linear, and the variation can then be scaled to a more
limited range. Other secondary variables (e.g. wind patterns, exosphere
temperature, aerosol depth, cloud patterns and cover fractions) are
marginalized in the 150 sets of simulations. Therefore, we have $3150=(4\times
5+1)150$ total simulations.

For each simulation, we measure several observable properites: differential
astrometry (non-ideal astrometric patterns), absolute photometry (overall
transmission), relative photometry (spectrally-dependent transmission), the
fraction of light in the PSF wing, overall nominal background level, moon
light intensity near the Moon, and Moon light intensity far from the Moon. We
define each metric as follows. Differential astronometry is measured by using
the difference in centroid of two sources with a monochromatic spectral energy
distribution at 350 nm and 550 nm. Absolute photometry is measured by using
the ratio of total counts of a source with a 550 nm SED. Relative photometry
is measured by the ratio of total counts from a source with a 350 nm SED and a
550 nm SED. The fraction of light in the PSF wing is measured by estimating
the source flux in rows with a maximum of less than 10 counts with the total
flux from a source with a 350 nm SED. The counts are collected in rows, since
the measurement involves highly saturated stars to see the wing pattern and
therefore the bleed trail is avoided. The nominal background is measured by
using the median counts using a V band filter. The Moon background is measured
in the same manner.

Table 2:  The systematic percent change of various image observables with the most important environmental variablesfor typical ranges of environmental variables. The error on each value is about 2%. The blank entries in the table correspond with a less than 2% correlation. Image Characteristic | Temperature | Humidity | Pressure | Altitude | Cloud Depth | Residual Variation
---|---|---|---|---|---|---
(Range) | (0,40∘C) | (0,100%) | (740,780mmHg) | (30,90∘) | (0,0.2mag) |
Differential Astrometry | -10 | 2 | 4 | -284 | - | 23
Absolute Photometry | 2 | -5 | - | 24 | -18 | 7
Relative Photometry | 8 | - | -3 | 64 | - | 59
Wing Fraction | - | - | - | - | 153 | 170
Background | -10 | -11 | - | -40 | -3 | 24
Near Moonlight | - | 24 | - | -47 | 3 | 30
Far Moonlight | - | 24 | - | -28 | 7 | 27

In Figures 20, 21, 22, 23, and 20, we the change in relative values of the
observable quantities with changes in the enviromental variables. In most
cases, the changes are approximately linear. In Table 2 we show the average
systematic percent change of observables with the range of the environmental
variables with accuracy to 2%. These would be the approximate predicted change
of observable quantities that an observer would likely encounter with a wide
range of conditions. We produce several well-known correlations that are well-
known: the correlation in altitude with absolute and relative photometry (24%
and 64%, absorption at higher airmass), the anti-correlation of altitude with
differential astrometry (-284%, atmospheric dispersion), the anti-correlation
of background with altitude (-40%), and the anti-correlation of absolute
photometry with cloud depth (-18%, basic cloud attenuation). However, several
subtle correlations are evident. One is that astrometry is affected by
temperature (-10%), humidity (2%), and pressure (4%) through the wavelength-
dependence of the index of refraction of air which is consistent with other
studies (Filippenko 1982, Stone 1996). Another is that absolute and relative
photometry can be affected by temperature (2% and 8%) by modifying the density
profile. This is also seen in the anti-correlation of relative photometry with
pressure (-3%). Humidity also reduce the absolute photometry through
absorption as expected (-5%).

The PSF wing fraction was quite strongly affected by cloud depth. This is
consistent with the estimate above that most of the PSF wing at arcminute
scales is due to water droplets, but the magnitude of the correlation is large
(153%). It may be possible to then improve photometry by estimating cloud
depth by looking at the PSF wing. The background level was affected by
temperature (-10%), humidity (-11%), and cloud depth (-3%) in addition to
altitude. It is difficult to provide a simple explanation for these
correlations, since the background level is an emergent property of the
calculation. In addition, scattered moonlight behaves differently than non-
lunar background and it varies with the angular distance to the moon. Here
there is a correlation in Lunar background with humidity (24%, 24%), anti-
correlation with altitude (-47%, and -28%), and correlation with cloud depth
(3%, 7%). These result implies that in moonless nights background will be
better at higher temperatures, higher humidities, and less clouds, but when
the Moon is visible lower humidities and more clouds are preferred.

There are also a significant number of subtle correlations of less important
variables implied by the physics implementation that are beyond the scope of
this work. The last column in Table 2 lists the residual variation of each
observable when the five most important environmental variables are held
fixed. This includes other variables implicit in the simulation. For instance,
we found that the twilight light intensity after sunset is most affected by
upper atmosphere scattering. The density in the upper atmosphere is affected
by the exosphere temperature, which is related to the solar cycle. Thus,
twilight and sunset intensities should follow the Solar cycle.

##  6 Conclusion, Validation, and Future Work

We presented a comprehensive photon Monte Carlo implementation of the
propagation of light through a self-consistent representation of the
atmosphere. The results are fully implemented in v5.7 of PhoSim which is
publically available at https://www.phosim.org. The new capabilities should be
sufficient for a number of applications that require realistic images that
have complex photometric, astrometric, PSF patterns that depend on observing
conditions. Some results will depend on the input data that includes the
climate historical data, the vertical thermodynamic profiles, the aerosol
patterns, and the the atomic data which will continue to be validated and
studied for non-astronomical applications. We have shown validation of the
cloud model with the structure function and opacity statistics and analytic
validation of differential chromatic refraction for low zenith angles.
Numerical validation of PhoSim is extensive and includes photon positional
errors to thousandths of pixels and transmission accuracies to one part in a
million. The emergent predicted results of the PSF wing, background levels at
various sites, transmission variation, and the correlations with environmental
variables can be studied by the astronomical community in a number of
different ways. Continued validation by the entire community will build trust
in the capabilities, and improve and make the results more robust. Future work
will refine the input data as future studies come to a greater understanding
of atmospheric details, and the code is constructed so users can input their
own atmosphere parameters to custom match specific observations or explore
dependencies. The capabilities included in this work should lead to a more
thorough understanding of the effect of the atmosphere on image properties and
many diverse measurements in astronomy.

We thank an anonymous referee for many helpful comments. JRP acknowledges
support from Purdue University. JRP acknowledges innumerable important and
helpful conversations with J. G. Jernigan. We acknowledge critical
conversations about layer-wise absorption with M. Wood-Vasey and cloud
photometry statistics with C. Stubbs. This work greatly benefited from the
AERONET, HITRAN, CIESEN, KGRD, and NCEP/NCAR databases.

## References

  * Ackermann et al. (2012) Ackermann et al. 2012, ApJS 203, 4.
  * AERONET (2022) AERONET (Aerosol Robotic Network) database, https://aeronet.gsfc.nasa.gov/
  * Anderson et al. (1986) Anderson, G. P., et al. 1986, Environmental Research Papers, 954.
  * Andersson, Peterson, & Madejski (2007) Andersson, K. A., Peterson, J. R., Madejski, G. 2007, ApJ 670, 1010.
  * Angeli et al. (2016) Angeli, G., Xin, B., Claver, C. et al., SPIE vol 9911, 18.
  * Aspnes & Studna (1983) Aspnes, D. E. & Studna, A. A. 1983, Phys. Rev. B. 27, 985.
  * Bard et al. (2016) Bard, D., Kratochvil, J. M., & Dawson, W. 2016, ApJ 819, 2.
  * Bard et al. (2013) Bard, D., Kratochvil, J. M., Chang, C. et al. 2013, ApJ 774, 1.
  * Bertin (2009) Bertin, E. 2009 Moemorie della Societa Astronmica Italiana 80, 422.
  * Beamer et al. (2015) Beamer, B., Nomerotski, A., Tsybychev, D. 2015, Jinst 10, 5027.
  * Bohren & Huffmann (1983) Bohren, C. F. & Huffman, D. R., “Absorption & Scattering of Light by Small Particles”, New York: Wiley.
  * Brasseur & Solomon (1984) Brasseur, G. & Solomon, S. 1984, Aeronomy of the Middle Atmosphere, D. Reidel Publishing, Dordrecht, Holland, 440.
  * Bretonniere (2022) Bretonniere, H., Kuchner, U., Huertas-Company, M. 2023, A&A 671, 102.
  * Britton (2004) Britton, M. C. 2004, Proc. SPIE 5497, 290.
  * Brock (1981) Brock, T. D. 1981, Ecological Modelling 14, 1.
  * Burke et al. (2019b) Burke, C., Aleo, P. Chen, Y., 2019, MNRAS 490, 3952.
  * Burke et al. (2019a) Burke, C., Peterson, J. R., Egami, E. 2019, JATIS 5, 3, 038002.
  * Carlsten et al. (2018) Carlsten, S. G., Strauss, M., Lupton, R. et al. 2018, MNRAS 479, 149.
  * Champion, Cole, & Kantor (1976) Champion, K. S. W., Cole, A. E., & Kantor, A. J. 1976, Standard and Reference Atmospheres, Handbook of Geophysics and Space Environments, Ed. A. S. Jursa, Chapter 14.
  * Chang et al. (2013a) Chang, C., Jarvis, M., Jain, B. et al 2013a, MNRAS 434, 2121.
  * Chang et al. (2013b) Chang, C., Kahn, S. M., Jernigan, J. G., 2013b, MNRAS 428, 2695.
  * Chang et al. (2012) Chang, C., Marshall, P. J., Jernigan, J. G., et al. 2012, MNRAS 427, 2572.
  * CIESIN (2018) Center for Internation Earth Science Information Network (CIESIN), Columbia University, 2018. Gridded Population of the World, Version 4: Population Revision 11, Palisades, NY, NASA SEDAC.
  * Cotterell et al. (2017) Cotterell, M. I., Willoughby, R. E. , Bzdek, B. R. et. al. 2017, Atmos. Chem. Phys. 17, 9837.
  * Davis et al. (2012) Davis, J. et al., SPIE 8443E, 1.
  * Dobke et al. (2010) Dobke, B. M., Johston, D. E., Massey, R. et al. 2010, PASP 122, 947.
  * Jaccia (1971) Jaccia, L. G. 1971, Smithsonia Astronomical Observatory Special Report, 332.
  * Göttsche & Olesen (2001) Göttsche, F. M. & Olesen, F. S. 2001, Remote Sensing Environment, 76, 337.
  * Ellerbroek (2002) Ellerbroek, B. L. 2002, Proceedings of the Beyond Conventional Optics, 229.
  * Fabritius (1878) Fabritius, W. 1878, Astron. Nachr. 93, 17.
  * Filippenko (1982) Filippenko, A. 1982, PASP 94, 715.
  * Goldstein & Nagel (1971) Goldstein, R. A. & Nagel, R., Simulation 16, 1.
  * Hale & Querry (1973) Hale, G. M. & Querry, M. R. 1973, Appl. Opt. 12, 555.
  * Hastings & Dunbar (1999) Hastings, D. A. & Dunbar, P. K., 1999. Global Land One-kilometer Base Elevation Digitial Elevation Model, Key to Geophysical Records Documentation (KGRD), NOAA, NGDC, Boulder, Co.
  * HITRAN (2014) High-Resolution Transmission Molecular Absorption Database (HITRAN) code 2014, http://www.cfa.harvard.edu/hitran.
  * Holben et al. (1998) Holben, B. N., Eck, T. F., Slutsker, I. et al. 1998, Rem. Sens. Environ. 66, 1.
  * Ivezic et al. (2007) Ivezic, Z. et al. 2007, AJ 134, 3, 973.
  * Jolissaint (2010) Jolissaint, L. 2010, JEOS Vol 5.
  * Kalnay et al. (1996) Kalnay, E. et al. 1996, Bulletin of the American Meterological Society 77, 473.
  * Lane, Glindemann & Dainty (1992) Lane, R. G., Glindemann, A, & Dainty, J. C. 2002, Waves in Random Media, Vol. 2, Issue 3.
  * Le Louarn (2002) Le Louarn, M. 2002, Proceedings of Beyond Conventional Optics, 217.
  * Li et al. (2016) Li, B., Li, G., Cheng, J. et al. 2016, RASS 16, 9.
  * Li (1976) Li., H. H. 1976, J. Phys. Chem. Ref. Data 5, 329.
  * Mandelbaum et al. (2014) Mandelbaum, R., Rowe, B., Bosch, J. et al. 2014, ApJ 212, 1.
  * Merlin et al. (2022) Merlin, C., Castellano, M., Bretonniere, H. 2023, A&A 671, 101.
  * Meyers & Burchat (2015) Meyers, J. E. & Burchat, P. R. 2015, JInst 10, 6.
  * Miknaitis et al. (2007) Miknaitis, G., Pignata, G., Rest, A. 2007, ApJ 666, 674.
  * Narayan et al. (2016) Narayan, G., Rest, A. Tucker, B. E., 2016, ApJS 224, 3.
  * Nie et al. (2021a) Nie., L, Li, G., Peterson, J. R. et al. 2021, MNRAS 503, 4436.
  * Nie et al. (2021b) Nie, L., Li, G., Fan, Z. et al. 2021, MNRAS 508, 3785.
  * Owens (1967) Owens, J. C. 1967, Appl. Opt. Vol 6., 1.
  * Peterson, Jernigan, & Kahn (2004) Peterson, J. R., Jernigan, J. G., Kahn, S. M. 2004, ApJ 615, 545.
  * Peterson, Marshall, & Andersson (2007) Peterson, J. R., Marshall, P. J., Andersson, K. 2007, ApJ 655, 109.
  * Peterson et al. (2020) Peterson, J. R. et al. 2020, ApJ 889, 2.
  * Peterson et al. (2019) Peterson, J. R. et al. 2019, ApJ 873, 98.
  * Peterson et al. (2015) Peterson, J. R. et al. 2015, ApJS 218, 14.
  * Phillip & Taft (1964) Phillip, H. R. & Taft, E. A., Phys. Rev. 136, 1445.
  * PubChem (2019) PubChem 2019, ncbi.nlm.nih.gov
  * Rothman et al. (2009) Rothman, L. S. et al. 2009, JQSQR 100, 533\.
  * Rothman et al. (2005) Rothman, L. S. et al. 2005, JQSQR 96, 139\.
  * Rothman et al. (2003) Rothman, L. S. et al. 2003, JQSQR 60, 665\.
  * Rothman et al. (1998) Rothman, L. S. et al. 1998, JQSQR 82, 5\.
  * Rothman et al. (1992) Rothman, L. S. et al. 1992, JQSQR 48, 469\.
  * Rothman et al. (1987) Rothman, L. S. et al. 1987, App. Opt. 26, 4058.
  * Rowe et al. (2015) Rowe, B. T. P., Jarvis, M., Mandelbaum et al. 2015, Astronomy & Computing 10, 121.
  * Sanchez et al. (2020) Sanchez J., Walter, C., Awan, H. 2020, MNRAS 497, 1.
  * Smith (1982) Smith, M. A. H. 1986, NASA Tech Mem 83289.
  * Stone (1996) Stone, R. C. 1996, PASP 108, 1051.
  * Thomas & Kahn (2018) Thomas, D. & Kahn, S. M., 2018, ApJ 868, 1.
  * Thomas et al. (2016) Thomas, S. J., Chandrasekharan, S., Lotz, P. et al. 2016, SPIE vol 9906, 3.
  * Tomasi (2005) Tomassi, C., Vitale, V., Petkov, B., Lupi, A., & Cacciari, A. 2005, Appl Opt 44, 16.
  * van de Hulst (1957) van de Hulst, H. C. 1957 , Light Scattering by Small Particles, New York: Wiley 1976, 176.
  * Walcek (1994) Walcek, C. 1994, Monthly Weather Review 122, 1024.
  * Walter (2015) Walter, C. 2015, JInst 10, 5015.
  * Wagner & Pruss (1993) Wagner, W. & Pruss, A., 1993, J. Phys. Chem. Ref. Data 22, 783.
  * Wagner, Saul, & Pruss (1994) Wagner, W., Saul, A., & Pruss, A. 1994, J. Phys. Chem. Ref. Data 23, 515.
  * WMO (1982) World Meterological Organization (WMO) 1982, Stratosphere 1981: Theory and Measurements, Report 11, NASA, Greenbelt, MD.
  * WMO (1986) World Meterological Organization (WMO) 1986, Atmospheric Ozone 1985, Report 16, Geneva, Switzerland.
  * Wood-Vasey et al. (2007) Wood-Vasey, W. M., Miknaitis, G., Stubbs, C. W. et al. 2007, ApJ 666, 694.
  * Xin et al. (2018) Xin, B., Ivezic, Z., Lupton, R. H. et al. 2018, AJ 156, 222.
  * Xin et al. (2015) Xin, B., Claver, C., Liang, M. et al. 2015, App. Opt., vol 54, issue 30, 9045.
  * Young (2004) Young, A. T. 2004, ApJ 127, 3622.