desihub / specsim

Quick simulations of spectrograph response
2 stars 9 forks source link

obsflux artifacts #19

Closed sbailey closed 9 years ago

sbailey commented 9 years ago

obsflux has non-physical artifacts that appear to come from edge effects at the edges of each camera coverage:

import specsim
import numpy as np

atmosphere = specsim.atmosphere.Atmosphere(
    skyConditions='dark', basePath=os.environ['DESIMODEL'])
qsim = specsim.quick.Quick(
    atmosphere=atmosphere, basePath=os.environ['DESIMODEL'])

wavemin, wavemax, wavestep = 5600, 6000, 1.0
qsim.setWavelengthGrid(wavemin, wavemax-2, wavestep)
wave = np.arange(wavemin, wavemax, wavestep)
flux = np.ones(len(wave))
inspec = specsim.spectrum.SpectralFluxDensity(wave, flux)
results = qsim.simulate(sourceType='elg', sourceSpectrum=inspec, 
    airmass = 1.25, expTime = 1000, downsampling=1)

axvline(qsim.cameras[0].wavelengthRange[-1], color='r', ls=':')
axvline(qsim.cameras[1].wavelengthRange[0], color='r', ls=':')
plot(results.wave, results.obsflux)

specsim_artifacts

Admittedly, the definition of "obsflux" is a bit ill-defined in the overlap region between two cameras since it depends upon the resolution of 2 different cameras. But a flat input spectrum should still result in a flat output spectrum without divots.

sbailey commented 9 years ago

The problem arises from specsim.quick around lines 353-354:

for camera in self.cameras:
    # [...]
    # Estimate this camera's contribution to the coadded observed source flux.
    self.observedFlux += camera.throughput*camera.sparseKernel.dot(self.sourceFlux)
    throughputTotal += camera.throughput

# Approximate the observed source flux spectrum as the throughput-weighted sum of
# the flux smoothed by each camera's resolution.
thruMask = throughputTotal > 0
self.observedFlux[thruMask] /= throughputTotal[thruMask]

At the edges of the cameras, flux gets lost by being convolved into the zero-throughput region, from where it is never recovered by the throughput normalization. This can be "fixed" using a resolution convolved throughput:

    self.observedFlux += camera.sparseKernel.dot(camera.throughput*self.sourceFlux)
    throughputTotal += camera.sparseKernel.dot(camera.throughput)

though I'm not sure that actually represents the physical order in which resolution and throughput actually happen to photons. For consideration.

julienguy commented 9 years ago

Fixed in branch fixedges In QuickCamera, extend sparseKernel computing at +- nhalf of wavelength range with positive throughput to fix obsflux edge artifact. Note however that this is not ideal : the throughput should be applied before the resolution convolution, because it is a function of the true wavelength.

Artefacts on the edges of the cameras positive throughput regions have disappeared. There are still edges effect on the edge of the input wavelength range, because there is nothing we can do about it (we don't know what is the flux outside this range). Flux is only conserved when integrated over a full PSF width.

image

dkirkby commented 9 years ago

Closed with #20