Coadd script

In [94]:

    

# imports
import glob
import pdb

import matplotlib.pyplot as plt

from scipy.interpolate import interp1d

from astropy.table import Table, vstack

from xastropy.xutils import xdebug as xdb


    

In [4]:

    

x1d_path = os.getenv('DROPBOX_DIR')+'/COS-LRG/tmp/'


    

In [6]:

    

x1d_files = glob.glob(x1d_path+'*x1d.fits')


    

In [9]:

    

tbl = Table.read(x1d_files[0])


    

    




WARNING: UnitsWarning: The unit 'angstrom' has been deprecated in the FITS standard. Suggested: 10**-1 nm. [astropy.units.format.utils]
WARNING: UnitsWarning: The unit 'erg' has been deprecated in the FITS standard. Suggested: cm2 g s-2. [astropy.units.format.utils]
WARNING: UnitsWarning: 'erg /s /cm**2 /angstrom' contains multiple slashes, which is discouraged by the FITS standard [astropy.units.format.generic]

In [15]:

    

tbl[0:1]


    

    
Out[15]:




<Table length=1>

SEGMENT
EXPTIME
NELEM
WAVELENGTH [16384]
FLUX [16384]
ERROR [16384]
GROSS [16384]
GCOUNTS [16384]
NET [16384]
BACKGROUND [16384]
DQ [16384]
DQ_WGT [16384]
s
Angstrom
erg / (Angstrom cm2 s)
erg / (Angstrom cm2 s)
ct / s
ct
ct / s
ct / s

str4
float64
int32
float64
float32
float32
float32
float32
float32
float32
int16
float32
FUVA
1128.192
16384
1185.08305117 .. 2496.14189625
0.0 .. 0.0
2.04881e-16 .. 6.06116e-17
0.0 .. 0.0
0.0 .. 0.0
0.0 .. 0.0
0.0 .. 0.0
128 .. 128
0.0 .. 0.0

In [20]:

    

# Load
sega_tbls = []
for x1d_file in x1d_files:
    tbl = Table.read(x1d_file)
    sega_tbls.append(tbl[0:1])


    

In [31]:

    

# Grab one wavelength array
wave = sega_tbls[0]['WAVELENGTH'][0,:].data
wave[0:5]


    

    
Out[31]:





array([ 1185.08305117,  1185.16248733,  1185.24192356,  1185.32135985,
        1185.40079622])

In [37]:

    

# Sum exposure time
total_time = np.zeros_like(wave)
for sega_tbl in sega_tbls:
    total_time += sega_tbl['DQ_WGT'][0,:]*sega_tbl['EXPTIME']
#xdb.xhist(total_time)


    

In [42]:

    

# Find DQmin for all exposures -- Why are we doing this step??
dqmin = np.ones_like(wave).astype(int) * 99999
for sega_tbl in sega_tbls:
    # Reset DQ
    dq = sega_tbl['DQ'][0,:].data
    reset_1024 = dq == 1024
    dq[reset_1024] = 2
    dqmin = np.minimum(dq, dqmin)


    

In [52]:

    

# Find DQ_WGT max for all exposures
DQWmax = np.zeros_like(wave)
for sega_tbl in sega_tbls:
    # Reset DQ
    dqw = sega_tbl['DQ_WGT'][0,:].data
    DQWmax = np.maximum(dqw, DQWmax)
#xdb.xhist(dqwmax)


    

In [63]:

    

# Generate calib values
total_counts = np.zeros_like(wave)
for sega_tbl in sega_tbls:
    # 
    total_counts += DQWmax * sega_tbl['GCOUNTS'][0,:]
xdb.xplot(wave, total_counts)


    

In [88]:

    

# Calibration
wave_calib, calib = [], []
for sega_tbl in sega_tbls:
    # 
    gddq = (sega_tbl['DQ'] > 0) & (sega_tbl['FLUX'] > 0)
    # Append
    wave_calib.append(sega_tbl['WAVELENGTH'][gddq].data.flatten())
    calib.append( (sega_tbl['NET'][gddq] / sega_tbl['FLUX'][gddq]).data)
#xdb.xhist(total_counts)


    

In [89]:

    

wave_calib = np.concatenate(wave_calib)
calib = np.concatenate(calib)
# sort
srt = np.argsort(wave_calib)
wave_calib = wave_calib[srt]
calib = calib[srt]
xdb.xplot(wave_calib, calib)


    

In [82]:

    

# Cut down
gdwv = wave_calib < 2100.  # Anything above that is junk


    

In [92]:

    

# Spline
f = interp1d(wave_calib[gdwv], calib[gdwv], bounds_error=False, fill_value=0.)  # cubic behaves badly


    

In [96]:

    

plt.clf()
ax = plt.gca()
ax.scatter(wave_calib[gdwv], calib[gdwv])
ax.plot(wave, f(wave))
plt.show()


    

In [ ]: