# A Spherical Shells Model of Atmospheric Absorption for Instrument
Calibration

Nicolas Donders  The University of Alabama in Huntsville, Huntsville AL  The
Department of Space Science and The Center for Space Plasma and Aeronomic
Research, Huntsville AL  Genevieve Vigil  NASA Marshall Space Flight Center,
Huntsville AL  Adam Kobelski  NASA Marshall Space Flight Center, Huntsville AL
Amy Winebarger  NASA Marshall Space Flight Center, Huntsville AL  Larry Paxton
Johns Hopkins Applied Physics Laboratory, Laurel MD  Charles Kankelborg
Montana State University, Bozeman MT  Gary Zank  The University of Alabama in
Huntsville, Huntsville AL  The Department of Space Science and The Center for
Space Plasma and Aeronomic Research, Huntsville AL

###### Abstract

We present a model for atmospheric absorption of solar ultraviolet (UV)
radiation. The initial motivation for this work is to predict this effect and
correct it in Sounding Rocket (SR) experiments. In particular, the Full-sun
Ultraviolet Rocket Spectrograph (FURST) is anticipated to launch in mid-2023.
FURST has the potential to observe UV absorption while imaging solar spectra
between 120-181 nm, at a resolution of $\mathcal{R}>2\cdot 10^{4}$ ($\Delta
V<\pm 15$ km/s), and at altitudes of between $\approx$ 110-255 km. This model
uses estimates for density and temperature, as well as laboratory measurements
of the absorption cross-section, to predict the UV absorption of solar
radiation at high altitudes. Refraction correction is discussed and partially
implemented but is negligible for the results presented. Absorption by
molecular Oxygen is the primary driver within the UV spectral range of our
interest. The model is built with a wide range of applications in mind. The
primary result is a method for inversion of the absorption cross-section from
images obtained during an instrument flight, even if atmospheric observations
were not initially intended. The potential to obtain measurements of
atmospheric properties is an exciting prospect, especially since sounding
rockets are the only method currently available for probing this altitude in-
situ. Simulation of noisy spectral images along the FURST flight profile is
performed using data from the High-Resolution Telescope and Spectrograph
(HRTS) SR and the FISM2 model for comparison. Analysis of these simulated
signals allows us to capture the Signal-to-Noise Ratio (SNR) of FURST and the
capability to measure atmospheric absorption properties as a function of
altitude. Based on the prevalence of distinct spectral features, our
calculations demonstrate that atmospheric absorption may be used to perform
wavelength calibration from in-flight data.

##  1 Introduction

Atmospheric absorption, in particular from Oxygen and Nitrogen molecules,
significantly limits our ability to observe the Ultraviolet (UV) spectrum from
the ground [Friedman et al., 1951]. Many teams use sounding rockets, weather
balloons, or high-altitude observatories for aeronomy or solar physics studies
to get above most of the atmosphere. Correcting for this absorption is
possible, however, some light will always be lost.

There are many models that predict the properties of Earth’s upper atmosphere.
The Naval Research Laboratory Mass Spectrometer and Incoherent Scatter
Exosphere 2000 (NRLMSISE-00) atmospheric model is typically used to obtain the
density of elements and their temperatures as they vary with altitude [Moussa
and Kuznetsova, 2022]. The most common absorption coefficient database, The
High-Resolution Transmission (HITRAN) molecular spectroscopic database is a
common resource for obtaining absorption coefficients, but it does not contain
absorption cross-section information for O2 in the UV spectral range [Gordon
et al., 2017]. The Max-Planck Institute for Chemistry in Mainz, Germany (MPI-
Mainz) Ultraviolet-Visible (UV-VIS) database [Keller-Rudek et al., 2014, 2013,
Sander et al., 2014] does contain this data, however, it requires significant
post-processing.

In our research, there are very few complete models for UV absorption and even
fewer that cover wavelengths below 300 nm [Smette et al., 2015, Kausch et al.,
2015, Jones et al., 2013, Noll et al., 2012, Lean, 1980]. To give an example,
He et al. [2019], Sun et al. [2018, 2022] aims to calculate the high-altitude
atmospheric emission effects due to “day-glow.” However, it does this by
looking through Earth’s upper atmospheric layers. This absorption model was
designed for the specific purpose of determining the optical depth between two
orbiting spacecraft and inferring the abundance of certain molecules in that
level of the upper atmosphere. Many aeronomy models account for the self-
absorption and re-emission of light, and more information can be found in
Degenstein [2000].

A pressing issue with all of these available models is that they are closed-
source and platform-specific programs. It would be extremely beneficial to the
atmospheric community if a general-purpose program could be developed that
pulled together these various resources. Thus, we have developed an open-
source platform-independent model with the primary goal of being more broadly
applicable. This goal is in line with other National Science Foundation (NSF)
goals such as the EarthCube project led by the National Center for Atmospheric
Research (NCAR) [Abernathey et al., 2017].

As a reference for building our model, a review paper by Meier (1991)
discusses many models of the time and presents a plot of the “Altitude of Unit
Optical Depth” (Meier [1991], Figure 2). Our initial aim was to match these
results by combining the NRL density model and the MPI-Mainz UV-VIS database,
providing a platform that can be adapted to fit the user’s needs.

The initial inspiration for this work was to understand atmospheric absorption
during sounding rocket flights, for example, the upcoming Full-sun Ultraviolet
Rocket Spectrograph (FURST). FURST was motivated by a distinct lack of
complete Vacuum Ultraviolet (VUV, 120-181 nm) spectra at high-resolution
($\mathcal{R}>2\cdot 10^{4}$, $\Delta V<15$ km/s) of the full solar disk. A
good example of the currently available data is SUMER, which has a good
resolution ($\mathcal{R}\approx 7-16\cdot 10^{4}$ from 67-161 nm) but did not
capture the full disk (1” x 300” slit) [Curdt et al., 2001, Peter, 1999, Peter
and Judge, 1999]. SUMER can perform a raster across the surface, but it takes
31 hours to do so, limiting the temporal resolution of such a measurement.
This and other slit-based spectrometers focus on narrow wavelength bands and
cannot capture the full disk of the Sun over small time scales. FURST will
provide not only an opportunity to observe the full disk of the Sun in such an
important and wide wavelength range but should observe UV absorption from
Earth’s upper atmosphere.

In addition, we intend to launch a set of instruments from higher latitudes.
The Hi-C Flare campaign will require waiting on the launch pad for hours in
order to launch at the opportune time to capture an active flaring region on
the Sun [Savage et al., 2021]. With this constraint, Poker Flat, Alaska is
required to be used. However, at this latitude and extended launch windows,
some atmospheric absorption will be inevitable. Predicting this effect will
assist our team in adjusting the launch windows (time of year/day) for this
exciting campaign.

In this paper, we discuss the development of an atmospheric absorption model
that can be used for various purposes. To test one application, we simulate
the potential impact on an instrument and utilize specific absorption features
as fiducial points to aid in spectral calibration. Additionally, given a high
spectral and radiometric resolution, we have the capability to analyze
spectral images of past sounding rockets to uncover the fundamental physical
properties of absorption-causing elements. In-situ measurements of the kind
and at these altitudes can only be done with sounding rockets.

In Section 2, we introduce the methods for building up the spherical shells
model, including refraction and atmospheric properties. We set up the model by
introducing the optical depth calculation for a curved atmosphere with
refraction. The model is then simplified for computation with imported
databases for atmospheric density and molecular absorption cross-section. We
finish the model by simulating a spectral image on our sounding rocket
instrument with these absorption effects. Finally, we discuss two possible
results from such an approach: spectral calibration and inversion of
atmospheric information from raw data. In Section 3 we will summarize these
results, and in Section 4 we discuss the potential uses for this technique and
future work.

##  2 Methods

###  2.1 Differential Path Length

To introduce the necessary components for building up our model, we start with
the optical depth. This term characterizes the amount of material that light
will have to pass through before arriving at the detector. One normally uses
Beer-Lambert’s Law to calculate the optical depth $\tau$, which requires
knowledge of the absorption cross-section $\sigma$ (cm2/molecule) and the
number density $\eta$ (molecules/cm3) as follows:

$\displaystyle I=I_{0}e^{-\tau},\
\tau(r,\lambda)=\int_{r}^{\infty}\sigma\left(r,\lambda\right)\eta\left(r^{\prime}\right)\frac{dr^{\prime}}{\mu}\implies\tau_{i,j}=\sum_{i^{\prime}=i}^{\infty}\sigma_{i^{\prime},j}\eta_{i^{\prime}}\Delta
L_{i},$ (1)

where $I$ is the intensity at the detector, $I_{0}$ is the intensity of the
source. The optical depth varies with height $r$ (km) and wavelength $\lambda$
(nm). Since we are developing a database-driven model, we switch to index
notation and label them as $i$ and $j$, respectively. We write the
transmission factor $T=I/I_{0}$ to measure how much light would be absorbed in
the specified altitudes and wavelengths.

