###########################################################################
# Gemini NIFS data reduction script - PYRAF version
# Reduction for:  GENERAL BASELINE CALIBRATIONS
#
# DESCRIPTION:
#
# This script is a pyraf-compatible version of the original NIFS scripts
# written by Tracy Beck. To use this script, do the following:
#
# 1) In the call to gemlist below, update the raw_dir, rootname and file
# numbers corresponding to your calibration files.
#
# 2) In the directory where the files reside, run this script within pyraf
# using the following command:
#
# ---> pyexecute("NIFS_Basecalib.py")
#
# The script should then run through to the end without user interventio.
#
# Data used are available via the Gemini web pages:
# http://www.gemini.edu/sciops/instruments/nifs/data-format-and-reduction
#
# or via the Gemin Science Archive:
# http://www1.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/gsa/
#
# File rootname: N20100410 (but science and arc was taken on 20100401)
#   362-367     Lamps on Flats
#   368-372     Lamps off Flats
#   373-374     Darks for arc
#   375-376     Ronchi Flats
#       181     Arcs
#
# VERSION:
#  1.0: Adapted from original cl scripts by Tracy Beck. R.McDermid, 05Sep2012
#
# AUTHOR: R.McDermid (rmcdermid@gemini.edu)
#
###########################################################################

###########################################################################
# Current limitations:  The NIFS Baseline calibration reductions have     #
# not been well tested on data that was obtained with non-standard        #
# wavelength configurations.  (i.e., a different central wavelength       #
# setting than the default Z, J, H and K-band values of 1.05, 1.25, 1.65  #
# and 2.20 microns).                                                      #
###########################################################################

###########################################################################
# STEP 1:  Prepare IRAF  		                                  #
###########################################################################
# Import some useful python utilities
import sys
import getopt
import os
import time
import shutil
# Import the pyraf module and relevant packages
from pyraf import iraf
iraf.gemini()
iraf.nifs()
iraf.gnirs()
iraf.gemtools()
import pyfits

# Unlearn the used tasks
iraf.unlearn(iraf.gemini,iraf.gemtools,iraf.gnirs,iraf.nifs)

# Create a log file and back up the previous one if it already exists
log = 'Basecalib.log'
if os.path.exists(log):
    t = time.localtime()
    app = "_"+str(t[0])+str(t[1]).zfill(2)+str(t[2]).zfill(2)+'_'+ \
        str(t[3]).zfill(2)+':'+str(t[4]).zfill(2)+':'+str(t[5]).zfill(2)
    shutil.move(log,log+app)

iraf.set(stdimage='imt2048')

# Prepare the package for NIFS
iraf.nsheaders("nifs",logfile=log)

###########################################################################
# STEP 2:  Define Variables and Reduction Lists                           #
###########################################################################
# Set clobber to 'yes' for the script. This still does not make the gemini tasks
# overwrite files, so you will likely have to remove files if you re-run the script.
user_clobber=iraf.envget("clobber")
iraf.reset(clobber='yes')

#-------------------------------------------------------------------------
# BEGIN EDITING HERE
#
# CHANGE THESE LINES TO YOUR DIRECTORY AND FILE NAMES
raw_data = "/Data/NIFS/NIFS_Tutorial_Data/Cals/pytest/"
iraf.gemlist('N20100410S','362-367',Stdout='flatlist')
iraf.gemlist('N20100410S','368-372',Stdout='flatdarklist')
iraf.gemlist('N20100401S','181',Stdout='arclist')
iraf.gemlist('N20100410S','373-374',Stdout='arcdarklist')
iraf.gemlist('N20100410S','375-376',Stdout='ronchilist')
#
# END EDITING HERE
#--------------------------------------------------------------------------

# Set the file names for the mane calibration outputs. Just use the first name
# in the list of relevant files, which we get from the file lists just created on the disk
calflat=str(open("flatlist", "r").readlines()[0]).strip()
flatdark=str(open("flatdarklist", "r").readlines()[0]).strip()
arc=str(open("arclist", "r").readlines()[0]).strip()
arcdark=str(open("arcdarklist", "r").readlines()[0]).strip()
ronchiflat=str(open("ronchilist", "r").readlines()[0]).strip()

###########################################################################
# STEP 3:  Determine the shift to the MDF file                            #
###########################################################################

