##^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009
#DAOFIND_HxRG
#
#PURPOSE:
#  To examine a slopefitted file from one of the HxRG detectors.  The 
#  script will generate a text file with the path 
#
#    HxRG.DataDir+HxRG.DetStr+HxRG.ElecSTr+HxRG.Date+HxRG.Object+'/Photometry'
#
#  The text file contains the positions of the stars found and can be used as 
#  input to DAOPhot_HxRG.py or RadProf_HxRG.py.
#
#  Also, if TvFlag==1, interactive mode is entered and you can use the 
#  examine functions from IRAF to looks at the FWHM of the stars.  If 
#  you do this, make sure ds9 is already open.
#
#INPUTS:
#  FitsFileName: string
#     The full path to the fits file that will be examined.
#  
#KEYWORDS:
#  TvFlag: int
#          0 - Don't plot or display
#          1 - Display and offer user imexamine
#  OverRide: int    
#	   0 - Don't overwrite the existing .coo or .mag files
#          1 - Overwrite the files
#  OverWrite: int   
#	   0 - Read the paramters from the existing .par files
#          1 - Rewrite the parameter files with the values in this file
#
#NOTE:
#  It seems that giving a small FWHM is preferable to a large one in order
#  to identify a wide range of stars.
#
#Example Calling Sequence:
#run DAOFind_HxRG.py '/nfs/slac/g/ki/ki03/lances/H2RG-32-147/ASIC/Reduced/07Dec18/NGC956/SaturationRamp_I_H2RG_SIPIN_5_Reads_Dec18_2007_20_22_37_I_Slope.fits' --TvFlag=0
#run DAOFind_HxRG.py '/nfs/slac/g/ki/ki03/lances/H2RG-32-147/ASIC/Reduced/07Dec18/NGC956/SaturationRamp_I_H2RG_SIPIN_10_Reads_Dec18_2007_20_44_42_I_RPS_Slope.fits' --OverRide=1 --OverWrite=1

##SETUP:
execfile('/afs/slac/u/ki/lances/python/python_HxRG/HxRG_Setup.py')
from HxRG_Class import *
import pdb

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

#Form the class that will contain the necessary info
HxRG = HxRG_C(FitsFileName=FitsFileName)

#Now parse the command line keyword arguments 
keywords = ['TvFlag=', 'OverRide=', 'OverWrite=' ]
OverWrite   = 0
OverRide    = 0
TvFlag      = 0

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 ['--OverWrite']:
    OverWrite = int(p)


#Creat the Photometry Directory
if not os.path.exists(HxRG.PhotDir): os.mkdir(HxRG.PhotDir)

#Read in the array
HxRG.Get_Raw_Header(FitsFileName)
Im = pyfits.getdata(FitsFileName)

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

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

#KEY PARAMETERS Paramters to use
FWHM      = 4.5                   #Full Width at Half Max 
SKY       = 0.8		  #Average value of the sky
DATAMIN   = 0.1                    #Minimum value for acceptable data
DATAMAX   = 15000		  #Maximum value for acceptable data
GAIN      = 5.0                  #e-/ADU
READNOISE = 0.15                     #Readnoise for a single
BACKSIG   = 3
NUMREADS  = 15                    #Number of reads in the ramp
NUMFILES  = 9.0                   #Number of files used to make dither 
GAIN      = NUMFILES*GAIN         #If N images are averaged, GAIN=N*GAIN
print 'Using FWHM of ' + str(FWHM)

#CENTERPARS #############################################################
CenterParsFileName = HxRG.PhotDir+HxRG.ObjectName+'_'+HxRG.Filter+\
  '_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 = HxRG.PhotDir+HxRG.ObjectName+'_'+HxRG.Filter+\
  '_FitSkyPar.par'
if OverWrite == 1 or not os.path.exists(FitSkyParsFileName) :
  iraf.fitskypars.setParam('salgorithm','mode')
  iraf.fitskypars.setParam('skyvalue', SKY)
  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 = HxRG.PhotDir+HxRG.ObjectName+'_'+HxRG.Filter+\
   '_FindPar.par'
if OverWrite == 1 or not os.path.exists(FindParsFileName) :
  iraf.findpars.setParam('threshold',  1.0)
  iraf.findpars.setParam('sharphi',    0.5)	#0.5-1.0 likely cosmic ray, bad
  iraf.findpars.setParam('roundhi',    0.8)
  iraf.findpars.setParam('roundlo',   -0.8)
  iraf.findpars.saveParList(filename=FindParsFileName)
else :
  iraf.findpars.setParList(ParList = FindParsFileName)

#DATAPARS -#############################################################
DataParsFileName = HxRG.PhotDir+HxRG.ObjectName+'_'+HxRG.Filter+\
   '_DataPar.par'
if OverWrite == 1 or not os.path.exists(DataParsFileName) :
  iraf.datapars.setParam('fwhmpsf', FWHM)       #FWHM of biggest star
  iraf.datapars.setParam('sigma',   BACKSIG)    #Std. Dev of sky pixels
  iraf.datapars.setParam('datamin', DATAMIN)    #Min good data value ~Sky-3Sigma
  iraf.datapars.setParam('datamax', DATAMAX)  #Max good data value
  iraf.datapars.setParam('epadu',   GAIN)     #e- per ADU; calculated above
  iraf.datapars.setParam('readnoise', READNOISE)#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.setParam('obstime','UTMID')     #UT at middle
  iraf.datapars.saveParList(filename=DataParsFileName)
else:
  iraf.datapars.setParList(ParList=DataParsFileName)

########################################################################
########################################################DAOFIND
#Set up daofind to go without prompting for input
iraf.daofind.setParam('image',FitsFileName)             #Set ImageName
iraf.daofind.setParam('verify','no')                    #Don't verify
iraf.daofind.setParam('interactive','no')               #Interactive
iraf.daofind.saveParList(filename=HxRG.PhotDir+HxRG.ObjectName+'_'+HxRG.Filter+\
     '_daofind.par')

#If the output file exists, overwrite it
DaoFindFile = HxRG.PhotDir+HxRG.RawObjectKey+'.coo'
if os.path.isfile(DaoFindFile) and OverRide==1:
   os.remove(DaoFindFile);
iraf.daofind.setParam('output',DaoFindFile)  #Output File

#Run DAOFIND once unless OverRide = 1
#iraf.daofind.setParam('graphics', 'stdgraph')
#iraf.daofind.setParam('display', 'display(image=\'$1\',frame=$2)')
print 'Running DAOFIND'
daostring=iraf.daofind(mode='h',Stdout=1)

#Mark the centroids on the image if TvFlag = 1
if TvFlag >= 1:
  MaxDisp = Im.mean()+3*Im.std()
  MaxDisp = 20
  iraf.set(stdimage=HxRG.StdImage)
  iraf.display(FitsFileName,1,z1=0,z2=MaxDisp,zrange='no')
  if os.path.exists(DaoFindFile):
     iraf.tvmark(1, DaoFindFile)
     iraf.imexamine()