It is the convention to denote the zenith angle $\theta$ in the form
$\mu\equiv\cos\theta$. In traditional models, the factor $\mu$ alone can
account for the change to optical depth due to changing zenith angle (see left
diagram in Figure 1). This is valid for a small angle $\theta$, however, we
want a more robust model that takes into account the curvature of the Earth’s
atmosphere, refraction, and changing altitude. With this in mind, we have
combined it with the model altitude resolution $\Delta r$ and will refer to it
as the “differential path length,” $\Delta L_{i}\equiv\Delta r_{i}/\mu_{i}$.

Solving for the differential path length is a primary contribution of this
model. The diagram for setting up the geometrical argument is shown in Figure
1. To add further complexity, we have considered the cross-section $\sigma$ to
be a function of temperature. This dependence becomes very important at the
altitude of unit optical depth, as our later results will highlight.
Temperature is of course dependent upon altitude, so we write
$\sigma\equiv\sigma\left(r,\lambda\right)$. We also attempt to account for the
impact of refraction to ensure we have covered the most significant factors.

Using these given variables, we can calculate the adjusted $\Delta L$ as
described in Figure 1. In this Figure, the spherical shells are shown on the
left, along with 3 triangles that highlight the relevant trigonometry.
Beginning at ground level, we solve for each $\Delta L$ along the path length
from the observer to the source. The triangles here indicate the necessary
angles and side lengths needed to solve for each $\Delta L_{i}$. As in
Equation 1, these non-uniform differential path lengths are what distinguish
this setup from the simple model with uniform $\Delta L=\Delta r/\mu$.

Figure 1: Cartoon of Geometry and Trigonometry to solve (not to scale). The
spherical shells have an initial zenith angle that depends upon the starting
point of the observation, $r_{0}$. The differential path length $\Delta L_{i}$
is labeled for a refraction/refraction-less setup. The next 3 diagrams show
the relevant geometries at this initial height, with the effect of refraction
exaggerated. The first triangle shows the necessary angles and side lengths in
order to solve for $\Delta L_{0}$. The next two show the progression in
solving for $\Delta L_{i}$ along each additional shell. For the proper optical
depth calculation in Equation 1, we must calculate all intermediate $\Delta
L_{i}$ along the optical path at each subsequent height $r_{i}$.

There are many ways to solve these triangles for $\Delta L_{i}$, but we have
done so in a way that requires the intermediate angle $\theta_{Li}$ in order
to prepare for the application of refraction. To solve for the differential
path length $\Delta L_{i}$, we employ simple trigonometric relations (see A)
to work out that $\Delta L_{0}=\sqrt{\left(r_{0}+\Delta
r_{0}\right)^{2}-r_{0}^{2}\left(1-\mu_{0}^{2}\right)}-r_{0}\mu_{0}$. This is
the solution for the first spherical shell at the initial altitude $r_{0}$.
Extending to the generalized solution for each subsequent shell $i$ in the
differential path length $\Delta L_{i}$:

$\displaystyle\Delta L_{i}$
$\displaystyle=\sqrt{\left(r_{0}+\sum_{i^{\prime}=0}^{i}\Delta
r_{i^{\prime}}\right)^{2}-\left(r_{0}+\sum_{i^{\prime}=0}^{i-1}\Delta
r_{i^{\prime}}\right)^{2}\left(1-\mu_{i}^{2}\right)}-\left(r_{0}+\sum_{i^{\prime}=0}^{i-1}\Delta
r_{i^{\prime}}\right)\mu_{i},$ (2) $\displaystyle\mu^{2}_{i}$
$\displaystyle=1-\left(\frac{r_{0}+\sum_{i^{\prime}=0}^{i-1}\Delta
r_{i^{\prime}}}{r_{0}+\sum_{i^{\prime}=0}^{i}\Delta
r_{i^{\prime}}}\right)^{2}\left(1-\mu_{i-1}^{2}\right),$ (3)

giving us a solution for the length of each discrete path between all shells.
This solution is fed into Equation 1 to yield the total absorption at each
particular wavelength for the starting altitude $r_{0}$. We can then adjust
the starting altitude $r_{0}$ and perform the calculations again to build up
the model for the absorption profile.

As mentioned in Section 1, a popular paper that describes atmospheric
absorption in these wavelengths is Meier [1991]. We set the date, time, and
location listed in Table 6 from that paper: March 21, 1980, 10:00 AM local
time, at the GPS coordinates for White Sands Missile Range (WSMR), NM. From
those details, we can query various databases and online resources for the
initial zenith angle $\mu_{0}\approx 0.7$ (from PySolar Stafford [2018]), and
the altitude and sea level for $r_{0}\approx 6379.4$ km (from USG [2021]). We
also load into the program the PyNRLMSISE-00 atmospheric model data at a
defined resolution $\Delta r=1$ km [Bender, 2022]. Higher spatial resolutions
can of course be specified as desired, but show no appreciable difference to
the results presented in this paper.

In the calculations above, $\theta_{Li}$ is solved at each step and depends on
the normal angle from the previous shell. This was done intentionally to
provide a natural progression to refraction correction.

###  2.2 Refraction

We refer to Born and Wolf [2013] for a rigorous derivation, showing that
refraction corrections are small for altitudes above 20 km and do not exceed
0.2∘. For completeness, we have developed the tools for applying refraction to
the intermediate angle $\theta_{Li}$ at the shell boundaries. For the
applications in this paper, it is not necessary to apply these corrections
Refer to A for more details.

From Snell’s law, there exists an equivalence on either side of the shell
“wall,” namely

$\displaystyle\mu_{i}^{2}=1-\left(\frac{n_{i-1}}{n_{i}}\right)^{2}\left(1-\mu_{i-1}^{2}\right),$
(4)

where n is the index of refraction of O2 and $\theta$ is the angle to the
normal of the shell boundary. Substitution into Equation 3 yields a correction
factor onto $\mu_{i}$. If the density is low ($\rho\leq 1$ g/cm-3), it can be
expected that the refractive index $n$ is near unity [Liu and Daum, 2008,
Phillips, 1920]. From the Lorenz-Lorentz relation [Buckingham and Graham,
1974], the refraction index can be expressed relative to the previous
refractive index as

$\displaystyle n_{i}=\frac{\eta_{i}}{\eta_{i-1}}\left(n_{i-1}-1\right)+1,$ (5)

where $\eta$ is the number density and $n$ is the refractive index. Using the
same low-density assumptions, it is known through the Cauchy formula [Cauchy,
1830, Born and Wolf, 2013] that the dispersion formula for the refractive
index as a function of wavelength can be expressed as

$\displaystyle n_{0}$ $\displaystyle=A\left(1+\frac{B}{\lambda^{2}}\right)+1,$
(6)

where constants A and B are tabulated from observations (see B).

A final complication related to refraction correction arises when we consider
the initial zenith angle compared with the “perceived” zenith angle. As shown
in Figure 1, the refracted path as perceived by the observer has a given angle
$\theta_{0}$ that is slightly higher than the actual zenith angle. As such, we
must correct the given initial zenith angle such that the curved refraction
path aligns with the location of the Sun. Since the correction for refraction
seems to be small for our applications in this paper, it is left for future
work. In the next section, we will put all of our functions together and
remove the recursive relationships.

###  2.3 Assembling the Model

In order to avoid computationally expensive functions, we simplify our
summations and nested expressions. Since we are using a uniform atmospheric
density model, apply a uniform $\Delta r$. Expanding the nested functions is
rather trivial and results in a much cleaner and faster computation (see B for
complete derivation).

$\displaystyle\Delta L_{i}$
$\displaystyle=\sqrt{\left[r_{0}+\left(i+1\right)\Delta
r_{0}\right]^{2}-\left(r_{0}+i\Delta
r_{0}\right)^{2}\left(1-\mu_{i}^{2}\right)}-\left(r_{0}+i\Delta
r_{0}\right)\mu_{i},$ (7) $\displaystyle\mu_{i}^{2}$
$\displaystyle=1-\left(\frac{r_{0}+\Delta r_{0}}{r_{0}+\left(i+1\right)\Delta
r_{0}}\right)^{2}\left(\frac{n_{0}}{n_{i}}\right)^{2}\left(1-\mu_{0}^{2}\right),$
(8) $\displaystyle n_{i}$
$\displaystyle=\frac{\eta_{i}}{\eta_{0}}\left(n_{0}-1\right)+1.$ (9)

This equation holds true for any $i\geq 0$. To avoid confusion with the
symbols, recall again that the refractive index is $n$ and the number density
is $\eta$. Among the initial conditions needed are the initial altitude
$r_{0}$ and the cosine of the zenith angle $\mu_{0}$, which are given based on
the location and time of day. As mentioned in Section 2.2, the refractive
index $n_{0}$ should be a function of wavelength. This correction is ignored
for now, and in fact refraction, in general, does not contribute to the
results presented in Section 3 and can be set to a constant of unity. In the
next sections, we discuss the models for number density $\eta$ and the
absorption cross-section $\sigma$.

###  2.4 Atmospheric Properties

