Skip to content

Commit

Permalink
Merge pull request #375 from lsst/u/yoachim/sed_twilight
Browse files Browse the repository at this point in the history 
force wavelenth arrays to be idenical for twilight solar spectra
  • Loading branch information
@yoachim
yoachim committed Oct 19, 2023
2 parents 7d18279 + f3ac82a commit c3affd5
Showing 1 changed file with 63 additions and 35 deletions.
98 changes: 63 additions & 35 deletions rubin_sim/skybrightness/interp_components.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,17 +78,21 @@ def load_spec_files(filenames, mags=False):
The ESO npz files contain the following arrays:
* filter_wave: The central wavelengths of the pre-computed magnitudes
* wave: wavelengths for the spectra
* spec: array of spectra and magnitudes along with the relevant variable inputs. For example,
airglow has dtype = [('airmass', '<f8'), ('solar_flux', '<f8'), ('spectra', '<f8', (17001,)),
('mags', '<f8', (6,)]
For each unique airmass and solar_flux value, there is a 17001 elements spectra and 6 magnitudes.
* spec: array of spectra and magnitudes along with the relevant
variable inputs. For example, airglow has dtype =
[('airmass', '<f8'), ('solar_flux', '<f8'), ('spectra',
'<f8', (17001,)),
('mags', '<f8', (6,)]
For each unique airmass and solar_flux value, there is a 17001 elements
spectra and 6 magnitudes.
"""

if len(filenames) == 1:
temp = np.load(filenames[0])
wave = temp["wave"].copy()
# Note old camelCase here. Might need to update if sav files regenerated
# Note old camelCase here. Might need to update if sav
# files regenerated
filter_wave = temp["filterWave"].copy()
if mags:
# don't copy the spectra to save memory space
Expand Down Expand Up @@ -171,7 +175,8 @@ def __init__(self, comp_name=None, sorted_order=["airmass", "nightTimes"], mags=
else:
self.spec_size = 0

# What order are the dimesions sorted by (from how the .npz was packaged)
# What order are the dimesions sorted by (from how
# the .npz was packaged)
self.sorted_order = sorted_order
self.dim_dict = {}
self.dim_sizes = {}
Expand All @@ -190,8 +195,8 @@ def __call__(self, intep_points, filter_names=["u", "g", "r", "i", "z", "y"]):

def indx_and_weights(self, points, grid):
"""
for given 1-D points, find the grid points on either side and return the weights
assume grid is sorted
for given 1-D points, find the grid points on
either side and return the weights assume grid is sorted
"""

order = np.argsort(points)
Expand All @@ -218,8 +223,9 @@ def indx_and_weights(self, points, grid):

def _weighting(self, interp_points, values):
"""
given a list/array of airmass values, return a dict with the interpolated
spectrum at each airmass and the wavelength array.
given a list/array of airmass values, return a
dict with the interpolated spectrum at each airmass
and the wavelength array.
Input interp_points should be sorted
"""
Expand All @@ -235,7 +241,8 @@ def _weighting(self, interp_points, values):

nextra = 1

# XXX--should I use the log spectra? Make a check and switch back and forth?
# XXX--should I use the log spectra?
# Make a check and switch back and forth?
results[in_range] = (
w_r[:, np.newaxis] * values[indx_r * nextra] + w_l[:, np.newaxis] * values[indx_l * nextra]
)
Expand Down Expand Up @@ -287,7 +294,8 @@ def __init__(self, comp_name="UpperAtm", mags=False):

class MergedSpec(BaseSingleInterp):
"""
Interpolate the spectra caused by the sum of the scattered starlight, airglow, upper and lower atmosphere.
Interpolate the spectra caused by the sum of the scattered
starlight, airglow, upper and lower atmosphere.
"""

def __init__(self, comp_name="MergedSpec", mags=False):
Expand Down Expand Up @@ -331,13 +339,20 @@ def _weighting(self, interp_points, values):

class TwilightInterp:
def __init__(self, mags=False, dark_sky_mags=None, fit_results=None):
"""
Read the Solar spectrum into a handy object and compute mags in different filters
mags: If true, only return the LSST filter magnitudes, otherwise return the full spectrum
dark_sky_mags = dict of the zenith dark sky values to be assumed. The twilight fits are
done relative to the dark sky level.
fit_results = dict of twilight parameters based on twilight_func. Keys should be filter names.
"""Read the Solar spectrum into a handy object and
compute mags in different filters
Parameters
----------
mags : `bool`
If True, only return the LSST filter magnitudes,
otherwise return the full spectrum
dark_sky_mags : dict
Dict of the zenith dark sky values to be assumed.
The twilight fits are done relative to the dark sky level.
fit_results : dict
Dict of twilight parameters based on twilight_func.
Keys should be filter names.
"""

if dark_sky_mags is None:
Expand Down Expand Up @@ -381,8 +396,11 @@ def __init__(self, mags=False, dark_sky_mags=None, fit_results=None):
bp_temp.resample_bandpass(
wavelen_min=self.solar_spec.wavelen.min(),
wavelen_max=self.solar_spec.wavelen.max(),
wavelen_step=0.2,
wavelen_step=self.solar_spec.wavelen[1] - self.solar_spec.wavelen[0],
)
# Force wavelengths to be identical so
# it doesn't try to resample again later
bp_temp.wavelen = self.solar_spec.wavelen
canon_filters[filter_name] = bp_temp

# Tack on the LSST filters
Expand All @@ -399,18 +417,22 @@ def __init__(self, mags=False, dark_sky_mags=None, fit_results=None):
self.filter_names.append(key)

# MAGIC NUMBERS from fitting the all-sky camera:
# Code to generate values in sims_skybrightness/examples/fitTwiSlopesSimul.py
# Which in turn uses twilight maps from sims_skybrightness/examples/buildTwilMaps.py
# Code to generate values in
# sims_skybrightness/examples/fitTwiSlopesSimul.py
# Which in turn uses twilight maps from
# sims_skybrightness/examples/buildTwilMaps.py
# values are of the form:
# 0: ratio of f^z_12 to f_dark^z
# 1: slope of curve wrt sun alt
# 2: airmass term (10^(arg[2]*(X-1)))
# 3: azimuth term.
# 4: zenith dark sky flux (erg/s/cm^2)

# For z and y, just assuming the shape parameter fits are similar to the other bands.
# For z and y, just assuming the shape parameter
# fits are similar to the other bands.
# Looks like the diode is not sensitive enough to detect faint sky.
# Using the Patat et al 2006 I-band values for z and modeified a little for y as a temp fix.
# Using the Patat et al 2006 I-band values for z and
# modeified a little for y as a temp fix.
if fit_results is None:
self.fit_results = {
"B": [
Expand Down Expand Up @@ -438,7 +460,8 @@ def __init__(self, mags=False, dark_sky_mags=None, fit_results=None):
"y": [2.0, 24.08, 0.3, 0.3, -666],
}

# XXX-completely arbitrary fudge factor to make things brighter in the blue
# XXX-completely arbitrary fudge factor to make things
# brighter in the blue
# Just copy the blue and say it's brighter.
self.fit_results["u"] = [
16.0,
Expand Down Expand Up @@ -471,10 +494,11 @@ def __init__(self, mags=False, dark_sky_mags=None, fit_results=None):

self.solar_wave = self.solar_spec.wavelen
self.solar_flux = self.solar_spec.flambda
# This one isn't as bad as the model grids, maybe we could get away with computing the magnitudes
# in the __call__ each time.
# This one isn't as bad as the model grids, maybe we could get
# away with computing the magnitudes in the __call__ each time.
if mags:
# Load up the LSST filters and convert the solarSpec.flabda and solarSpec.wavelen to fluxes
# Load up the LSST filters and convert the
# solarSpec.flabda and solarSpec.wavelen to fluxes
through_path = through_path = os.path.join(get_data_dir(), "throughputs", "baseline")
self.lsst_filter_names = ["u", "g", "r", "i", "z", "y"]
self.lsst_equations = np.zeros(
Expand All @@ -495,7 +519,8 @@ def __init__(self, mags=False, dark_sky_mags=None, fit_results=None):
temp_b = Bandpass()
temp_b.set_bandpass(bp["wave"], bp["trans"])
self.lsst_eff_wave.append(temp_b.calc_eff_wavelen()[0])
# Loop through the parameters and interpolate to new eff wavelengths
# Loop through the parameters and interpolate to new
# eff wavelengths
for i in np.arange(self.lsst_equations[0, :].size):
interp = InterpolatedUnivariateSpline(self.eff_wave, fits[:, i])
self.lsst_equations[:, i] = interp(self.lsst_eff_wave)
Expand Down Expand Up @@ -540,8 +565,9 @@ def interp_mag(
filter_names=["u", "g", "r", "i", "z", "y"],
):
"""
Originally fit the twilight with a cutoff of sun altitude of -11 degrees. I think it can be safely
extrapolated farther, but be warned you may be entering a regime where it breaks down.
Originally fit the twilight with a cutoff of sun altitude of
-11 degrees. I think it can be safely extrapolated farther,
but be warned you may be entering a regime where it breaks down.
"""
npts = len(filter_names)
result = np.zeros((np.size(interp_points), npts), dtype=float)
Expand Down Expand Up @@ -587,7 +613,8 @@ def interp_spec(self, interp_points, max_am=3.0, limits=(np.radians(15.0), np.ra
& (interp_points["airmass"] >= 1.0)
)[0]

# Compute the expected flux in each of the filters that we have fits for
# Compute the expected flux in each of the filters that
# we have fits for
fluxes = []
for filter_name in self.filter_names:
out_of_range = np.where(interp_points["sunAlt"] > np.max(limits))[0]
Expand Down Expand Up @@ -679,8 +706,9 @@ def _weighting(self, interp_points, values):

class ZodiacalInterp(BaseSingleInterp):
"""
Interpolate the zodiacal light based on the airmass and the healpix ID where
the healpixels are in ecliptic coordinates, with the sun at ecliptic longitude zero
Interpolate the zodiacal light based on the airmass
and the healpix ID where the healpixels are in ecliptic
coordinates, with the sun at ecliptic longitude zero
"""

def __init__(self, comp_name="Zodiacal", sorted_order=["airmass", "hpid"], mags=False):
Expand Down

0 comments on commit c3affd5

Please sign in to comment.