iraf.nfprepare(calflat,rawpath=raw_data,outpref="s", shiftx='INDEF',  
               shifty='INDEF',fl_vardq='no',fl_corr='no',fl_nonl='no', logfile=log)

###########################################################################
# STEP 4:  Make the Flat field and BPM                                    #
###########################################################################

iraf.nfprepare("@flatlist",rawpath=raw_data,shiftim="s"+calflat,
                     fl_vardq='yes',fl_int='yes',fl_corr='no',fl_nonl='no', logfile=log)

iraf.nfprepare("@flatdarklist",rawpath=raw_data,shiftim="s"+calflat,
                     fl_vardq='yes',fl_int='yes',fl_corr='no',fl_nonl='no', logfile=log)

iraf.gemcombine("n//@flatlist",output="gn"+calflat,fl_dqpr='yes',
                fl_vardq='yes',masktype="none",logfile=log)
iraf.gemcombine("n//@flatdarklist",output="gn"+flatdark,fl_dqpr='yes',
                fl_vardq='yes',masktype="none",logfile=log)

[iraf.nsreduce("gn"+calflat,fl_nscut='yes',fl_nsappw='yes',fl_vardq='yes',
                  fl_sky='no',fl_dark='no',fl_flat='no',logfile=log)
iraf.nsreduce("gn"+flatdark,fl_nscut='yes',fl_nsappw='yes',fl_vardq='yes',
                  fl_sky='no',fl_dark='no',fl_flat='no',logfile=log)


# creating flat image, final name = rnN....._sflat.fits
iraf.nsflat("rgn"+calflat,darks="rgn"+flatdark,flatfile="rgn"+calflat+"_sflat",
            darkfile="rgn"+flatdark+"_dark",fl_save_dark='yes',process="fit",
            thr_flo=0.15,thr_fup=1.55,fl_vardq='yes',logfile=log) 


#rectify the flat for slit function differences - make the final flat.
iraf.nsslitfunction("rgn"+calflat,"rgn"+calflat+"_flat",
                    flat="rgn"+calflat+"_sflat",dark="rgn"+flatdark+"_dark",combine="median",
                    order=3,fl_vary='no',logfile=log)


###########################################################################
# STEP 5:  Reduce the Arc and determine the wavelength solution           #
###########################################################################

iraf.nfprepare("@arclist", rawpath=raw_data, shiftimage="s"+calflat,
               bpm="rgn"+calflat+"_sflat_bpm.pl",fl_vardq='yes',fl_corr='no',fl_nonl='no', logfile=log)

iraf.nfprepare("@arcdarklist", rawpath=raw_data, shiftimage="s"+calflat,
               bpm="rgn"+calflat+"_sflat_bpm.pl",fl_vardq='yes',fl_corr='no',fl_nonl='no', logfile=log)

# Determine the number of input arcs and arc darks so that the
# routine runs automatically for single or multiple files.
nfiles = len(open("arclist").readlines())
if nfiles > 1: 
   iraf.gemcombine("n//@arclist",output="gn"+arc,
      fl_dqpr='yes',fl_vardq='yes',masktype="none",logfile=log)
else:
   iraf.copy("n"+arc+".fits","gn"+arc+".fits")

nfiles = len(open("arcdarklist").readlines())
if nfiles > 1:
   iraf.gemcombine("n//@arcdarklist",output="gn"+arcdark,
      fl_dqpr='yes',fl_vardq='yes',masktype="none",logfile=log)
else:
   iraf.copy("n"+arcdark+".fits","gn"+arcdark+".fits")

iraf.nsreduce("gn"+arc,outpr="r",darki="gn"+arcdark,flati="rgn"+calflat+"_flat",
   fl_vardq='no', fl_nscut='yes', fl_nsappw='yes', fl_sky='no', fl_dark='yes',fl_flat='yes', 
   logfile=log)

###########################################################################
#  DATA REDUCTION HINT -                                                  # 
# For the nswavelength call, the different wavelength settings            #
# use different vaues for some of the parameters. For optimal auto        #
# results, use:                                                           #
#                                                                         #
# K-band: thresho=50.0, cradius=8.0   -->  (gives rms of 0.1 to 0.3)      # 
# H-band: thresho=100.0, cradius=8.0  -->  (gives rms of 0.05 to 0.15)    #
# J-band: thresho=100.0               -->  (gives rms of 0.03 to 0.09)    # 
# Z-band: Currently not working very well for non-interactive mode        #
#                                                                         # 
# Note that better RMS fits can be obtained by running the wavelength     #
# calibration interactively and identifying all of the lines              #
# manually.  Tedious, but will give more accurate results than the        #
# automatic mode (i.e., fl_inter-).  Use fl_iner+ for manual mode.        #
#                                                                         # 
###########################################################################

# Determine the wavelength of the observation and set the arc coordinate
# file.  If the user wishes to change the coordinate file to a different
# one, they need only to change the "clist" variable to their line list
# in the coordli= parameter in the nswavelength call.
hdulist = pyfits.open("rgn"+arc+".fits")
band = hdulist[0].header['GRATING'][0:1]

if band == "Z":
    clist="nifs$data/ArXe_Z.dat"
    my_thresh=100.0
elif band == "K":
    clist="nifs$data/ArXe_K.dat"
    my_thresh=50.0
else:
    clist="gnirs$data/argon.dat"
    my_thresh=100.0

iraf.nswavelength("rgn"+arc, coordli=clist, nsum=10, thresho=my_thresh, 
   trace='yes',fwidth=2.0,match=-6,cradius=8.0,fl_inter='no',nfound=10,nlost=10,
   logfile=log)

###########################################################################
# STEP 6:                                                                 #
#  Trace the spatial curvature and spectral distortion in the Ronchi flat #
###########################################################################
iraf.nfprepare("@ronchilist",rawpath=raw_data, shiftimage="s"+calflat,
               bpm="rgn"+calflat+"_sflat_bpm.pl", fl_vardq='yes',fl_corr='no',fl_nonl='no',logfile=log)

# Determine the number of input Ronchi calibration mask files so that
# the routine runs automatically for single or multiple files.
nfiles = len(open("ronchilist").readlines())
if nfiles > 1:
    iraf.gemcombine("n//@ronchilist",output="gn"+ronchiflat,fl_dqpr='yes',
                    masktype="none",fl_vardq='yes',logfile=log)
else:
    iraf.copy("n"+ronchiflat+".fits","gn"+ronchiflat+".fits")

iraf.nsreduce("gn"+ronchiflat, outpref="r",  dark="rgn"+flatdark+"_dark",
              flatimage="rgn"+calflat+"_flat", fl_nscut='yes', fl_nsappw='yes', 
              fl_flat='yes', fl_sky='no', fl_dark='yes', fl_vardq='no',
              logfile=log)

iraf.nfsdist("rgn"+ronchiflat, 
             fwidth=6.0, cradius=8.0, glshift=2.8, 
             minsep=6.5, thresh=2000.0, nlost=3, 
             fl_inter='no',logfile=log)

###########################################################################
# Reset to user defaults                                                  #
###########################################################################
if user_clobber == "no":
    iraf.set(clobber='no')

###########################################################################
#		End of the Baseline Calibration reduction                 #
###########################################################################
#	                                                                  #
#  The final output files created from this script for later science      #
#  reduction have prefixes and file names of:                             #
#     1. Shift reference file:  "s"+calflat                               #
#     2. Flat field:  "rn"+calflat+"_flat"                                #
#     3. Flat BPM (for DQ plane generation):  "rn"+calflat+"_flat_bpm.pl" #
#     4. Wavelength referenced Arc:  "wrgn"+arc                            #
#     5. Spatially referenced Ronchi Flat:  "rgn"+ronchiflat               #
#     For this reduction,                                                 #
#        Shift ref. file =   sN20100410S0362.fits                         #
#        Flat field      =  rgnN20100410S0362_flat.fits                    #
#        Flat BPM        =  rgnN20100410S0362_sflat_bpm.pl                 #
#        Arc frame       =  wrgnN20100401S0181.fits                        #
#        Ronchi flat     =  rgnN20100410S0375.fits                         #
#	                                                                  #
#  NOTE:  Other important information for reducing the science data is    #
#    included in the "database" directory that is created and edited      #
#    within the above "nswavelength" and "nfsdist" IRAF calls. For a      #
#    proper science reduction to work (particularly the "nsfitcoords"     #
#    step), the science data must either be reduced in the same directory #
#    as the calibrations, or the whole "database" directory created by    #
#    this script must be copied into the working science reduction        #
#    directory.                                                           #
#                                                                         #
###########################################################################