In order to compute the expected absorption from the atmosphere, we need a
model of the density of the material causing the absorption. We could begin by
setting up the number density to follow the general exponential trend.
However, a complex data-driven model has already been developed by the Naval
Research Laboratory: the NRLMSISE-00 Model [Moussa and Kuznetsova, 2022,
Bender, 2022]. This model provides users with the temperature and density for
Hydrogen, Helium, Oxygen, O2, Nitrogen, N2, and Argon for specified dates and
locations as a function of altitude.

Figure 2: Densities of H, He, O, O2, N, N2, and Ar. These data show good
agreement with Meier [1991], Figure 30.

The data plotted in Figure 2 is retrieved at the time of year and for the
specified launch location from Meier [1991]. We can adjust the time of day as
well to see the effect this has on the final model. Similarly, the zenith
angle can be found in online tools such as the Solar Position Calculator
provided by NOAA [Cornwall et al., 2022, NOAA, 2022] or the python integrated
package “PySolar” [Stafford, 2018]. For this model, at solar noon $\mu\approx
0.7$ at ground level.

Figure 3: Absorption cross section for O2 between 120-181 nm (left) and
varying with Altitude between 120-130.4 nm. These data show agreement with
Meier [1991], Figure 27.

The final piece of the optical depth calculation is the absorption cross-
section. We show in Figure 3 the values for O2 available from the MPI-Mainz
UV-VIS Spectral Database [Keller-Rudek et al., 2013, Sander et al., 2014,
Keller-Rudek et al., 2014] for the specified spectral range. Cross-sections
are given in terms of cm2 per molecule. This plot roughly matches the order of
magnitude of Meier [1991], Figure 27, as well as a much older paper (Watanabe
et al. [1953], Figure 10).

The data set itself is based on laboratory measurements and detailed at
various wavelengths and temperatures. As such, the data are sparse and some
interpolation is needed to set up the continuous 2D array shown in Figure 3
(right). For this interpolation, we retain the given values while ensuring
smooth variation between any missing points. This is the simplest method that
preserves the general trend seen in the data. As noted later in Section 2.7,
the lack of complete data is another motivation for detailing this model. It
may be possible to fill in this missing information by looking at raw spectral
data taken in-situ and determining the absorption cross-section from that
measurement.

Note many interesting features here in Figure 3 (left). Firstly, the continuum
is nicely visible, with a few small peaks near the middle part of the range.
Given a large enough Signal-to-Noise Ratio (SNR), these stand-out peaks could
aid in calibration. The highest value, nearing the -16 order at around 125 nm,
will be useful as well. Finally, there are the Schumann-Runge bands at the far
end, beyond 175 nm. This messy region could provide fiducials for spectral
calibration [Yoshino et al., 1984].

Putting the cross-section data together with the number density is the final
piece to calculating the optical depth as it varies with height and wavelength
$\tau(r,\lambda)$. A useful representation of this result is the altitude of
“unit optical depth.” This is the location where $\tau=1$ in the exponential
term of Equation 1, and is shown as a function of altitude in Figure 4.

Figure 4: Altitude of Unit Optical Depth as calculated by Equation 1 and 7.

Shows agreement with Meier [1991], Lean [1987]

This result is in good agreement with Meier [1991], which also matched
previous results from Lean [1987]. This plot makes sense when compared with
the sharp drop in the density of O2 in Figure 2. Additionally, the right plot
in Figure 3 shows that this altitude has a significant shift in absorption
cross-section. These effects combine to produce a sharp transition at around
80-120 km. The sharp peaks near the start of the range will certainly provide
some calibration points, but the Schumann-Runge bands will most likely be too
low for this specific test. In the next section, we simulate a set of noisy
camera images for testing and analysis of this model.

###  2.5 Simulated Spectra

After calculating the optical depth, we can simulate an expected CCD signal by
combining the transmission value $T$ with the estimated photon flux $I_{0}$.
In this section, we will explain the method for generating a simulated solar
spectrum by using data from the High-Resolution Telescope and Spectrograph
(HRTS) as a baseline for solar activity.

The HRTS spectrum is well known and is one of the better databases for the UV
range [Brekke, 1993]. Other sources exist such as AIA and LISARD, which
present solar monitoring information for the general public [LIS, 2021]. As
the name suggests, HRTS has a high spectral resolution (0.005 nm). This is
necessary if we want to use it to predict what FURST will see, which is
expected to have a similar or better resolution.

The Flare Irradiance Spectral Model (FISM2) is also used. Described in
Chamberlin et al. [2020] and available on LIS [2021], this model puts together
actual measurements alongside an empirical model of solar variability to
produce a more complete data set of UV spectral irradiance at 0.1 nm
resolution.

The intensity of the HRTS Quiet Sun (QS) signal is shown in Figure 5. Since
the noise floor for HRTS is high, we also use the FISM2 model as a reference
from the same date (March 1980). The HRTS spectra shown here are originally
given in units of erg/s/sr/Å/cm2. We convert to ph/s/nm/cm2 by using the size
of the full disk ($9.352\cdot 10^{-3}$ sr), but an additional factor of
$\approx$ 1.4x is used in Figure 5 to scale the intensity to match the FISM2
model.

Figure 5: The solid line shows the HRTS spectra for the Quiet Sun (QS) scaled
by 1.4x to match the FISM2 model shown by the dashed line.

This difference is relatively small but could be attributed to the scaling of
HRTS to the full disk. Additionally, the QS spectra do not include additional
irradiance from Explosive Events (EEs), Active Regions (ARs), etc. Even so, we
prefer to use the HRTS spectra over the FISM2 model in this application
because it is given at a much higher spectral resolution. The best way to
obtain a full-disk high-resolution spectrum would be to measure it directly,
which is one of the primary motivations for developing an instrument like
FURST.

To simulate the spectral signals we need to fold both measurements into the
instrument functions for FURST. This conversion from the raw-data values to
Data Numbers (DNs) is done by having knowledge of the instrument, cameras,
filters, and other effects such as noise, and can be found in the C. Of note,
there will be a spectral plate scale of $\approx 0.005$ nm per pixel. This
defines the 1-dimensional linear wavelength-to-pixel mapping function used in
this paper to demonstrate the model’s use for calibration.

Additionally, specifying an exposure time requires us to implement the
expected flight profile, since a single exposure occurs at various altitudes.
For this setup, we have a 10-second exposure window at the apogee near 254 km,
where there is effectively no absorption, and another at the bottom end of the
flight between 115-132 km, where there is much more absorption.

Figure 6: Top: A simulated 2D noisy image from a 10-second exposure between
115-132 km. Bottom: The median signal level across each row (upper/solid black
line) with the estimated solar continuum (middle/red dashed line) and the
estimated noise floor (lower/solid red line) for the FURST instrument.

With these conversions in place, a noisy image is generated for a FURST camera
signal. The top image in Figure 6 is the 2D simulated signal for the window
between 115-132 km. The noise can be seen here, but the spectral features are
clear in the bottom plot. The top solid line is the FURST signal averaged
across all rows from the simulated 2D image. The red dashed line beneath it
contains the FISM2 simulated signal as described in Figure 5.

An important comparison to this available data is the estimated FURST noise
floor. To calculate this value, we take the various sources of noise from
Table 1 in C and work out the highest resulting output signal from an input
signal of zero photons. This works out to an output of approximately 2.5 DN
shown as the bottom solid red line in Figure 6.

On average, the Standard Error of the Mean (SEM) for these signals is much
less than 1% so it is not displayed here. It is also encouraging to see that
both signals are many orders of magnitude higher than the estimated noise
floor for FURST. This helps to point out the clear improvements being made
with the instrumental capabilities FURST will provide.

With the appropriate conversions and noise accounted for, we now have a method
for computing atmospheric absorption at the altitudes along a sounding rocket
flight profile. This will allow us to see if the instruments have enough SNR
and spectral resolution to be useful for calibration or inversion of
atmospheric properties.

###  2.6 Calibration

Before a spectral image from an instrument can be properly analyzed for solar
physics, it must first be calibrated. There are two different types of
calibration of interest to us here: radiometric and spectral. In this context,
radiometric calibration refers to determining the conversion of the signal at
each pixel from the computer signal in terms of data number (DN) back into
physical units (W cm-2 s-1 sr-1 Å-1). Likewise, spectral calibration involves
determining the fitting parameters for recasting the pixel numbers in terms of
their known wavelength values. One common parameter is the spectral plate
scale (in nm per pixel), although typically mapping functions are 2D nonlinear
functions.

Figure 7: Simulated signal from 10-second exposures, row-averaged. The left
plot shows the full wavelength range for one channel. The top right is the
transmission factor of a zoomed-in region, calculated as shown from the ratio
between the simulated data at the higher altitude ($S_{H}$) and the lower
altitude ($S_{L}$). The lower right plot compares this segment with a known
peak in absorption cross-section.

