# GMOS Data Reduction Cookbook companion script to the chapter:
#   "Reduction of Images with IRAF" 
#
# IRAF CL command script to: 
# Process GMOS images for Messier 8, in program GS-2006B-Q-18.
#    2016-Mar-29  shaw@noao.edu
#
# The names for the relevant header keywords and their expected values are described 
# in the DRC chapter entitled "Supplementary Material"
#
# Perform the following starting in the parent work directory:
#   cd /path/to/work_directory
# 
# Fetch the Static BPM file from the tutorial, place it in your work directory, and uncompress: 
#   bpm_gmos-s_EEV_v1_2x2_img_MEF.fits
# 
# To run this script in an IRAF session, load the gemini, gemtools & gmos packages, 
# then enter the following at the IRAF prompt: 
#   cl < gmos_img_proc.cl

print ("### Begin Processing GMOS/MOS Images ###")
print ("###")

cd ./raw

print ("===Exosure Selection===")
print (" --Selecting Bias Exposures--")
#
## Select bias exposures. 
# Bias exposures with a common observation class, RoI, and CCD binning:
s1 = "obstype=='BIAS' && obsclass=='dayCal' && detrO1ys>1024 && ccdsum=='2 2'"

# Select bias exposures within a few weeks of the target observations:
s2 = "&& @'date-obs' > '2006-09-01' && @'date-obs' < '2006-10-30'"
string biasSelect
biasSelect = s1 // s2

# Select exposures using information from both PHDU and HDU-1
hselect ("*.fits[1,inherit=yes]", "$I", biasSelect, > "bias_tmp.txt")

# Now remove the IRAF index and kernel parameters from the file names
string omit = "exten,index,kernel"
gemextn ("@bias_tmp.txt", omit=omit, outfile="biasFiles.txt")

### Flat-fields
#
print (" --Selecting Flat-field Expsoures--")

# Flat-fields must match the Object and observation class.
# Note: IRAF refers to the OBJECT keyword with the special name: "i_title" 
s1 = "i_title?='Twilight' && obsclass?='dayCal' "

# Flat-fields must also match the RoI, CCD binning, grating and aperture:
s2 = "&& detro1ys>1024 && ccdsum?='2 2' && grating?='MIRROR' && maskname?='None' "

# Select flats obtained contemporaneously with the observations
string date
date = "&& @'date-obs' > '2006-09-10' && @'date-obs' < '2006-10-10'"
string flt
flt = s1 // s2 // date

