Merge pull request #1 from desihub/sjb/randseed

merging sjb/randseed to master

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:

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,

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)))