To use the identified absorption features for spectral calibration, we use a
spherical shells model to properly account for the levels of expected noise
and test whether or not we have sufficient resolution to be able to fit these
absorption features. In Figure 7, we show an example of this. Here we see the
simulated signals from a higher altitude $S_{H}$ and a lower altitude $S_{L}$.
In the right plots, we compare the ratio of those signals with the absorption
cross-section plot, using the corresponding altitude.

A simple Gaussian is used to map out corresponding pixel and wavelength pairs.
From just two of these pairs, we recalculate the spectral plate scale as to
within $\approx$0.5% error. More sophisticated wavelength mapping programs are
to be used in practice to develop our calibration beyond just solving for the
linear parameter. This is a promising result that shows we indeed have
sufficient SNR for a viable in-situ calibration, as well as that these
features are ubiquitous across the spectrum.

###  2.7 Inversion

Looking back to the models used to generate the spectral images in the
previous sections, the NRL density model is well-trusted but the absorption
cross-section database is sparse. This motivated an investigation into a
specific application for sounding rockets. It may be possible to invert the
atmospheric absorption from spectral images taken during a flight.

From Equation 1, the density and absorption cross-section is folded into the
integral that yields optical depth. If spectral images are taken over a
relatively small range of altitudes, temperature variations will be small and
the cross-section $\sigma$ can be considered independent of height. Thus, the
term can be separated following

$\displaystyle T_{i,j}$
$\displaystyle=e^{-\sigma_{j}\sum_{i^{\prime}=i}^{\infty}\eta_{i^{\prime},j}\Delta
L_{i}}\implies\frac{-\ln{T_{i,j}}}{\sum_{i^{\prime}=i}^{\infty}\eta_{i^{\prime},j}\Delta
L_{i}},=\sigma_{j},$ (10)

where the transmission factor $T$ must be found at altitude $i$ and wavelength
$j$. If a spectral image is taken near apogee then absorption is effectively
negligible. Taking the ratio of an observation from a lower height with one
taken at the maximum height will give a proxy of $T$ at that lower altitude
$i$. Then we can apply Equation 10 to invert the absorption cross-section
using our model.

We can test this with the simulated results from FURST described in Section
2.6 and shown in Figure 7 (left). We take the simulated noisy spectral images
taken at the corresponding highest ($\approx$ 254 km) and lowest ($\approx$
115-132 km) altitudes. Aside from the noise, the main difference between these
two spectral images can be attributed to the atmospheric absorption effect. We
then divide the lower image from the upper images, as in Equation 10, and
recover a “measured” optical depth.

The results of this test are shown in Section 3, Figure 9. As an alternative
test, we could place more trust in the cross-section data and attempt to solve
for a density profile. The results of this experiment are shown in D.

##  3 Results

We first show a comparison of the differential path length as calculated by
the “simple” model (where $\Delta L=\Delta r/\mu$) with the result from the
spherical-shells model (where $\Delta L$ is calculated as in Equation 7). An
example is shown in Figure 8, generated at varying initial zenith angles
($\mu_{0}=$ 0.2, 0.4, and 0.6) to highlight the effect. We can see the
difference between the use of uniform and non-uniform differential path
lengths from the simple model and the shells model, respectively.

Figure 8: Comparison of the Differential Path Length $\Delta L$ at different
initial zenith angles $\mu$ in shades of gray. The bottom red dashed line
indicates the altitude resolution $\Delta r$ of 1 km in this example. The
simple model would calculate $\Delta L$ as $\Delta r/\mu$, which corresponds
to the initial point on the left, at 0 km.

The simple model would estimate a uniform differential path length as $\Delta
r/\mu$, which corresponds to the initial points (at 0 km) in Figure 8.
However, the spherical shells model allows for this to be a non-uniform value.
At higher shells, the angle between the shell and the path approaches normal
($\mu\longrightarrow 1$), forcing the value of the differential path length to
approach the shell spacing ($\Delta L\longrightarrow\Delta r$) (bottom/red
dashed line). To help understand the effect of this, we can compare the sum of
the $\Delta L$ values between each model. When $\mu\approx 0.2$, the shells
model puts this value at about half of the simple model. As we approach a
zenith angle greater than $0.6$ the difference drops to less than 10%. This
model is still in-line with the previous work by Meier [1991] since the time
of day provided a high zenith angle. These results, however, predict slightly
less absorption overall for our flight path than simpler models.

Next, as can be seen in Figures 6 and 7, once the signal is transmitted
through the instrument we have enough SNR to tell where the absorption lines
are. These fiducials provide calibration points to verify the pixel-to-
wavelength mapping of the instrument to within 0.5%, even though we only
considered 2 lines in this basic example. Furthermore, this calibration only
considered one wavelength channel for FURST. By coupling more fiducials and a
broader range we hope to add significantly to the wavelength calibration of
the data through this in-situ method.

Finally, we solve for the absorption cross-section $\sigma$ required to
produce the spectra seen by an instrument like FURST. As described in Section
2.7, this inversion required taking the ratio of the spectral image taken at
apogee with one from a lower height (see Figure 7) and is formulated in
Equation 10.

Figure 9: Recalculated absorption cross-section from simulated images (red)
compared with the known data (black). Both are smoothed and scaled for easier
comparison.

In Figure 9, the black and red lines represent the known and the computed
absorption cross-section $\sigma$, smoothed and scaled for ease of comparison.
We clearly reproduce the values quite well, but it is also apparent where
noise has made an impact. There are a few deep lines (around 122 nm and 127
nm) that are a result of this.

Additionally, the absolute values are a bit off. From our testing, this can be
attributed to the temperature change at the altitudes where this simulated
signal occurs There is a sharp change in the temperature, and thus the
absorption cross-section values. To better reproduce the laboratory values, we
can employ a more nuanced calculation. For example, we have utilized a
weighted average based on the flight profile. However, more work must be done
to understand what is missing.

Regardless, the primary features of interest are the peak locations rather
than their absolute values, and this result shows that they are almost
completely recovered. There is also a tendency for higher values to have less
error. This is because there is less noise associated with the spectral image
signal from which it is calculated. It may not be possible to perform this
calculation from instruments with a low SNR. This further stresses the
importance of radiometric calibration, and the added value such an effort
brings.

This method shows a broadly applicable technique for the inversion of
atmospheric properties from existing data. It may help significantly to fill
gaps in the predicted absorption points in our databases. We have set up this
model to look at molecular oxygen (O2) since this is the primary contributor
to UV absorption. However, other molecules can be explored with this model
that may be of greater interest to atmospheric scientists (e.g. Ozone,
Methane, Aerosols, etc.).

##  4 Discussion

This updated model has been shown to accurately reproduce the expected
atmospheric absorption predicted by previous models using modern techniques.
Our results match previous reports well while accounting for the curvature of
the Earth and its effects on a non-uniform differential path length. For
zenith angles above $\mu>0.6$, a simple uniform differential path length will
change the optical depth by a factor of less than 10%. This model is
especially suited to applications with a low zenith angle, such as
heliophysics tests performed outside of solar noon, at high latitudes, or with
a low flight apogee. With knowledge of the proper absorption effects, this
model can assist in extending the launch capabilities for such missions.

This model’s biggest benefit for our team, though, is in improving our
understanding of the UV absorption effect in FURST data products. It will also
enable us to set limitations or requirements for future launch conditions.
Missions with a desire to measure or avoid UV absorption can benefit from such
understanding.

There are of course many improvements that can be made. Extrapolating the
given density data set out above 1000 km, and increasing the resolution, could
improve the resolved density profile or better predict the absorption amount.
There also remain other atmospheric effects that we anticipate being useful
for the calibration and inversion of the FURST data products. Hydrogen and
Oxygen resonant absorption adds additional absorption features that can be
utilized for calibration and cross-section inversions. A thermally expanded
atmosphere due to solar sub-storms, or other such responses to flare events,
will have the effect of increasing the altitude of unit optical depth [Jacchia
and Slowey, 1966].

The largest improvement to be made is with the sparseness of the absorption
cross-section data. Around the upper heights of the rocket flight, the
temperature does not vary all that much. However, the NRL model shows
significant variation in temperature at the altitudes corresponding to the
beginning and end of the shutter open time. We have attempted to use the
available data and have adjusted our code to accept whatever new data becomes
available. Rather than rely on laboratory measurements, we could look into
calculating the theoretical spectral width of the O2 absorption features. Not
only should spectral widths change with temperature, but some lines may only
appear at certain temperature ranges. Understanding the physics behind
calculating the resonant absorption lines will likely lead to an understanding
of how to solve this problem.

Additional work on FURST involves improved radiometric calculations for
tracking errors throughout the calibration system and flight instrument and
has been discussed in a separate paper [Vigil et al., 2021]. With these high-
resolution capabilities, the noise floor of the FURST instrument will be well
below the expected solar continuum, making the absorption profile easier to
resolve. Following our spherical shells model, in-situ wavelength calibration
and radiometric inversion of the atmospheric absorption cross-section should
be possible.

##  5 Acknowledgements

