##^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009
#! /usr/bin/env python
#
#This script will run DAOFind and DAOPhot on the images in FitsFileList.
#It will generate the .coo and .mag files and put them in the same directory
#that the file list is placed in.   
#
#Some of the keywords in the fits header are assumed to be set.  They are
#the following: 
# 
#AIRMASS- the average airmass according to Stetson's formula using Simpson's rule
#ITIME  - the integration time 
#UTMID  - the universal time at the middle of the dither sequence 
#FILTER - the filter used 
#
#KEYWORDS:
#  TvFlag      = 0 - Don't plot or display
#	         1 - Display and offer user imexamine
#  OverRide    = 0 - Don't overwrite the existing .coo or .mag files
#	         1 - Overwrite the files
#  OverWrite   = 0 - Read the paramters from the existing .par files
#		 1 - Rewrite the parameter files with the values in this file
#  SkipPSF     = 1 - Skip the PSF fitting
#  SkipNStar   = 1 - Skip the NStar
#  SkipDAOFind = 1 - Skip the DAOFind routine
#  SkipAllStar = 1 - Skip the AllStar routine
#  SkipPSF     = 1 - Skip the first psf fitting
#  SkipPSF2    = 1 - Skip the fitting of the psf after subtraction of neighbors
#  Comment     = ''- Comment to describe comditions in psf 
#Example Calling Sequence:
#run ~/python/H4RGAnly/BatchDAOFind.py '/nfs/slac/g/ki/ki09/lances/Reduced/07Apr27/Landolt110956/Astrometry/Landolt110956ISlopes.lst'

import sys,os,getopt
from pyraf import iraf 
import pyfits
from pylab import *
import numpy
import numdisplay

#The filename will be specified as the first argument. These files have
#the form [Object Name]_[Filter Used]_SlopeFlatFielded...fits
FitsFileList = sys.argv[1]

#Now parse the command line keyword arguments 
keywords = ['TvFlag=', 'SkipDAOFind=', 'SkipDAOPhot=', 'SkipPSF=', \
	    'SkipNStar=', 'SkipSubStar=','SkipAllStar=','OverRide=', \
            'SkipPSF2=', 'UsePSFList=','OverWrite=','Comment=' ]
OverWrite   = 0
OverRide    = 0
TvFlag      = 0
SkipDAOFind = 0
SkipDAOPhot = 0
Comment     = ''

opts, extraparams = getopt.getopt(sys.argv[2:],'',keywords)
for o,p in opts:
  if o in ['--TvFlag']:
    TvFlag=int(p)
  elif o in ['--OverRide']:
    OverRide=int(p)
  elif o in ['--SkipDAOFind']:
    SkipDAOFind = int(p)
  elif o in ['--SkipDAOPhot']:
    SkipDAOPhot = int(p)
  elif o in ['--OverWrite']:
    OverWrite = int(p)
  elif o in ['--UsePSFList']:
    UsePSFList = int(p)
  elif o in ['--Comment']:
    Comment = p

#Import the good packages from IRAF
iraf.digiphot(_doprint=0)
iraf.daophot(_doprint=0)
iraf.apphot(_doprint=0)

#Extract various keys from the full path
BaseName     = os.path.basename(FitsFileList)
DirName      = os.path.dirname(FitsFileList)
OutDir       = DirName+'/'
os.chdir(DirName)

#Get the filenames as a list of strings 
iraf.tproject.setParam('intable',FitsFileList)
iraf.tproject.setParam('outtable','STDOUT')
iraf.tproject.setParam('columns','c1')
FitsFileNames = iraf.tproject(Stdout=1, mode ='h')

#The names have a [0] appended to them
FitsFileName  = (FitsFileNames[1].split('['))[0]
FitsHDU    = pyfits.open(FitsFileName)
FitsHeader = FitsHDU[0].header
ThisFilter = FitsHeader['FILTER']

