###########################################################################
# Gemini NIFS data reduction script
# Reduction for:  TELLURIC STANDARD 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_Telluric.py")
#
# Notes:
# - In this script, the telluric standard star was observed in a pattern of
# ABABABABA, where A = an exposure with the star in the NIFS field, and
# B = blank sky. If your telluric was observed with a different patter,
# you will have to change the telluriclist and skylist files.
#
# - A sky frame is constructed by median combining sky frames.  Editing
# the way the sky subtraction is done should be easy if telluric data
# were obtained by offsetting to the sky (In this case, just look at the
# way the NIFS Science reduction is constructed).
#
# - This reduction does the extraction of the telluric spectra
# in the "nfextract" step interactively, and requires user input at this
# stage. You will also need to have an image viewer running (e.g. ds9).
# The rest should run without intervention. You can also
# automate the 1D extraction if you provide the star centers to nfextract
# (see nfextract help for details).
#
# - 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: N20100401 (but we use daytime calibrations from 20100410)
# Files 138-146 : K-band telluric star observed in ABABABABA pattern
#
# VERSION:
#  1.0: Adapted from original cl scripts by Tracy Beck. R.McDermid, 05Sep2012
#
# AUTHOR: R.McDermid (rmcdermid@gemini.edu)
#
###########################################################################

###########################################################################
# 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 = 'Telluric.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                           #
#                                                                         #
#        CHANGE THESE LINES TO YOUR DIRECTORY AND FILE NAMES              #
###########################################################################
# 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 the "raw_data" and "cal_data" to correspond to the
# directory paths for your raw science frames and processed 
# calibrations,respectively.  If all of your data is in your 
# working directory, change both values to "" (null strings).
raw_data = "/Data/NIFS/NIFS_Tutorial_Data/HD141004/K/Telluric/pytest/"
cal_data   = "/Data/NIFS/NIFS_Tutorial_Data/Cals/pytest/"

# Provide here the processed calibration files, which were created by
# the NIFS_Basecalib reduction script:
calflat    = "N20100410S0362"
arc        = "N20100401S0181"
ronchiflat = "N20100410S0375"

# Here we generate the telluric and sky list files. The example observations
# had a pattern ABABABABA, where A=telluric, B=sky. Check the gemlist help for
# other ways to generate the lists. You can also look at the 'POFFSET' and 'QOFFSET'
# keywords to determine which frames were on source, since these are usually small (<2")
# for targets, and larger (>3") for sky frames.
iraf.gemlist('N20100401S','138-146x2',Stdout='telluriclist')
iraf.gemlist('N20100401S','139-145x2',Stdout='skylist')
#
# END EDITING HERE
#--------------------------------------------------------------------------

# Use the first telluric frame as the base name for the combined telluric spectrum
telluric=str(open("telluriclist", "r").readlines()[0]).strip()

###########################################################################
# STEP 3:  Get the Calibrations for the Reduction                         #
###########################################################################

# Copy required files and transformation database into the current working
# directory
iraf.copy(cal_data+"rgn"+ronchiflat+".fits",output="./")
iraf.copy(cal_data+"wrgn"+arc+".fits",output="./")
if not os.path.isdir("./database"):
    os.mkdir('./database/')
iraf.copy(cal_data+"database/*",output="./database/")

###########################################################################
# STEP 4:  Reduce the Telluric Standard                                   #
###########################################################################

# Prepare the data
iraf.nfprepare("@telluriclist",rawpath=raw_data,shiftim=cal_data+"s"+calflat,
  bpm=cal_data+"rgn"+calflat+"_sflat_bpm.pl",fl_vardq='yes',
  fl_int='yes',fl_corr='no',fl_nonl='no')

iraf.nfprepare("@skylist",rawpath=raw_data,shiftim=cal_data+"s"+calflat,
  bpm=cal_data+"rgn"+calflat+"_sflat_bpm.pl",fl_vardq='yes',
  fl_int='yes',fl_corr='no',fl_nonl='no')

# Make a median combined sky from the offset frames, which uses the telluric
# frame as the name, with _sky appended.
iraf.gemcombine("n//@skylist",output="n"+telluric+"_sky",
  fl_dqpr='yes',fl_vardq='yes',masktype="none",combine="median",logfile=log)

# Do the sky subtraction on all the individual frames. Read the list (then
# get rid of stupid '\n' character returns) first.
telluriclist=open("telluriclist", "r").readlines()
telluriclist=[word.strip() for word in telluriclist]
for image in telluriclist:
    iraf.gemarith ("n"+image, "-", "n"+telluric+"_sky", "sn"+image, fl_vardq="yes", 
                     logfile=log)

#reduce and flat field the data
iraf.nsreduce("sn@telluriclist",outpref="r", 
  flatim=cal_data+"rgn"+calflat+"_flat",
  fl_nscut='yes',fl_nsappw='no',fl_vardq='yes',fl_sky='no',fl_dark='no',fl_flat='yes',logfile=log)

#fix bad pixels from the DQ plane
iraf.nffixbad("rsn@telluriclist",outpref="b",logfile=log)

#derive the 2D to 3D spatial/spectral transformation
iraf.nsfitcoords("brsn@telluriclist",outpref="f", fl_int='no',
    lamptr="wrgn"+arc, sdisttr="rgn"+ronchiflat,logfile=log,lxorder=4,syorder=4)

#apply the transformation determined in the nffitcoords step
iraf.nstransform("fbrsn@telluriclist",outpref="t",logfile=log)

# Extract 1D spectra from the 2D data
# NOTE: In order to run nfextract interactively you need an image 
# viewer (like ds9) open. To run interactively, hit any to mark the aperture
# position and display the extracted spectra.  Hit "q" to continue
# to the next spectrum.
iraf.nfextract("tfbrsn@telluriclist",outpref="x",diameter=0.5, fl_int='yes',
   fl_zval='yes',z1=0,z2=10000.0, logfile=log)

iraf.onedspec()
iraf.mkspec("blackbody", "Blacbody", ncols=2048, nlines=1, func=3,
   start_wave=20000, end_wave=25000, temp=8000)
iraf.hedit('blackbody', 'CTYPE1', 'linear', add='yes', ver='no')
iraf.hedit('blackbody', 'CD1_1', 2.1340415478, add='yes', ver='no')
iraf.hedit("blackbody", "CRVAL1", 20005.32031, add='yes', ver='no')
iraf.hedit("blackbody", "WAT1_001", "wtyle=linear axtype=wave", add='yes', ver='no')

# Combine all the 1D spectra to one final output file
iraf.gemcombine("xtfbrsn//@telluriclist",output="gxtfbrsn"+telluric,
   statsec="[*]", combine="median",logfile=log,masktype="none",
   fl_vardq='yes')


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

###########################################################################
#          End of the Telluric Calibration Data Reduction                 #
#                                                                         #
#  The output of this reduction script is a 1-D spectrum used for         #
# telluric calibration of NIFS science data.  For this particular         #
# reduction the output file name is "gxtfbrsn"+telluric, or:              #
# gxtfbrsnN20100401S0138. The file prefixes are described below.          #
#                                                                         #
# g = gemcombined/gemarithed   n=nfprepared  s=skysubtracted              #
# r=nsreduced  b = bad pixel corrected  f= run through nffitcoords        # 
# t = nftransformed   x = extracted to a 1D spectrum                      #
#                                                                         #
# This script is meant to be a guideline as a method of a typical data    #
# reduction for NIFS frames.  Of course, NIFS PIs can add or skip steps   #
# in this reduction as they deem fit in order to reduce their particular  #
# datasets.                                                               #
#                                                                         #
###########################################################################