The primary author wishes to thank their advisors, mentors, and colleagues
that helped with countless meetings, advice, and support along the way: Amy,
Charles, Larry, Gen, Adam, Gary, Michael, Keyvan, Clayton, Juan, all of the
staff at UAH for the funding support, his wife Crystal, but most importantly,
God (Matthew 6:33-34). This material is based upon work supported by the NSF
EPSCoR RII-Track-1.2a (Non-invasive plasma diagnostics for LTP) Cooperative
Agreement OIA-1655280. Any opinions, findings, conclusions, or recommendations
expressed in this material are those of the authors and do not necessarily
reflect the views of the National Science Foundation.

## References

  * Friedman et al. [1951] H. Friedman, S. W. Lichtman, and E. T. Byram.  Photon counter measurements of solar x-rays and extreme ultraviolet light.  _Phys. Rev._ , 83:1025–1030, Sep 1951.  10.1103/PhysRev.83.1025.  URL link.aps.org/doi/10.1103/PhysRev.83.1025.
  * Moussa and Kuznetsova [2022] M Moussa and M Kuznetsova.  Nrlmsis atmosphere model.  Technical Report Version 1.0 and 2.0, NASA Community Coordinated Model Center, 2022.  kauai.ccmc.gsfc.nasa.gov/instantrun/msis.
  * Gordon et al. [2017] Iouli E Gordon, Laurence S Rothman, Christian Hill, Roman V Kochanov, Y Tan, Peter F Bernath, Manfred Birk, V Boudon, Alain Campargue, KV Chance, et al.  The hitran2016 molecular spectroscopic database.  _Journal of Quantitative Spectroscopy and Radiative Transfer_ , 203:3–69, 2017.
  * Keller-Rudek et al. [2014] H Keller-Rudek, G K. Moortgat, R Sander, and R Sörensen.  The mpi-mainz uv/vis spectral atlas of gaseous molecules of atmospheric interest.  Technical report, Earth System Science Data, 2014.  uv-vis-spectral-atlas-mainz.org/uvvis/cross_sections/Oxygen/O2.spc.
  * Keller-Rudek et al. [2013] Hannelore Keller-Rudek, GK Moortgat, Rolf Sander, and R Sörensen.  The mpi-mainz uv/vis spectral atlas of gaseous molecules of atmospheric interest.  _Earth System Science Data_ , 5(2):365–373, 2013\.
  * Sander et al. [2014] Rolf Sander, Hannelore Keller-Rudek, Geert Moortgat, and Rüdiger Sörensen.  The mpi-mainz uv/vis spectral atlas of gaseous molecules of atmospheric interest.  In _EGU General Assembly Conference Abstracts_ , page 1753, 2014.
  * Smette et al. [2015] Alain Smette, Hugues Sana, S Noll, H Horst, W Kausch, S Kimeswenger, M Barden, C Szyszka, AM Jones, A Gallenne, et al.  Molecfit: A general tool for telluric absorption correction-i. method and application to eso instruments.  _Astronomy & Astrophysics_, 576:A77, 2015.
  * Kausch et al. [2015] W Kausch, S Noll, A Smette, S Kimeswenger, M Barden, C Szyszka, AM Jones, Hugues Sana, H Horst, and F Kerber.  Molecfit: A general tool for telluric absorption correction-ii. quantitative evaluation on eso-vlt/x-shooterspectra.  _Astronomy & Astrophysics_, 576:A78, 2015.
  * Jones et al. [2013] Amy Jones, Stefan Noll, Wolfgang Kausch, Cezary Szyszka, and Stefan Kimeswenger.  An advanced scattered moonlight model for cerro paranal.  _Astronomy & Astrophysics_, 560:A91, 2013.
  * Noll et al. [2012] S Noll, W Kausch, M Barden, AM Jones, C Szyszka, S Kimeswenger, and J Vinther.  An atmospheric radiation model for cerro paranal-i. the optical spectral range.  _Astronomy & Astrophysics_, 543:A92, 2012.
  * Lean [1980] Judith L Lean.  _Atmospheric ultraviolet absorption spectroscopy_.  PhD thesis, University of Adelaide, Dept. of Physics, 1980.
  * He et al. [2019] Weiwei He, Kuijun Wu, Yutao Feng, Di Fu, Zhenwei Chen, and Faquan Li.  The radiative transfer characteristics of the o2 infrared atmospheric band in limb-viewing geometry.  _Remote Sensing_ , 11(22):2702, 2019.
  * Sun et al. [2018] Kang Sun, Iouli E Gordon, Christopher E Sioris, Xiong Liu, Kelly Chance, and Steven C Wofsy.  Reevaluating the use of o2 a 1$\delta$ g band in spaceborne remote sensing of greenhouse gases.  _Geophysical Research Letters_ , 45(11):5779–5787, 2018.
  * Sun et al. [2022] Kang Sun, Mahdi Yousefi, Christopher Chan Miller, Kelly Chance, Gonzalo González Abad, Iouli E Gordon, Xiong Liu, Ewan O’Sullivan, Christopher E Sioris, and Steven C Wofsy.  An optimal estimation-based retrieval of upper atmospheric oxygen airglow and temperature from sciamachy limb observations.  _Atmospheric Measurement Techniques Discussions_ , pages 1–37, 2022\.
  * Degenstein [2000] Douglas Arthur Degenstein.  _Atmospheric volume emission tomography from a satellite platform._ The University of Saskatchewan, 2000.
  * Abernathey et al. [2017] Ryan Abernathey et al.  Collaborative proposal: Earthcube integration: Pangeo: An open source big data climate science platform.  _GSTDTAP_ , 2017.
  * Meier [1991] RR Meier.  Ultraviolet spectroscopy and remote sensing of the upper atmosphere.  _Space Science Reviews_ , 58(1):1–185, 1991.
  * Curdt et al. [2001] W Curdt, P Brekke, U Feldman, K Wilhelm, BN Dwivedi, U Schühle, and P Lemaire.  The sumer spectral atlas of solar-disk features.  _Astronomy & Astrophysics_, 375(2):591–613, 2001.
  * Peter [1999] H Peter.  Analysis of transition-region emission-line profiles from full-disk scans of the sun using the sumer instrument on soho.  _The Astrophysical Journal_ , 516(1):490, 1999\.
  * Peter and Judge [1999] Hardi Peter and PG Judge.  On the doppler shifts of solar ultraviolet emission lines.  _The Astrophysical Journal_ , 522(2):1148, 1999\.
  * Savage et al. [2021] S. Savage, A. Winebarger, L. Glesener, L. Golub, P. Chamberlin, Hi-C Flare Team, Foxsi-4 Team, and Snifs Team.  The First Solar Flare Sounding Rocket Campaign and Its Potential Impacts for High Energy Solar Instrumentation.  In _American Astronomical Society Meeting Abstracts_ , volume 53 of _American Astronomical Society Meeting Abstracts_ , page 313.15, June 2021\.
  * Stafford [2018] Brandon Stafford.  pysolar, September 2018.  URL doi.org/10.5281/zenodo.1461066.
  * USG [2021] Nov 2021.  URL nationalmap.gov/epqs/.
  * Bender [2022] Stefan Bender.  Pynrlmsise-00, January 2022.  URL pynrlmsise00.readthedocs.io.
  * Born and Wolf [2013] Max Born and Emil Wolf.  _Principles of optics: electromagnetic theory of propagation, interference and diffraction of light_.  Elsevier, 2013.
  * Liu and Daum [2008] Yangang Liu and Peter H Daum.  Relationship of refractive index to mass density and self-consistency of mixing rules for multicomponent mixtures like ambient aerosols.  _Journal of Aerosol Science_ , 39(11):974–986, 2008.
  * Phillips [1920] P Phillips.  The relation between the refractivity and density of carbon dioxide.  _Proceedings of the Royal Society of London_ , 97(684):225–240, 1920.
  * Buckingham and Graham [1974] Amyand David Buckingham and C Graham.  The density dependence of the refractivity of gases.  _Proceedings of the Royal Society of London_ , 337(1609):275–291, 1974.
  * Cauchy [1830] L. Cauchy.  Sur la dispersion de la lumiere.  _bull. des. sc. maht._ , 14:9, 1830.  URL cir.nii.ac.jp/crid/1570572699363130880.
  * Cornwall et al. [2022] C Cornwall, A Horiuchi, and C Lehman.  Noaa esrl solar position calculator.  Technical report, NOAA Earth System Research Laboratories Global Monitoring Laboratory, 2022.  gml.noaa.gov/grad/solcalc/azel.html.
  * NOAA [2022] NOAA.  Noaa esrl solar calculator.  Technical report, NOAA Earth System Research Laboratories Global Monitoring Laboratory, 2022.  gml.noaa.gov/grad/solcalc/.
  * Watanabe et al. [1953] Kenichi Watanabe, Murray Zelikoff, and Edward C Inn.  Absorption coefficients of several atmospheric gases.  Technical report, Air Force Cambridge Research Labs Hansom AFB MA, 1953\.
  * Yoshino et al. [1984] K Yoshino, DE Freeman, and William H Parkinson.  Atlas of the schumann–runge absorption bands of o2 in the wavelength region 175–205 nm.  _Journal of physical and chemical reference data_ , 13(1):207–227, 1984.
  * Lean [1987] Judith Lean.  Solar ultraviolet irradiance variations: A review.  _Journal of Geophysical Research: Atmospheres_ , 92(D1):839–868, 1987.  doi.org/10.1029/JD092iD01p00839.  URL agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/JD092iD01p00839.
  * Brekke [1993] P Brekke.  An ultraviolet spectral atlas of the sun between 1190-1730 a.  _The Astrophysical Journal Supplement Series_ , 87:443–450, 1993.
  * LIS [2021] Lasp interactive solar irradiance data center (lisird), 2021.  URL https://lasp.colorado.edu/lisird/.
  * Chamberlin et al. [2020] P. C. Chamberlin, F. G. Eparvier, V. Knoer, H. Leise, A. Pankratz, M. Snow, B. Templeman, E. M. B. Thiemann, D. L. Woodraska, and T. N. Woods.  The flare irradiance spectral model-version 2 (fism2).  _Space Weather_ , 18(12):e2020SW002588, 2020.  doi.org/10.1029/2020SW002588.  URL agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2020SW002588.  e2020SW002588 10.1029/2020SW002588.
  * Jacchia and Slowey [1966] L. G. Jacchia and J. Slowey.  The Shape and Location of the Diurnal Bulge in the Upper Atmosphere.  _SAO Special Report_ , 207, April 1966.
  * Vigil et al. [2021] Genevieve D Vigil, Amy Winebarger, Laurel Rachmeler, Nicolas Donders, P Subramania Athiray, Ken Kobayashi, and Charles Kankelborg.  Design for a portable calibration system for the full-sun uv rocket spectrometer instrument.  _Journal of Astronomical Telescopes, Instruments, and Systems_ , 7(3):035009, 2021.