BaseName      = os.path.basename(FitsFileNames[1])
BaseParts     = BaseName.split('_')
ObjectName    = BaseParts[0]
BaseKey       = BaseParts[0]

#Get some of the keywords from the keywords file 
if os.path.isfile('/afs/slac/u/ki/lances/python/H4RGAnly/keywords.py'):
  execfile('/afs/slac/u/ki/lances/python/H4RGAnly/keywords.py');

#GET THE DATE
Fullpath=os.getcwd()

#PARAMETERS in the various files
datapars   = iraf.datapars.getParList()
daopars    = iraf.daopars.getParList()
findpars   = iraf.findpars.getParList()
centerpars = iraf.centerpars.getParList()
skypars    = iraf.fitskypars.getParList()
photpars   = iraf.photpars.getParList()

#KEY PARAMETERS Paramters to use
FWHM      = 8.75		  #Full Width at Half Max 
GAIN      = 5.0 		  #e-/ADU
READNOISE = 0.45		  #Readnoise for a single
NUMREADS  = 10			  #Number of reads in the ramp
NUMFILES  = 1 			  #Number of files used to make dither 
GAIN      = NUMFILES*GAIN         #If N images are averaged, GAIN=N*GAIN
print 'Using FWHM of ' + str(FWHM) 

#PHOTPARS #############################################################
PhotParsFileName = OutDir+ObjectName+'_'+ThisFilter+Comment+'PhotPar.par'
if OverWrite == 1 or not os.path.exists(PhotParsFileName) :
  iraf.photpars.setParam('apertures',FWHM)
  iraf.photpars.saveParList(filename=PhotParsFileName)
else :
  iraf.photpars.setParList(ParList = PhotParsFileName)

#CENTERPARS #############################################################
CenterParsFileName = OutDir+ObjectName+'_'+ThisFilter+Comment+\
	'CenterPar.par'
if OverWrite == 1 or not os.path.exists(CenterParsFileName) :
  iraf.centerpars.setParam('calgorithm','none') #Recommended for first pass
  iraf.centerpars.setParam('cbox',4.5*FWHM)     #Centroiding box
  iraf.centerpars.setParam('cthreshold','0')
  iraf.centerpars.saveParList(filename=CenterParsFileName)
else :
  iraf.centerpars.setParList(ParList=CenterParsFileName)

#SKYPARS ################################################################
FitSkyParsFileName = OutDir+ObjectName+'_'+ThisFilter+Comment+\
	'FitSkyPar.par'
if OverWrite == 1 or not os.path.exists(FitSkyParsFileName) :
  iraf.fitskypars.setParam('salgorithm','mode')
  iraf.fitskypars.setParam('skyvalue',2.1)
  iraf.fitskypars.setParam('annulus',4.*FWHM)
  iraf.fitskypars.setParam('dannulus',2*FWHM)
  iraf.fitskypars.saveParList(filename=FitSkyParsFileName)
else:
  iraf.fitskypars.setParList(ParList=FitSkyParsFileName)

#FINDPARS ##############################################################
FindParsFileName =  OutDir+ObjectName+'_'+ThisFilter+Comment+\
	'FindPar.par'
if OverWrite == 1 or not os.path.exists(FindParsFileName) :
  iraf.findpars.setParam('threshold',10.0)
  iraf.findpars.setParam('sharphi',  3)
  iraf.findpars.setParam('roundhi', 1.3)
  iraf.findpars.setParam('roundlo', -1.3)
  iraf.findpars.saveParList(filename=FindParsFileName)
else :
  iraf.findpars.setParList(ParList = FindParsFileName)

#DATAPARS - http://stsdas.stsci.edu/cgi-bin/gethelp.cgi?datapars.hlp
DataParsFileName =  OutDir+ObjectName+'_'+ThisFilter+Comment+\
	'DataPar.par'
