If you are viewing this notebook through notebookviewer online and wish to install and use AtmoBuilder then please clone the github repository to get started.
AtmoBuilder requires the following to installed and appropriately sourced (including all pre-requisite software):
- LSST Software Stack
- sims_photUtils
- sims_sed_library
- AstroMLAtmoBuilder's wavelength gridding is limited by the availabe MODTRAN data. As a result this requires the PhysicalParameters.py in sims_photUtils to be altered as follows:
Stack Default:
- minwavelen = 300
- maxwavelen = 1150.1
- wavelenstep = 0.1
Required by AtmoBuilder:
- minwavelen = 300
- maxwavelen = 1100
- wavelenstep = 0.5
In the future this step may not be necessary, keep an eye out on updates to the github repository!
Atmospheric extinction can be modeled by 5 individual components: molecular (Rayleigh) scattering, aerosol (Mie) scattering, molecular absorption by ozone (O$_3$), water vapor and combined O$_2$ trace species (see Section 5.2 in LSST document LSE-180):
$S^{atm}(\lambda) = \Pi_{k=1}^5 \, {\rm e}^{- t_k \, \tau_k^{std}(\lambda, X)} = {\rm e}^{-\sum_{k=1}^5 t_k \, \tau_k^{std}(\lambda, X)}$
Here $X$ is airmass, $t_k$ is the ratio of the column density (or optical depth) of the $k$-th component and its value for the standard atmosphere, and $\tau_k^{std}(\lambda, X)$ is the wavelength-dependent optical depth of the $k$-th component for the standard atmosphere (with $t_k=1$) viewed at airmass $X$.
Aerosol (Mie) scattering data is not provided with MODTRAN data and has been modeled with a power-law relationship:
$ \tau_{aerosol}(\lambda, X, \alpha) = 0.1 \, X \, \left({550 \, {\rm nm} \over \lambda} \right)^\alpha $
The power-law exponent ($\alpha$) is chosen to be the 6th free parameter when modeling atmospheres. In summation, we can define a single parameter array of $t_k$ values for each of the 5 extincton components and 1 for the power-law exponent:
$ [t_{H_2O}, t_{O_2}, t_{O_3}, t_{Rayleigh}, t_{Aerosol}, t_{\alpha}] $
The standard atmosphere is defined by the following parameter array and airmass:
- P: [1.0,1.0,1.0,1.0,1.0,1.7]
- X: 1.2
In order to create an atmosphere for use with AtmoBuilder, a parameter array and airmass like the one above must be given. Details on the specific procedure can be found in the following sections.
Importing the AtmoBuilder class:
In [1]:
from AtmoBuilder import AtmoBuilder
% matplotlib inline
Initializing the class:
In [2]:
ab = AtmoBuilder()
When initialized AtmoBuilder will output the code above to notify you that it has found and read MODTRAN data, and that it has found filter and hardware data from the LSST software stack.
Retrieving MODTRAN data:
In [3]:
# List of read airmasses
ab.airmasses
Out[3]:
In [4]:
# Standard transmission profile keyed to airmass and component
ab.transDict[1.2]['Rayleigh']
Out[4]:
Plotting the filter and hardware response curves (hardware includes filters):
In [5]:
ab.filterPlot()
In [6]:
ab.hardwarePlot()
AtmoBuilder utilizes an independent atmo object which keeps track of several useful elements. Let's start by creating two atmospheres:
In [7]:
atmo = ab.buildAtmo([1.0,1.0,1.0,2.2,1.4,1.7],1.4)
atmo_std = ab.buildAtmo([1.0,1.0,1.0,1.0,1.0,1.7],1.2)
The line of code above has two important keyword parameters:
- P: parameter array containing 6 values, one for each of the atmospheric extinction components
- X: airmassRetrieving atmo attributes:
In [8]:
# Parameters and airmass
atmo.P, atmo.X
Out[8]:
In [9]:
# Total transmission [array]
atmo.sb
Out[9]:
In [10]:
# Individual component transmission profiles [component-keyed dictionary]
atmo.sbDict
Out[10]:
AtmoBuilder has a lot of various plotting functions including one for atmospheres. The following function will take an atmo object and plot its transmission profile:
In [11]:
ab.transPlot(atmo)
AtmoBuilder can also plot the individual components:
In [12]:
ab.transPlot(atmo, includeComponents=True)
Want to compare the atmosphere to another atmosphere?
In [13]:
ab.transPlot(atmo, atmo2=atmo_std)
Of course, keeping track of a standard atmosphere like the one we created can become laborious. To avoid this, internally in AtmoBuilder there are class constants which define the parameters and airmass for a standard atmosphere. For all major functionality you need not worry about keeping track of the the standard atmosphere. For example:
In [14]:
ab.transPlot(atmo,includeStdAtmo=True)
In order to run regressions and to use many of the basic operations and functions in AtmoBuilder we need SED data to be read in. This is not done when the class is initialized because it can take a fair amount of time depending on how many SEDs you want to read in.
By default AtmoBuilder uses six SED types:
- Kurucz Main Sequence Stars (read from LSST's SED library) ['mss']
- White Dwarf Stars (read from LSST's SED library) ['wds']
- Galaxies (generated at different redshifts from LSST's SED library) ['gals']
- MLT Dwarf Stars (read from LSST's SED library) ['mlts']
- Quasars (generated at different redshifts from AtmoBuilder repository) ['qsos']
- Supernovae (generated at different redshifts from AtmoBuilder repository) ['sns']Reading in Kurucz MS stars:
In [15]:
ab.readMSs()
You can call similar functions for each of the six SED types to read them. If you want to use all of the available SED types (or a large subset) you can call the following (this will take a little bit of time):
In [16]:
ab.readAll()
If you want to run some simple operations it may be helpful to retrieve the SED data you read in. With AtmoBuilder it is remarkably easy to do so as it is all stored as class attributes. Lets retrieve some quasar data:
In [17]:
qso = ab.qsos
qso_z = ab.qsoRedshifts
We can see what the SED attributes are set when SED data is read in two ways:
In [18]:
ab?
... or we can look at each of the individual SED read functions:
In [19]:
ab.readMSs?
Once you have some SED data read in you can use them to calculate magnitudes, delta magnitudes, and so on...
Lets combine hardware data and an atmosphere object to get the combined throughput:
In [20]:
throughput = ab.combineThroughputs(atmo)
throughput_std = ab.combineThroughputs(atmo_std)
We can then also plot the combined throughput:
In [21]:
ab.throughputPlot(throughput)
Like transPlot, we can plot two throughputs or include the standard throughput:
In [22]:
ab.throughputPlot(throughput, bpDict2=throughput_std)
ab.throughputPlot(throughput, includeStdAtmo=True)
A characteristic function we care about is the normalized bandpass response function for a given throughput. The phi function can be retrieved from a filter-keyed bandpass dictionary as follows (example for u filter):
In [23]:
throughput['u'].phi
Out[23]:
Because these throughputs are bandpass objects, you can also access the wavelength range and the overall throughput:
In [24]:
throughput['u'].sb
Out[24]:
In [25]:
throughput['u'].wavelen
Out[25]:
AtmoBuilder can also plot phi functions:
In [26]:
ab.phiPlot(throughput)
Of course it might be handy to compare one phi function to another:
In [27]:
ab.phiPlot(throughput,throughput_std)
AtmoBuilder can also plot the change in normalized bandpass response function:
In [28]:
ab.dphiPlot(throughput, throughput_std)
We might want to compare one dphi to another:
In [29]:
throughput_temp = ab.combineThroughputs(ab.buildAtmo([2.3,2.1,1.2,0.9,1.3,1.8],1.4))
ab.dphiPlot(throughput, throughput_std, throughput_temp)
Sometimes its hard to really see the difference between two phi functions as is often the case when plotting regression data. We can use the ddphiPlot to plot the difference between to dphis in the form of a subtraction between them:
In [30]:
ab.ddphiPlot(throughput,throughput_temp,throughput_std)
Using the throughput we can calculate the magnitudes of stars we read in (defaults to using Kurucz MS SEDs)
In [31]:
mags = ab.mags(throughput)
mags_std = ab.mags(throughput_std)
Or we can calculate the magnitudes and change in magnitudes of other SEDs:
In [32]:
qso_mags = ab.mags(throughput, seds=qso, sedkeylist=qso_z)
qso_mags_std = ab.mags(throughput_std, seds=qso, sedkeylist=qso_z)
In [33]:
dmags = ab.dmags(mags, mags_std)
qso_dmags = ab.dmags(qso_mags, qso_mags_std)
Both mags and dmags are returned as filter-keyed dictionaries of magnitudes:
In [34]:
mags['u']
dmags['u']
Out[34]:
Suppose we want to see the change in magnitudes between two atmospheres for a single SED type, we can use the following function:
In [35]:
ab.dmagPlot(throughput,throughput_std,'qsos')
Or another SED type example:
In [36]:
ab.dmagPlot(throughput,throughput_std,'gals')
Regressions (currently limited to two varying parameters) are the most integral part of AtmoBuilder and are run using a single function.
In [37]:
ab.computeAtmoFit?
Making use of the atmospheres we created at the beginning of this tutorial, lets run some regressions. Note that we only changed two parameters relative to the standard atmosphere when we created the atmosphere. These were the parameters for the Rayleigh and Aerosol components of atmospheric extinction, so let us regress those two components. (I have added the pickleString keyword argument which saves the pickle dump files and plots with an additional string component in its file name)
In [38]:
ab.computeAtmoFit('Rayleigh', 'Aerosol', atmo, pickleString='Tutorial')
Sometimes the contour plots do not adequately show the log-likelihood matrix, we can directly plot the logL matrix instead of the regression contours using 'useLogL':
In [39]:
ab.computeAtmoFit('Rayleigh', 'Aerosol', atmo, useLogL=True, pickleString='Tutorial')
It is also possible to plot both the contours and the log-likelihood matrix on the same plot:
In [40]:
ab.computeAtmoFit('Rayleigh','Aerosol',atmo, plotBoth=True, pickleString='Tutorial')
A colorbar for the log-likelihood column might be handy and can be included using a keyword argument. Unfortunately, it does resize the middle column slightly.
In [41]:
ab.computeAtmoFit('Rayleigh', 'Aerosol', atmo, useLogL=True, includeColorBar=True, pickleString='Tutorial')
By default regressions will run using the Kurucz main sequence star SED data. However, we can also specify a different regression SED: (I have also set generate Dphi to False to condense the notebook a little)
In [42]:
ab.computeAtmoFit('Rayleigh', 'Aerosol', atmo, regressionSed='qsos', generateDphi=False, pickleString='Tutorial')
Like the regression SEDs, we can also set which comparison SEDs to use. By default, computeAtmoFit will always use the remaining SEDs (all of them except the original regression SED).
In [43]:
ab.computeAtmoFit('Rayleigh', 'Aerosol', atmo, regressionSed='qsos', comparisonSeds=['gals','sns'], generateDphi=False, pickleString='Tutorial')
The default percent error value for dmags in mmags is set to be 0.5%, however we can also set it to other values. For example, we can set it to be 5%:
In [44]:
ab.computeAtmoFit('Rayleigh', 'Aerosol', atmo, err=0.05, generateDphi=False, pickleString='Tutorial')
One additional feature of interest is the deltaGrey keyword parameter. It is essentially an additional extinction factor implemented to add some offset representative of cloudy and potentially grey skies.
It works in the following way:
- deltaGrey = 0.0 : default, no effect
- deltaGrey >= 0.0 : subtract deltaGrey value from dmags (dG should be in mmags, for example: 20)
- deltaGrey < 0.0 : subtract median dmag value from dmags
In [45]:
ab.computeAtmoFit('Rayleigh', 'Aerosol', atmo, deltaGrey=20, generateDphi=False, pickleString='Tutorial')
In [46]:
ab.computeAtmoFit('Rayleigh', 'Aerosol', atmo, deltaGrey=-1, generateDphi=False, pickleString='Tutorial')
In [ ]: