from __future__ import absolute_import, division, print_function
import sys, os
import numpy as np
from astropy.table import Table
import astropy.units as u
import specsim.simulator
from desispec.frame import Frame
import desispec.io
from desispec.resolution import Resolution
import desisim.io
import desisim.simexp
from desisim.util import dateobs2night
import desisim.specsim
#-------------------------------------------------------------------------
def parse(options=None):
import argparse
parser = argparse.ArgumentParser(usage = "{prog} [options]")
parser.add_argument("--simspec", type=str, help="input simspec file")
parser.add_argument("--outdir", type=str, help="output directory")
parser.add_argument("--firstspec", type=int, default=0,
help="first spectrum to simulate")
parser.add_argument("--nspec", type=int, default=5000,
help="number of spectra to simulate")
parser.add_argument("--cframe", action="store_true",
help="directly write cframe")
parser.add_argument("--dwave", type=float, default=0.8, help="output wavelength step, in Angstrom")
if options is None:
args = parser.parse_args()
else:
args = parser.parse_args(options)
return args
[docs]def main(args=None):
'''
Converts simspec -> frame files; see fastframe --help for usage options
'''
#- TODO: use desiutil.log
if isinstance(args, (list, tuple, type(None))):
args = parse(args)
print('Reading files')
simspec = desisim.io.read_simspec(args.simspec, readphot=False)
if simspec.flavor == 'arc':
print('arc exposure; no frames to output')
return
fibermap = simspec.fibermap
obs = simspec.obsconditions
night = simspec.header['NIGHT']
expid = simspec.header['EXPID']
firstspec = args.firstspec
nspec = min(args.nspec, len(fibermap)-firstspec)
print('Simulating spectra {}-{}'.format(firstspec, firstspec+nspec))
wave = simspec.wave
flux = simspec.flux
ii = slice(firstspec, firstspec+nspec)
if simspec.flavor == 'science':
sim = desisim.simexp.simulate_spectra(wave, flux[ii],
fibermap=fibermap[ii], obsconditions=obs, dwave_out=args.dwave,
psfconvolve=True)
elif simspec.flavor in ['arc', 'flat', 'calib']:
x = fibermap['FIBERASSIGN_X']
y = fibermap['FIBERASSIGN_Y']
fiber_area = desisim.simexp.fiber_area_arcsec2(x, y)
surface_brightness = (flux.T / fiber_area).T
config = desisim.simexp._specsim_config_for_wave(wave, dwave_out=args.dwave)
# sim = specsim.simulator.Simulator(config, num_fibers=nspec)
sim = desisim.specsim.get_simulator(config, num_fibers=nspec,
camera_output=True)
sim.observation.exposure_time = simspec.header['EXPTIME'] * u.s
sbunit = 1e-17 * u.erg / (u.Angstrom * u.s * u.cm ** 2 * u.arcsec ** 2)
xy = np.vstack([x, y]).T * u.mm
sim.simulate(calibration_surface_brightness=surface_brightness[ii]*sbunit,
focal_positions=xy[ii])
else:
raise ValueError('Unknown simspec flavor {}'.format(simspec.flavor))
sim.generate_random_noise()
for i, results in enumerate(sim.camera_output):
results = sim.camera_output[i]
wave = results['wavelength']
scale=1e17
if args.cframe :
phot = scale*(results['observed_flux'] + results['random_noise_electrons']*results['flux_calibration']).T
ivar = 1./scale**2*results['flux_inverse_variance'].T
else :
phot = (results['num_source_electrons'] + \
results['num_sky_electrons'] + \
results['num_dark_electrons'] + \
results['random_noise_electrons']).T
ivar = 1.0 / results['variance_electrons'].T
R = Resolution(sim.instrument.cameras[i].get_output_resolution_matrix())
Rdata = np.tile(R.data.T, nspec).T.reshape(
nspec, R.data.shape[0], R.data.shape[1])
assert np.all(Rdata[0] == R.data)
assert phot.shape == (nspec, len(wave))
for spectro in range(10):
imin = max(firstspec, spectro*500) - firstspec
imax = min(firstspec+nspec, (spectro+1)*500) - firstspec
if imax <= imin:
continue
xphot = phot[imin:imax]
xivar = ivar[imin:imax]
xfibermap = fibermap[ii][imin:imax]
camera = '{}{}'.format(sim.camera_names[i], spectro)
meta = simspec.header.copy()
meta['CAMERA'] = camera
if args.cframe :
units = '1e-17 erg/(s cm2 Angstrom)'
else :
#
# We want to save electrons per angstrom and not electrons per bin
# to be consistent with the extraction code (specter.extract.ex2d).
# And to be FITS-compliant, we call electrons "counts".
#
units = 'count/Angstrom'
dwave=np.gradient(wave)
xphot /= dwave
xivar *= dwave**2
meta['BUNIT']=units
meta['DETECTOR'] = 'SIM'
if camera[0] == 'b':
meta['CCDSIZE'] = '4162,4232'
else:
meta['CCDSIZE'] = '4194,4256'
readnoise = sim.instrument.cameras[i].read_noise.value
meta['OBSRDNA'] = readnoise
meta['OBSRDNB'] = readnoise
meta['OBSRDNC'] = readnoise
meta['OBSRDND'] = readnoise
frame = Frame(wave, xphot, xivar, resolution_data=Rdata[0:imax-imin],
spectrograph=spectro, fibermap=xfibermap, meta=meta)
if args.cframe :
outfile = desispec.io.findfile('cframe', night, expid, camera,
outdir=args.outdir)
else :
outfile = desispec.io.findfile('frame', night, expid, camera,
outdir=args.outdir)
print('writing {}'.format(outfile))
desispec.io.write_frame(outfile, frame, units=units)