if OverWrite == 1 or not os.path.exists(DataParsFileName) :
  iraf.datapars.setParam('fwhmpsf',FWHM)    #FWHM of biggest star
  iraf.datapars.setParam('sigma','0.34')    #Std. Dev of sky pixels
  iraf.datapars.setParam('datamin','8')   #Min good data value ~Sky-3Sigma
  iraf.datapars.setParam('datamax','4500')  #Max good data value
  iraf.datapars.setParam('epadu'  , GAIN)   #e- per ADU; calculated above
  iraf.datapars.setParam('readnoise','.43') #readnoise
  #KEYWORDS IN FITSHEADER
  iraf.datapars.setParam('exposure','ITIME')	#Integration time
  iraf.datapars.setParam('airmass', 'AIRMASS')	#Airmass at middle
  iraf.datapars.setParam('filter', 'FILTER')	#Filter used
  iraf.datapars.saveParList(filename=DataParsFileName)
else:
  iraf.datapars.setParList(ParList=DataParsFileName)

#DAOPARS ###############################################################
DAOParsFileName =  OutDir+ObjectName+'_'+ThisFilter+Comment+\
	'DAOPar.par'
if OverWrite == 1 or not os.path.exists(DAOParsFileName) :
  iraf.daopars.setParam('psfrad',4.*FWHM)
  iraf.daopars.setParam('fitrad',1.*FWHM)
  iraf.daopars.setParam('function','auto')
  iraf.daopars.setParam('nclean',2)
  iraf.daopars.setParam('sannulus',1.5*FWHM)
  iraf.daopars.setParam('wsannulus',2*FWHM)
  iraf.daopars.setParam('fitsky','no')    #Refit the sky
  iraf.daopars.setParam('varorder',0)      #Use an analytic model
  iraf.daopars.setParam('recenter','yes')
  iraf.daopars.saveParList(filename=DAOParsFileName)
else: 
  iraf.daopars.setParList(ParList = DAOParsFileName)

########################################################################
########################################################DAOFIND
#Set up daofind to go without prompting for input
iraf.daofind.setParam('image','@'+FitsFileList)		#Set ImageName
iraf.daofind.setParam('verify','no')			#Don't verify
iraf.daofind.setParam('interactive','no')		#Interactive
iraf.daofind.saveParList(filename=OutDir+ObjectName+\
        '_'+ThisFilter+'daofind.par')
iraf.daofind.setParam('output','default')               #Output File

#Run DAOFIND once unless OverRide = 1
iraf.daofind.setParam('graphics', 'stdgraph')
iraf.daofind.setParam('display', 'display(image=\'$1\',frame=$2)')
if OverRide == 1 and SkipDAOFind != 1: 
   print 'Running DAOFIND'
   iraf.daofind(mode='h')
else :
   print 'Skipping DAOFind'

#Mark the centroids on the image if TvFlag = 1
if TvFlag >= 1:
  MaxDisp = Im.mean()+3*Im.std()
  iraf.set(stdimage='imt4096')
  iraf.display(FitsFileNames[1],1,z1=1,z2=MaxDisp,zrange='no')
  if os.path.exists(DaoFindFile):
     iraf.tvmark(1, DaoFindFile)

#########################################################DAOPHOT
#Run DAOPHOT once unless OverRide = 1
iraf.phot.setParam('image',FitsFileList)
iraf.phot.setParam('coords','default')
iraf.phot.setParam('verify','no')
iraf.phot.setParam('output','default')
iraf.phot.setParam('interactive','no')
iraf.phot.saveParList(filename=OutDir+ObjectName+\
        '_'+ThisFilter+'daophot.par')

if OverRide ==1 and SkipDAOPhot !=1 : 
   print 'Running DAOPhot'
   photstring = iraf.phot(mode='h',Stdout=1)
   DaoPhotTxtFile=OutDir+BaseKey+Comment+'.mag.txt'
else :
   print 'Skipping DAOPhot'