##  Appendix A Details of Differential Path Length Formula

The following shows more details in the derivation of the differential path
length calculations. Refer to Figure 1 as needed. To solve for the
differential path length $\Delta L_{i}$, we first employ the Law of Sines.
Firstly, since $\sin\left(\pi-\theta_{0}\right)=\sin\theta_{0}$, then

$\displaystyle\frac{\sin\theta_{0}}{r_{0}+\Delta r_{0}}$
$\displaystyle=\frac{\sin\left(\theta_{L0}\right)}{r_{0}}\implies\theta_{L0}=\sin^{-1}\left(\frac{r_{0}}{r_{0}+\Delta
r_{0}}\sin\theta_{0}\right).$ (A1)

Then we use the 180° rule to solve for $\theta=\theta_{0}-\theta_{L0}$.
Finally, after some algebra, we use the Law of Sines again to solve for
$\Delta L_{0}$.

$\displaystyle\frac{\sin\theta_{0}}{r_{0}+\Delta
r_{0}}=\frac{\sin\theta}{\Delta L_{0}}\implies\Delta L_{0}$
$\displaystyle=\sqrt{\left(r_{0}+\Delta
r_{0}\right)^{2}-r_{0}^{2}\left(1-\mu_{0}^{2}\right)}-r_{0}\mu_{0},$ (A2)

where we allow the convention of $\mu_{0}=\cos\theta_{0}$. Extending to each
subsequent shell, we can see the pattern for $\Delta L_{i}$ as

$\displaystyle\Delta L_{0}$ $\displaystyle=\sqrt{\left(r_{0}+\Delta
r_{0}\right)^{2}-r_{0}^{2}\left(1-\mu_{0}^{2}\right)}-r_{0}\mu_{0}$ (A3)
$\displaystyle\Delta L_{1}$ $\displaystyle=\sqrt{\left(r_{0}+\Delta
r_{0}+\Delta r_{1}\right)^{2}-\left(r_{0}+\Delta
r_{0}\right)^{2}\left(1-\mu_{1}^{2}\right)}-\left(r_{0}+\Delta
r_{0}\right)\mu_{1}$ (A4) $\displaystyle\Delta L_{2}$
$\displaystyle=\sqrt{\left(r_{0}+\Delta r_{0}+\Delta r_{1}+\Delta
r_{2}\right)^{2}-\left(r_{0}+\Delta r_{0}+\Delta
r_{1}\right)^{2}\left(1-\mu_{2}^{2}\right)}-\left(r_{0}+\Delta r_{0}+\Delta
r_{1}\right)\mu_{2}$ (A5) $\displaystyle\vdots$ $\displaystyle\Delta L_{i}$
$\displaystyle=\sqrt{\left(r_{0}+\sum_{i^{\prime}=0}^{i}\Delta
r_{i^{\prime}}\right)^{2}-\left(r_{0}+\sum_{i^{\prime}=0}^{i-1}\Delta
r_{i^{\prime}}\right)^{2}\left(1-\mu_{i}^{2}\right)}-\left(r_{0}+\sum_{i^{\prime}=0}^{i-1}\Delta
r_{i^{\prime}}\right)\mu_{i},$ (A6) $\displaystyle\mu_{i}$
$\displaystyle=\sqrt{1-\left(\frac{r_{0}+\sum_{i^{\prime}=0}^{i-1}\Delta
r_{i^{\prime}}}{r_{0}+\sum_{i^{\prime}=0}^{i}\Delta
r_{i^{\prime}}}\right)^{2}\left(1-\mu_{i-1}^{2}\right)},\ i>0,$ (A7)
$\displaystyle\mu_{0}$ $\displaystyle=\cos\theta_{0}.$ (A8)

From Snell’s law of refraction, there exists an equivalence on either side of
the shell “wall,” namely

$\displaystyle n_{i}\sin\theta_{i}$
$\displaystyle=n_{i-1}\sin\theta_{i-1}\implies\mu_{i}^{2}=1-\left(\frac{n_{i-1}}{n_{i}}\right)^{2}\left(1-\mu_{i-1}^{2}\right),$
(A9)

where n is the index of refraction of O2 and $\theta$ is the angle to the
normal of the shell boundary. Upon substitution into Equation A7, we end up
with a correction factor onto $\mu_{i}$ given by the ratio of the index of
refraction at each subsequent shell (${n_{i-1}}/{n_{i}}$). The relationship
between density and refractive index can be expressed linearly if $\rho\leq 1$
g/cm-3 (see Liu and Daum [2008], Phillips [1920], and others). For O2, this
corresponds to a number density $\eta_{O2}\leq 10^{22}$ atoms/cm3). If the
density is low, it can be expected that the refractive index $n\rightarrow 1$.
From the Lorenz-Lorentz relation Buckingham and Graham [1974] we can calculate
this as

$\displaystyle\frac{n^{2}-1}{n^{2}+2}\rho^{-1}$
$\displaystyle=\frac{N_{A}\alpha}{3\epsilon_{0}}\quad\overset{n\to
1}{\Longrightarrow}\quad n=\frac{N\alpha_{0}}{2\epsilon_{0}}\rho+1.\implies
n_{i}=\frac{\eta_{i}}{\eta_{i-1}}\left(n_{i-1}-1\right)+1.$ (A10)

In Equation A10, $N_{A}$ and $\epsilon_{0}$ are constants. The polarizability
$\alpha$ is highly dependant upon density and wavelength. However, under
linearity approximation at low density, the dependency can be expressed as a
function of number density $\eta$. Taking the ratio of refractive indices of
subsequent shells removes those constants, and only a single value for the
refractive index at the instrument starting location $n_{0}$ is needed.

Using those same low-density assumptions, it is known through the Cauchy
formula [Cauchy, 1830] and Born and Wolf [2013] that the dispersion formula
for the refractive index as a function of wavelength can be expressed as

$\displaystyle n_{0}$ $\displaystyle=A\left(1+\frac{B}{\lambda^{2}}\right)+1,$
(A11)

where constants A and B are tabulated from laboratory measurements to be
$2.663\cdot 10^{-4}$ and $5.07\cdot 10^{3}$ nm2 respectively. A more precise
dispersion relation can be formulated using polarizability functions, but this
is not necessary for a first-order solution. Since we have shown that the
refractive index can be expressed as a function of wavelength, this implies
that the differential path length is as well.

##  Appendix B Simplifying Functions

There is a fair amount of simplification that occurs when we see how the
nested functions begin to cancel out terms. For the zenith angle, the
recursive functions terminate at $i=0$ as follows.