# Generate lists of calibration filenames using above criteria, per filter.
hselect ("S*.fits[1,inherit=yes]", "$I", (flt // "&& filter2?='r_G0326'"), > "flt_tmp.txt")
gemextn ("@flt_tmp.txt", omit=omit, outfile="flat_r.txt")

hselect ("S*.fits[1,inherit=yes]", "$I", (flt // "&& filter2?='i_G0327'"), > "flt_tmp.txt")
gemextn ("@flt_tmp.txt", omit=omit, outfile="flat_i.txt")

hselect ("S*.fits[1,inherit=yes]", "$I", (flt // "&& filter2?='Ha_G0336'"), > "flt_tmp.txt")
gemextn ("@flt_tmp.txt", omit=omit, outfile="flat_ha.txt")

hselect ("S*.fits[1,inherit=yes]", "$I", (flt // "&& filter2?='HaC_G0337'"), > "flt_tmp.txt")
gemextn ("@flt_tmp.txt", omit=omit, outfile="flat_haC.txt")

hselect ("S*.fits[1,inherit=yes]", "$I", (flt // "&& filter2?='SII_G0335'"), > "flt_tmp.txt")
gemextn ("@flt_tmp.txt", omit=omit, outfile="flat_s2.txt")

### Process science exposures in the parent directory
print (" --Selecting Science Exposures--")

## Select science exposures via obs. type & class, and...
string sci
s1 = "obstype?='OBJECT' && obsclass?='science' "

# ...also by CCD RoI and binning, and by grating & (no)mask:
s2 = "&& detro1ys>1024 && ccdsum?='2 2' && grating?='MIRROR' && maskname?='None'"
sci = s1 // s2

# Build list of science exposures 
hselect ("S*.fits[1,inherit=yes]", "$I", (sci // "&& filter2?='r_G0326'"), > "sci_tmp.txt")
gemextn ("@sci_tmp.txt", omit=omit, outfile="sci_r.txt")

hselect ("S*.fits[1,inherit=yes]", "$I", (sci // "&& filter2?='i_G0327'"), > "sci_tmp.txt")
gemextn ("@sci_tmp.txt", omit=omit, outfile="sci_i.txt")

hselect ("S*.fits[1,inherit=yes]", "$I", (sci // "&& filter2?='Ha_G0336'"), > "sci_tmp.txt")
gemextn ("@sci_tmp.txt", omit=omit, outfile="sci_ha.txt")

hselect ("S*.fits[1,inherit=yes]", "$I", (sci // "&& filter2?='HaC_G0337'"), > "sci_tmp.txt")
gemextn ("@sci_tmp.txt", omit=omit, outfile="sci_haC.txt")

hselect ("S*.fits[1,inherit=yes]", "$I", (sci // "&& filter2?='SII_G0335'"), > "sci_tmp.txt")
gemextn ("@sci_tmp.txt", omit=omit, outfile="sci_s2.txt")

# Clean up.
dele ("*_tmp.txt")
rename ("biasFiles.txt,flat*.txt,sci*.txt", "../")

print ("===Creating MasterCals===")
#
## Create Bias Residual, including VAR and DQ arrays.
# Use primarily the default task parameters.
print (" --Creating Bias Residual MasterCals--")
cd ..

# Use primarily the default task parameters.
unlearn gbias
unlearn gemextn    # Disarm a bug in gbias
gbias.log = "gbiasLog.txt"
gbias.rawpath="./raw"
gbias.verbose=no

# Create the VAR and DQ arrays as well.
gbias ("@biasFiles.txt", "MCbias.fits", fl_vardq+)

## Process the raw flat-field exposures
print (" --Creating Flat-field MasterCals--")

unlearn giflat
giflat.logfile="giflatLog.txt"
# Caution: need to add a trailing slash to giflat.rawpath
giflat.rawpath="./raw/"
# The input BPM appears to make no difference. 
#giflat.bpm="bpm_gmos-s_EEV_v1_2x2_img_MEF.fits"
giflat.verbose=no

# Create the VAR and DQ arrays as well.
print ("  - H-alpha flat -")
giflat ("@flat_ha.txt",  "MCflat_ha.fits",  bias="MCbias.fits", fl_vardq+)
print ("  - H-alpha continuum flat -")
giflat ("@flat_haC.txt", "MCflat_haC.fits", bias="MCbias.fits", fl_vardq+)
print ("  - [S_II] flat -")
giflat ("@flat_s2.txt",  "MCflat_s2.fits",  bias="MCbias.fits", fl_vardq+)
print ("  - r flat -")
giflat ("@flat_r.txt",   "MCflat_r.fits",   bias="MCbias.fits", fl_vardq+)
print ("  - i flat -")
giflat ("@flat_i.txt",   "MCflat_i.fits",   bias="MCbias.fits", fl_vardq+)

# After reviewing the processing and output, delete intermediate files.
imdele ("rgS*.fits,gS*.fits")

### Perform basic processing
print ("===Processing Science Files===")
print (" --Performing Basic Processing--")

# Use primarily the default task parameters.
unlearn gireduce
gireduce.logfile="gireduceLog.txt"
gireduce.rawpath="./raw"
gireduce.bpm="bpm_gmos-s_EEV_v1_2x2_img_MEF.fits"
gireduce.verbose=no

# Perform basic reductions on exposures, 1 filter at a time.
print (" --- H-alpha ---")
gireduce ("@sci_ha.txt", bias="MCbias", fl_vardq+, flat1="MCflat_ha")
print (" --- H-alpha continuum---")
gireduce ("@sci_haC.txt", bias="MCbias", fl_vardq+, flat1="MCflat_haC")
print (" --- [S_II] ---")
gireduce ("@sci_s2.txt", bias="MCbias", fl_vardq+, flat1="MCflat_s2")
print (" --- r ---")
gireduce ("@sci_r.txt", bias="MCbias", fl_vardq+, flat1="MCflat_r")
print (" --- i ---")
gireduce ("@sci_i.txt", bias="MCbias", fl_vardq+, flat1="MCflat_i")
imdele ("gS*.fits")

print (" --Mosaicing extensions--")
# Now mosaic the extensions within each file
# Use primarily the default task parameters.
unlearn gemextn   # disarms a bug in gmosaic
unlearn gmosaic
gmosaic.logfile="gmosaicLog.txt"
gmosaic.verb=no

# gmosaic does not appear to work reliably with @filelist input. 
#
# The following approaches have failed (intermitantly):
# Create the list of images to process, and mosaic the extensions per file
#concatenate ("sci_*.txt", "sci_files.txt")
#sections rg//@sci_files.txt >rgsci_files.txt
#gmosaic ("@rgsci_files.txt", geointer="nearest", logfile="gmosaicLog.txt")
# Nor does this:
#gmosaic ("rg@sci_ha.txt", geointer="nearest", fl_vardq+, fl_fulldq+)

# So, we use explicit loops: 
sections ("rg//@sci_ha.txt", >"rgsci_ha.txt")
list = "rgsci_ha.txt"
while (fscan (list, s1) != EOF) {
    gmosaic (s1, geointer="nearest", fl_vardq+, fl_fulldq+)
}
sections ("rg//@sci_haC.txt", >"rgsci_haC.txt")
list = "rgsci_haC.txt"
while (fscan (list, s1) != EOF) {
    gmosaic (s1, geointer="nearest", fl_vardq+, fl_fulldq+)
}
sections ("rg//@sci_s2.txt", >"rgsci_s2.txt")
list = "rgsci_s2.txt"
while (fscan (list, s1) != EOF) {
    gmosaic (s1, geointer="nearest", fl_vardq+, fl_fulldq+)
}
sections ("rg//@sci_r.txt", >"rgsci_r.txt")
list = "rgsci_r.txt"
while (fscan (list, s1) != EOF) {
    gmosaic (s1, geointer="nearest", fl_vardq+, fl_fulldq+)
}
sections ("rg//@sci_i.txt", >"rgsci_i.txt")
list = "rgsci_i.txt"
while (fscan (list, s1) != EOF) {
    gmosaic (s1, geointer="nearest", fl_vardq+, fl_fulldq+)
}
#imdele ("rgS*.fits,gS*.fits")

### Advanced processing
print (" --Performing Advanced Processing--")

## Image Stacking
# Select the files per object region, and per filter. 
# WARNING: create all of these lists first before processing to avoid name conflicts.
s1 = "i_title?='M8-1' && "
hselect ("mrg*.fits[0]", "$I", (s1 // "filter2?='r_G0326'"), > "mrg_tmp.txt")
gemextn ("@mrg_tmp.txt", omit=omit, outfile="m8-1_r.txt")
hselect ("mrg*.fits[0]", "$I", (s1 // "filter2?='i_G0327'"), > "mrg_tmp.txt")
gemextn ("@mrg_tmp.txt", omit=omit, outfile="m8-1_i.txt")
hselect ("mrg*.fits[0]", "$I", (s1 // "filter2?='Ha_G0336'"), > "mrg_tmp.txt")
gemextn ("@mrg_tmp.txt", omit=omit, outfile="m8-1_ha.txt")
hselect ("mrg*.fits[0]", "$I", (s1 // "filter2?='HaC_G0337'"), > "mrg_tmp.txt")
gemextn ("@mrg_tmp.txt", omit=omit, outfile="m8-1_haC.txt")
hselect ("mrg*.fits[0]", "$I", (s1 // "filter2?='SII_G0335'"), > "mrg_tmp.txt")
gemextn ("@mrg_tmp.txt", omit=omit, outfile="m8-1_s2.txt")

s1 = "i_title?='M8-2' && "
hselect ("mrg*.fits[0]", "$I", (s1 // "filter2?='r_G0326'"), > "mrg_tmp.txt")
gemextn ("@mrg_tmp.txt", omit=omit, outfile="m8-2_r.txt")
hselect ("mrg*.fits[0]", "$I", (s1 // "filter2?='i_G0327'"), > "mrg_tmp.txt")
gemextn ("@mrg_tmp.txt", omit=omit, outfile="m8-2_i.txt")
hselect ("mrg*.fits[0]", "$I", (s1 // "filter2?='Ha_G0336'"), > "mrg_tmp.txt")
gemextn ("@mrg_tmp.txt", omit=omit, outfile="m8-2_ha.txt")
hselect ("mrg*.fits[0]", "$I", (s1 // "filter2?='HaC_G0337'"), > "mrg_tmp.txt")
gemextn ("@mrg_tmp.txt", omit=omit, outfile="m8-2_haC.txt")
hselect ("mrg*.fits[0]", "$I", (s1 // "filter2?='SII_G0335'"), > "mrg_tmp.txt")
gemextn ("@mrg_tmp.txt", omit=omit, outfile="m8-2_s2.txt")

s1 = "i_title?='M8-3' && "
hselect ("mrg*.fits[0]", "$I", (s1 // "filter2?='r_G0326'"), > "mrg_tmp.txt")
gemextn ("@mrg_tmp.txt", omit=omit, outfile="m8-3_r.txt")
hselect ("mrg*.fits[0]", "$I", (s1 // "filter2?='i_G0327'"), > "mrg_tmp.txt")
gemextn ("@mrg_tmp.txt", omit=omit, outfile="m8-3_i.txt")
hselect ("mrg*.fits[0]", "$I", (s1 // "filter2?='Ha_G0336'"), > "mrg_tmp.txt")
gemextn ("@mrg_tmp.txt", omit=omit, outfile="m8-3_ha.txt")
hselect ("mrg*.fits[0]", "$I", (s1 // "filter2?='HaC_G0337'"), > "mrg_tmp.txt")
gemextn ("@mrg_tmp.txt", omit=omit, outfile="m8-3_haC.txt")
hselect ("mrg*.fits[0]", "$I", (s1 // "filter2?='SII_G0335'"), > "mrg_tmp.txt")
gemextn ("@mrg_tmp.txt", omit=omit, outfile="m8-3_s2.txt")

dele ("mrg_tmp.txt")

## Co-add the images:
print ("  -Begin image co-addition-")

# Use primarily the default task parameters.
unlearn imcoadd
imcoadd.logfile="imcoaddLog.txt"
imcoadd.fwhm=3
imcoadd.datamax=6.e4
imcoadd.geointer="nearest"

imcoadd ("@m8-1_r.txt",   outimage="m8-1_r.fits")
imcoadd ("@m8-1_i.txt",   outimage="m8-1_i.fits")
imcoadd ("@m8-1_ha.txt",  outimage="m8-1_ha.fits")
imcoadd ("@m8-1_haC.txt", outimage="m8-1_haC.fits")
imcoadd ("@m8-1_s2.txt",  outimage="m8-1_s2.fits")

imcoadd ("@m8-2_r.txt",   outimage="m8-2_r.fits")
imcoadd ("@m8-2_i.txt",   outimage="m8-2_i.fits")
imcoadd ("@m8-2_ha.txt",  outimage="m8-2_ha.fits")
imcoadd ("@m8-2_haC.txt", outimage="m8-2_haC.fits")
imcoadd ("@m8-2_s2.txt",  outimage="m8-2_s2.fits")

imcoadd ("@m8-3_r.txt",   outimage="m8-3_r.fits")
imcoadd ("@m8-3_i.txt",   outimage="m8-3_i.sfits")
imcoadd ("@m8-3_ha.txt",  outimage="m8-3_ha.fits")
imcoadd ("@m8-3_haC.txt", outimage="m8-3_haC.fits")
imcoadd ("@m8-3_s2.txt",  outimage="m8-3_s2.fits")

# Inspect the co-added images for quality. Tweak task parameters as necessary. 
# Delete the intermediate files, unless they are needed to reprocess
dele ("*_trn*,*_pos,*_cen")
#imdele ("*badpix.pl,*_med.fits,*_mag.fits")
#imdele ("mrgS*.fits")

# r-band backgrounds need to be adjusted. About: -140, 0, -50 for pos. 1,2,3