diff --git a/bin/pixsim-desi b/bin/pixsim-desi index 0163670e5..868376c2f 100755 --- a/bin/pixsim-desi +++ b/bin/pixsim-desi @@ -42,12 +42,14 @@ parser = optparse.OptionParser( ) parser.add_option("--newexp", action="store_true", help="Create new exposure") +parser.add_option("--tileid", type=int, help="tile id") parser.add_option("--expid", type=int, help="exposure id") parser.add_option("--night", type=str, help="YEARMMDD") +parser.add_option("--airmass", type=float, help="Airmass [%default]", default=1.0) ### parser.add_option("--simfile", type=str, help="input sim file") parser.add_option("--camera", type=str, help="camera (e.g. b0, r5, z9)") parser.add_option("--verbose", action="store_true", help="print more stuff") -parser.add_option("--randseed", type=int, help="random number seed [%default]", default=0) +parser.add_option("--randseed", type=int, help="random number seed") parser.add_option("--nspec", type=int, help="Number of spectra to simulate [%default]", default=500) parser.add_option("--ncpu", type=int, help="Number of cpu cores to use [%default]", default=multiprocessing.cpu_count() / 2) @@ -72,12 +74,20 @@ if opts.camera is None: else: opts.camera = opts.camera.split(',') +#- Initialize tileid +if opts.tileid is None: + opts.tileid = obs.get_next_tileid() + #- Initialize random seeds +if opts.randseed is None: + if opts.tileid >= 0: + opts.randseed = opts.tileid + random.seed(opts.randseed) np.random.seed(opts.randseed) if opts.newexp: - obs.new_exposure(opts.nspec, expid=opts.expid) + obs.new_exposure(opts.nspec, expid=opts.expid, tileid=opts.tileid, airmass=opts.airmass) sys.exit(0) for camera in opts.camera: diff --git a/py/desisim/io.py b/py/desisim/io.py index a79bb84a9..13a1ae45f 100644 --- a/py/desisim/io.py +++ b/py/desisim/io.py @@ -93,7 +93,7 @@ def write_simspec(meta, truth, expid, night): Write $DESI_SPECTRO_SIM/$PIXPROD/{night}/simspec-{expid}.fits Args: - meta : metadata table to write to "META" HDU + meta : metadata table to write to "METADATA" HDU truth : dictionary with keys: FLUX - 2D array [nspec, nwave] in erg/s/cm2/A WAVE - 1D array of vacuum wavelengths [Angstroms] @@ -136,7 +136,7 @@ def write_simspec(meta, truth, expid, night): O2FLUX = 'erg/s/cm2', ) - write_bintable(outfile, meta, header=None, extname="META", + write_bintable(outfile, meta, header=None, extname="METADATA", comments=comments, units=units) #- Write object photon and sky photons for each channel @@ -160,7 +160,15 @@ def write_simspec(meta, truth, expid, night): return outfile - +#------------------------------------------------------------------------- +#- Parse header to make wavelength array +def load_wavelength(filename, extname): + hdr = fits.getheader(filename, extname) + wave = hdr['CRVAL1'] + np.arange(hdr['NAXIS1'])*hdr['CDELT1'] + if hdr['LOGLAM'] == 1: + wave = 10**wave + return wave + #------------------------------------------------------------------------- #- desimodel #- These should probably move to desimodel itself, diff --git a/py/desisim/obs.py b/py/desisim/obs.py index bde40089d..80b36a3b7 100644 --- a/py/desisim/obs.py +++ b/py/desisim/obs.py @@ -62,13 +62,16 @@ def _dict2ndarray(data, columns=None): # # return result -def new_exposure(nspec=5000, expid=None): +def new_exposure(nspec=5000, expid=None, tileid=None, airmass=1.0): """ Create a new exposure and output input simulation files. Does not generate pixel-level simulations or noisy spectra. Args: nspec (optional): integer number of spectra to simulate + expid (optional): positive integer exposure ID + tileid (optional): tile ID + airmass (optional): airmass, default 1.0 Writes: $DESI_SPECTRO_SIM/$PIXPROD/{night}/fibermap-{expid}.fits @@ -80,10 +83,12 @@ def new_exposure(nspec=5000, expid=None): """ if expid is None: expid = get_next_expid() + + if tileid is None: + tileid = get_next_tileid() dateobs = time.gmtime() night = get_night(utc=dateobs) - tileid = get_next_tileid() fibermap, truth = get_targets(nspec, tileid=tileid) flux = truth['FLUX'] @@ -105,14 +110,16 @@ def new_exposure(nspec=5000, expid=None): #- Project flux to photons phot = thru.photons(wave[ii], flux[:,ii], units='1e-17 erg/s/cm2/A', - objtype=truth['OBJTYPE'], exptime=params['exptime']) + objtype=truth['OBJTYPE'], exptime=params['exptime'], + airmass=airmass) truth['PHOT_'+channel] = phot truth['WAVE_'+channel] = wave[ii] #- Project sky flux to photons - skyphot = thru.photons(wave[ii], skyflux[ii], units='1e-17 erg/s/cm2/A/arcsec2', - objtype='SKY', exptime=params['exptime']) + skyphot = thru.photons(wave[ii], skyflux[ii]*airmass, + units='1e-17 erg/s/cm2/A/arcsec2', + objtype='SKY', exptime=params['exptime'], airmass=airmass) #- 2D version ### truth['SKYPHOT_'+channel] = np.tile(skyphot, nspec).reshape((nspec, len(ii)))