$\displaystyle\mu_{i}^{2}$
$\displaystyle=1-\left(\frac{r_{0}+\sum_{i^{\prime}=0}^{i-1}\Delta
r_{i^{\prime}}}{r_{0}+\sum_{i^{\prime}=0}^{i}\Delta
r_{i^{\prime}}}\right)^{2}\left(\frac{n_{i-1}}{n_{i}}\right)^{2}\left(1-\mu_{i-1}^{2}\right)$
(B1) $\displaystyle=1-\left(\frac{r_{0}+\sum_{i^{\prime}=0}^{i-1}\Delta
r_{i^{\prime}}}{r_{0}+\sum_{i^{\prime}=0}^{i}\Delta
r_{i^{\prime}}}\right)^{2}\left(\frac{n_{i-1}}{n_{i}}\right)^{2}\left(\frac{r_{0}+\sum_{i^{\prime}=0}^{i-2}\Delta
r_{i^{\prime}}}{r_{0}+\sum_{i^{\prime}=0}^{i-1}\Delta
r_{i^{\prime}}}\right)^{2}\left(\frac{n_{i-2}}{n_{i-1}}\right)^{2}\left(1-\mu_{i-2}^{2}\right)$
(B2) $\displaystyle=1-\left(\frac{r_{0}+\sum_{i^{\prime}=0}^{i-1}\Delta
r_{i^{\prime}}}{r_{0}+\sum_{i^{\prime}=0}^{i}\Delta
r_{i^{\prime}}}\right)^{2}\left(\frac{r_{0}+\sum_{i^{\prime}=0}^{i-2}\Delta
r_{i^{\prime}}}{r_{0}+\sum_{i^{\prime}=0}^{i-1}\Delta
r_{i^{\prime}}}\right)^{2}\dots$
$\displaystyle\quad\quad\quad\left(\frac{r_{0}+\Delta r_{0}}{r_{0}+\Delta
r_{0}+\Delta
r_{1}}\right)^{2}\left(\frac{n_{i-1}}{n_{i}}\right)^{2}\left(\frac{n_{i-2}}{n_{i-1}}\right)^{2}\dots$
$\displaystyle\quad\quad\quad\left(\frac{n_{0}}{n_{1}}\right)^{2}\left(1-\mu_{0}^{2}\right)$
(B3) $\displaystyle=1-\left(\frac{r_{0}+\Delta
r_{0}}{r_{0}+\sum_{i^{\prime}=i}^{i}\Delta
r_{i^{\prime}}}\right)^{2}\left(\frac{n_{0}}{n_{i}}\right)^{2}\left(1-\mu_{0}^{2}\right).$
(B4)

Likewise, the formula for refractive index reduces to Equations 9. To obtain
the final expression for $\Delta L$, we invoke the uniform shell-spacing
$\Delta r_{i}=\Delta r$. Thus, summations $\sum_{i^{\prime}=0}^{i}\Delta
r_{i}$ become $(i+1)\Delta r$. This returns Equations 7 and 8.

##  Appendix C Conversion of Signal to Data Numbers

A spectral signal is often given in units of ergs/s/Å/sr/cm2. Table 1 shows
the conversion process to units Data Numbers per exposure per pixel. Firstly,
we convert the energy from ergs to photons using the photon energy. Then we
have to account for the spectral plate scale for each wavelength channel in
FURST, which is $\approx 10.4$ nm across a width of 2048 pixels. Since we are
observing the full disk of the sun ($\approx$32 arc-minutes), then its size in
our field of view will be roughly constant at the value shown.

The signal is then attenuated by the instrument optics and the physical size
of the CCD. Finally, the photon count rate is spread across the height of the
CCD, which if the optics are aligned properly should match the number of rows
of pixels. In folding the photon signal into the instrument response function,
we take into account the charge spreading of the signal across adjacent
pixels. This width $\sigma$ is given from the FURST design and is around a
pixel in size depending on the wavelength. We compute the line-spread function
using a standard Gaussian,
$I_{0}\cdot\exp(-(\lambda-\lambda_{0})/2\sigma^{2})$, where $\lambda_{0}$ is
taken as the HRTS signal mapped onto the defined $\lambda$ space.

The camera sends this signal along the readout channels after a set exposure
time, which for this testing we have set to 10 seconds. In order to count up
the total number of photons during each exposure, we use the flight profile to
fold in the hanging transmission factor, as calculated by Equation 1 and shown
in the Results in Figure 4. In either case, it is at this point that the
signal is saved as an integer number of photons per exposure per pixel.

Variable | Source | Values | Units
---|---|---|---
HRTS Signal | Raw Data | $\approx 414^{+11912}_{-307}$ | erg/s/Å/sr/cm2
Photon Energy | Eph = h $\cdot$ c $\cdot$ $\lambda^{-1}$ | $\approx 1.32\cdot 10^{-11}$ | erg/ph
Spectral Plate Scale | Optical Setup | $5.078\cdot 10^{-3}$ | Å/x-pixel
Solar Disk Area | Observation | $9.352\cdot 10^{-3}$ | sr
Effective Area | GA $\cdot$ Refl. $\cdot$ Trans. $\cdot$ QE | $\approx 4.176\cdot 10^{-5}$ | cm2
Number of Rows | Variable | 1024 | y-pixel
Exposure Time | Variable | 10 | s/exp
Atmospheric Absorption | Model | 0.0-1.0 | unit-less
Photon Signal | Converted | $\approx 2070^{+91130}_{-1740}$ | ph/exp/pixel
Photon Noise | Poisson | $\approx\pm 110.5\ (3\sigma)$ | ph/pixel
Electrons per Photon | E${}_{\text{ph}}$/E${}_{\text{eh}}$ | $\approx 2.26$ | e-/ph
Electron Noise | Readout, Fano, Shot, etc. | $\approx\pm 25\ (3\sigma)$ | e-/pixel
Gain | Variable | 1.0 | DN/e-/pixel
Gain Bias | Variable | 0 (will be 3000) | DN/pixel
Camera Signal | Converted | $\approx 5645^{+258129}_{-4811}$ | DN/exp/pixel
Table 1: A summary of the unit conversions and sources of noise in the signal.
Beyond the Photon Signal step, the values are interpreted as integers by the
camera and readout systems.

Since there can be inherent variability in the count rate, we add Poisson
noise. Likewise, the CCD converts the incoming photons to electrons based on
the electron hole-pair Energy in Silicon (Eeh $=5.847944\cdot 10^{-12}$
erg/e-). The CCD reads out the total electron count from each pixel in terms
of voltage, with readout noise, gain, and gain bias. A separate calibration
code interprets the voltage signal to the computer system as a Data Number
based on this gain and gain-bias value. It is not necessary to include a gain-
bias value for this testing, since this value doesn’t fluctuate and is only
present to prevent negative values in the computer processing system.

##  Appendix D Inversion of Density Model

If we have enough confidence in the absorption cross-section, we can in
practice recalculate the density model instead. Starting as before, we can
take the ratio of the apogee image with any lower-altitude image as a proxy of
the transmission factor. We then use a forward model to invert the problem and
produce a potential solution for the atmospheric density profile. To do this,
we set up an optimization routine to vary the parameters of an assumed
functional form of the density profile.

This profile is carried through the spherical-shells model and compared with
the simulated “measured” optical depth. The profile used here was formulated
by observation of the NRL model to be a combination of a linear and parabolic
line. Thus, using a conventional approach often used in setting up a broken
power law (insert citation), we define the functional form of the density to
roughly follow

$\displaystyle y$
$\displaystyle=\frac{f_{1}-f_{2}}{\pi}\arctan\left(k\left(x_{0}-x\right)\right)+\frac{f_{1}+f_{2}}{2},$
(D1)

where y is the density, x is the altitude, and
$f_{1}=a(x-x_{0})^{2}+b(x-x_{0})+c$ and $f_{2}=b(x-x_{0})+c$ are the standard
parabolic and linear functions, respectively. The coefficients a, b, c, and x0
are adjusted by the optimization routine. The value for k adjusts the level of
sharpness around the break-point between the two functions and is simply set
to a relatively high number such as 1000.

Figure 10: Comparing the density model as used to generate the simulated
signal (green line) and the recomputed values through the optimization routine
(red line).

An initial guess for the parameters is given based on a fitting of the NRL
model. We then apply those parameters as the initial guess for computing the
density within equation 1. This result is compared with the signal data, and
the optimization routing minimizes the difference between the two by adjusting
those initial guesses. As can be seen in Figure 10, the fit is not perfect but
it is close. From these simulated noisy images we can reproduce the density
curve from the NRL model between the observed altitude and 1000km with a Root
Mean Squared Error (RMSE) of about 16%.

We do not consider this routine robust enough, as the initial guess is based
on required prior knowledge and complete trust in the sparse absorption cross-
section data set. This step will continue to be improved upon though, and
combining observations from various shutters and timings may improve this
routine. Additionally, we may be able to use the solved absorption cross-
section from one image to inform this inversion of the data from another
image.


obeypunctuation=true]University of Cambridge, UK

obeypunctuation=true]University of Cambridge, UK

obeypunctuation=true]University of Cambridge, UK

obeypunctuation=true]lowRISC CIC, UK

# Muntjac – Open Source Multicore RV64 Linux-capable SoC

