<jupyter_start><jupyter_text>Coadd script<jupyter_code># 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<jupyter_output><empty_output><jupyter_text>Testing<jupyter_code>x1d_path = os.getenv('DROPBOX_DIR')+'/COS-LRG/tmp/'
x1d_files = glob.glob(x1d_path+'*x1d.fits')<jupyter_output><empty_output><jupyter_text>Check an x1d file<jupyter_code>tbl = Table.read(x1d_files[0])
tbl[0:1]<jupyter_output><empty_output><jupyter_text>Segment A<jupyter_code># Load
sega_tbls = []
for x1d_file in x1d_files:
    tbl = Table.read(x1d_file)
    sega_tbls.append(tbl[0:1])
# Grab one wavelength array
wave = sega_tbls[0]['WAVELENGTH'][0,:].data
wave[0:5]
# 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)
# 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)
# 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)
# 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)<jupyter_output><empty_output><jupyter_text>Calibration<jupyter_code># 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)
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)
# Cut down
gdwv = wave_calib < 2100.  # Anything above that is junk
# Spline
f = interp1d(wave_calib[gdwv], calib[gdwv], bounds_error=False, fill_value=0.)  # cubic behaves badly
plt.clf()
ax = plt.gca()
ax.scatter(wave_calib[gdwv], calib[gdwv])
ax.plot(wave, f(wave))
plt.show()<jupyter_output><empty_output>