##^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009
#DAOFIND_HxRG
#
#PURPOSE:
#  To fit PSFs to a list of stars generated by DAOFind
#  1) Run DAOFIND_HxRG.py to form a list of coordinates
#  2) Run this script to fit radial profiles
#
#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 RadProf_HxRG.py '/nfs/slac/g/ki/ki03/lances/H2RG-32-147/ASIC/Reduced/07Dec14/SAO81134/SAO81134_H2RG_SIPIN_220_Reads_Dec15_2007_04_57_19_I_RPS_Slope.fits' --TvFlag=0

##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'

#Run the RadProf module
RadProfFile    = HxRG.PhotDir+HxRG.RawObjectKey+'.prf'
RadProfTxtFile = HxRG.PhotDir+HxRG.RawObjectKey+'.prf.txt'

if os.path.isfile(RadProfFile) and OverWrite==1:
   os.remove(RadProfFile);
iraf.daofind.setParam('output',DaoFindFile)  #Output File
iraf.digiphot()
iraf.digiphot.apphot()
iraf.radprof.setParam('image',FitsFileName)
iraf.radprof.setParam('radius',20)
iraf.radprof.setParam('step',0.5)
iraf.radprof.setParam('coords', DaoFindFile)
iraf.radprof.setParam('output', RadProfFile)
iraf.radprof.setParam('verify', 'no')
iraf.radprof.setParam('interactive','yes')
iraf.radprof(mode='h')

#Dump the relevant quantities to a txt file
print 'Extracting Parameters to file : ' + RadProfTxtFile
if os.path.exists(RadProfTxtFile) and OverWrite==1: 
  os.remove(RadProfTxtFile)
 
iraf.txdump(RadProfFile,\
  'XCENTER,YCENTER,PFWHM,INORM,MSKY,TINORM', Stdout=RadProfTxtFile, mode='h')