Xuan Guo<EMAIL_ADDRESS>[  ,  Daniel Bates<EMAIL_ADDRESS>[  ,  Robert Mullins<EMAIL_ADDRESS>[  and  Alex Bradbury
<EMAIL_ADDRESS>[

(2022)

††journalyear: 2022††copyright: ccby††conference: First Workshop on Open-
Source Computer Architecture Research; June 18, 2022; New York City††doi:

##  1\. Introduction

Muntjac is an open-source111Available at https://github.com/lowRISC/muntjac
under Apache-2.0 licence. collection of components which can be used to build
a multicore, Linux-capable system-on-chip. This includes a 64-bit RISC-V core,
a cache subsystem, and TileLink (SiFive, 2020) interconnect allowing cache-
coherent multicore configurations. Each component is easy to understand,
verify, and extend, with most being configurable enough to be useful across a
wide range of applications.

Figure 1. Example multi-core Muntjac system. All solid-border components are
provided by Muntjac.

Inspired by the lowRISC’s Ibex project (lowRISC, 2021), we prioritise
verification and documentation so users can make rapid progress and be
confident that Muntjac’s behaviour adheres to all appropriate standards. We
focus on having clean, well-tested designs with clear routes to further
customisation and improvements. Muntjac’s performance is competitive with
other open-source projects, and we aim for a power/performance/area compromise
that maximises the value of Muntjac as a baseline design for others to use.

We anticipate Muntjac being valuable in the following situations, among
others:

  * •

Education: the complete code for a high-quality modular processor is
available, with documentation, verification support and examples, and lends
itself to self-contained extensions.

  * •

Research: Muntjac can serve as a solid baseline for experimentation, removing
the need to build and test huge swathes of general-purpose infrastructure
needed to support a single novel component.

  * •

Industry: many real-world applications value “time to solution”. Muntjac is a
reliable, configurable, complete starting point, allowing engineers to focus
on the unique selling points of their target system.

Why use Muntjac? At the project’s inception, there were a number of similar
RISC-V projects available, but none had the right mix of features for us. We
believe we have now reached a unique point in terms of performance,
documentation, verification and ease of understanding. We have a consistent
coding style, all components have been designed together to complement each
other, and we are able to explain the design better to new contributors.
Accessibility is aided by a relatively small and well-structured codebase, and
usage of popular, well-supported hardware description languages
(SystemVerilog) and standards (RISC-V, TileLink) to maximise interoperability
with other tools and IP. Everything “just works” out-of-the-box, and can be
configured to suit new designs.

We aim to grow the community of Muntjac users, and together we will extend and
improve the range of components we offer.

##  2\. Core Overview

The Muntjac core is single-issue, in-order, scalar, and supports the RV64GC
instruction set with machine and supervisor ISA. Floating-point extensions (F
and D) are optional and can be configured with SystemVerilog parameters.

The Muntjac core is designed with modularity in mind. It aims to be easy to
understand, verify and extend. We aim to make components as loosely coupled as
possible: in general, we choose distributed stall signals over a global stall
signal or a stall-free design. Valid-ready signals are widely used both within
components and between components, so each component and each pipeline stage
can be mostly self-contained. Skid buffers are used when distributed stall
signals start to cause timing issues.

##  3\. Pipeline Overview

Figure 2. Muntjac frontend design

Muntjac’s frontend design is illustrated in Figure 2. There is a separate
PCGEN stage that generates the next PC, without peeking into the fetched
instruction bytes.

Muntjac supports the compressed instruction (C) extension for code size
reduction. We do not provide a parameter to turn if off; this eliminates the
misaligned instruction exception, which is the only exception that can
possibly be generated by jump instructions, simplifying the backend. The
complexity of compressed instructions and misaligned instructions is entirely
handled within the frontend; the instruction cache only needs to support
accesses aligned to word (32-bit) boundaries.

The branch target buffer (BTB), return address stack (RAS) and branch
predictor are used to predict the next PC. Currently, the branch predictor is
a simple bi-modal 2-bit saturating counter, but like other components, the
interface between frontend and branch predictor is well-defined and a
different implementation can be easily swapped in.

Figure 3. Muntjac backend design

The backend is depicted in Figure 3. The issue logic manages control, data and
structural hazards, and issues instructions into one of the functional units.

The arithmetic logic unit (ALU) and the branching unit complete in a single
cycle, while the latency of other functional units may be multiple and
variable number of cycles. To hide the latency of data cache accesses (usually
2 cycles when cache hits), there are two execution stages, EX1 and EX2. All
instructions progress from EX1 to EX2 before their result is written back.

CSR accesses and other system instruction executions are handled outside the
normal data flow; as they can affect global states they can only execute while
the pipeline is empty.

An optional FPU is available with Muntjac. Unlike many other implementations,
we opted not to use a recoded format internally, so that NaN-boxing behaviour
is more easily implemented and hard-to-test bugs related to recoding can be
avoided. We provide an option that enables FP registers (and FP load/store
instructions) but disables the FPU. This mode allows efficient emulation of a
RV64GC core without bearing the cost of the FPU: M-mode firmware can trap and
emulate the floating-point operations, but trapping is not needed when the
operating system (OS) merely saves or restores FP contexts without performing
actual computations.

##  4\. Cache Overview

Muntjac provides cache implementations with interfaces that comply with the
TileLink protocol. Both instruction and data caches communicate with the
pipeline via generic valid-ready interfaces, so they have to manage the
request replays themselves without help from the frontend/backend.

Muntjac’s data cache implements the TileLink Cached (TL-C) protocol for its
memory-facing interface and is therefore multi-core-capable. The cache is set-
associative, and for each way tags and data are stored in separate single-port
SRAMs. Tag SRAM, data SRAM and TLB lookup all take place in parallel when an
access request arrives from the backend. The cache has an AMOALU independent
from backend’s ALU for atomic operations. Forward progress is guaranteed by
locking cache lines for 16 cycles after a refill.

Muntjac’s instruction cache is similar to the data cache but with narrow
access and write support removed, and width reduced to 32-bit. It implements
the TileLink Uncached Heavy (TL-UH) protocol.

Muntjac provides a range of TileLink IPs allowing individual Muntjac cores to
be aggregated to a multi-core system. For simple designs, a broadcaster can be
used to ensure data L1s are coherent. A reference directory-based L2 cache
design is also provided. This cache uses TL-C for both the CPU-facing link and
memory-facing link, so despite its name, it could be cascaded and used as L3
caches as well. This cache design can also be used as individual banks of a
larger cache. Various adapter and bridge IPs are provided to convert TL-C to
TL-UH and from TL-UH to AXI.

##  5\. Utilisation and Performance

For utilisation and performance numbers, we use the default configuration
parameters. The L1 data and instruction cache are 16 KiB each, 4-way
associative. The L1 data and instruction TLBs are 32 entry each and 4-way
associative. The L2 cache is 64KiB in size and also 4-way associative.

When targeting Xilinx Kintex 7 with Synopsys Synplify, with FPU on, the core
uses 17420 LUTs and 6683 registers (13663 LUTs and 4538 registers for the
pipeline) at 89MHz; with FPU off, the core uses 9902 LUTs and 5266 registers
(6022 LUTs and 3098 registers for the pipeline alone) at 97MHz; with FP
registers but no FPUs, the core uses 9924 LUTs and 5281 registers (6111 LUTs
and 3121 registers for the pipeline) at 95MHz. Muntjac attempts to be both
FPGA and ASIC friendly and we did not employ FPGA-specific optimisations.

On GCC 9.2.0, with optimisation and tuning turned on, Muntjac achieves
Dhrystone score of 2.17 DMIPS/MHz and CoreMark score of 3.01/MHz.

##  6\. Verification

We have put considerable effort into verification given its importance for
encouraging adoption of the IP and growing its user base; we believe it can
empower users by simplifying the process of validating modifications and
extensions.

We validate the Muntjac core using a combination of systematic tests (riscv-
tests (RISC-V International, 2021)) and random test generation (riscv-dv
(Google, 2021)). We also regularly run ad-hoc system tests such as booting
Debian Linux and running parallel benchmarks such as PARSEC (Bienia, 2011). We
test our TileLink IP using an assertion suite, a set of functional
coverpoints, and a custom random traffic generator.

##  7\. Acknowledgement

This work was kindly supported by lowRISC CIC.

## References

  * (1)
  * Bienia (2011) Christian Bienia. 2011\.  _Benchmarking Modern Multiprocessors_.  Ph. D. Dissertation. Princeton University.
  * Google (2021) Google. 2021.  RISCV-DV: open-source instruction generator for RISC-V processor verification.  https://github.com/google/riscv-dv.  Accessed: 2021-11-08.
  * lowRISC (2021) lowRISC. 2021\.  Ibex RISC-V Core.  https://github.com/lowRISC/ibex.  Accessed: 2021-09-24.
  * RISC-V International (2021) RISC-V International. 2021\.  riscv-tests: unit tests for RISC-V processors.  https://github.com/riscv-software-src/riscv-tests.  Accessed: 2021-11-08.
  * SiFive (2020) SiFive. 2020.  TileLink Specification 1.8.1.  https://sifive.cdn.prismic.io/sifive/7bef6f5c-ed3a-4712-866a-1a2e0c6b7b13_tilelink_spec_1.8.1.pdf.  Accessed: 2021-09-24.