##^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009
#! /usr/bin/env python
#
#This script will run DAOPHOT routines on a given fits file (FitsFileName) 
#specified on the command line. The structure of the filename is specific to 
#observations taken with the H4RG detector at the KPNO 2.1m.   
#
#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/DAOPhot.py '/nfs/slac/g/ki/ki09/lances/Reduced/07Apr27/Landolt110956/Landolt110956_Y_SlopeFlatFieldedCenteredDitheredMean.fits' --TvFlag=2 --SkipDAOFind=0 --SkipDAOPhot=0 --SkipPSF=1 --OverRide=1 --SkipPSF2=1 --SkipAllStar=1 --SkipNStar=1 SkipSubStar=1 --OverWrite=1 --Comment='' --UsePSFList=0
#
#run ~/python/H4RGAnly/DAOPhot.py '/nfs/slac/g/ki/ki09/lances/Reduced/07Apr28/Landolt/Landolt_I_SlopeFlatFieldedCenteredDitheredMean.fits' --TvFlag=2 --SkipDAOFind=0 --SkipDAOPhot=0 --SkipPSF=1 --OverRide=1 --SkipPSF2=1 --SkipAllStar=1 --SkipNStar=1 SkipSubStar=1 --OverWrite=1 --Comment='' --UsePSFList=0
#
#run ~/python/H4RGAnly/DAOPhot.py '/nfs/slac/g/ki/ki09/lances/Reduced/07Apr26/M13/M13_G_SlopeFlatFieldedCenteredDitheredMean.fits' --TvFlag=2 --SkipDAOFind=0 --SkipDAOPhot=0 --SkipPSF=0 --OverRide=1 --SkipPSF2=0 --SkipAllStar=0 --SkipNStar=0 --SkipSubStar=0 --OverWrite=1 --Comment='VarPSF' --UsePSFList=0

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
FitsFileName = 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
SkipPSF     = 0
SkipPSF2    = 0
SkipNStar   = 0
SkipSubStar = 0
SkipAllStar = 0
UsePSFList  = 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 ['--SkipPSF']:
    SkipPSF = int(p)
  elif o in ['--SkipNStar']:
    SkipNStar = int(p)
  elif o in ['--SkipSubStar']:
    SkipSubStar = int(p)
  elif o in ['--SkipAllStar']:
    SkipAllStar = int (p)
  elif o in ['--SkipPSF2']:
    SkipPSF2 = int(p)
  elif o in ['--OverWrite']:
    OverWrite = int(p)
  elif o in ['--UsePSFList']:
    UsePSFList = int(p)
  elif o in ['--Comment']:
    Comment = p

#Extract various keys from the full path
BaseName     = os.path.basename(FitsFileName)
DirName      = os.path.dirname(FitsFileName)
BaseParts    = BaseName.split('_')
ObjectName   = BaseParts[0]
ThisFilter   = BaseParts[1] 
BaseParts    = BaseName.split('.')
BaseKey      = BaseParts[0]

#Create a directory for Photometry output
OutDir = DirName+'/Photometry/'
if not os.path.exists(OutDir): os.mkdir(OutDir)

#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()

#Read in the array
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()
skypars    = iraf.fitskypars.getParList()
photpars   = iraf.photpars.getParList()

#KEY PARAMETERS Paramters to use
FWHM      = 15		  #Full Width at Half Max 
GAIN      = 20.0#5.0 		  #e-/ADU
READNOISE = 0.45		  #Readnoise for a single
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) 

#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',27.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',12.0)
  iraf.findpars.setParam('sharphi',  1)
  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','1.8')   #Std. Dev of sky pixels
  iraf.datapars.setParam('datamin','27.1')  #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.setParam('obstime','UTMID')	#UT at middle
  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('varorder',2)	  #Try using variable PSF
  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('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',FitsFileName)		#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')

#If the output file exists, overwrite it
DaoFindFile = OutDir+BaseKey+Comment+'.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)')
if OverRide == 1 and SkipDAOFind != 1: 
   print 'Running DAOFIND'
   daostring=iraf.daofind(mode='h',Stdout=1)
else :
   print 'Skipping DAOFind'
   print 'Using : '+DaoFindFile

#Mark the centroids on the image if TvFlag = 1
if TvFlag >= 1:
  MaxDisp = Im.mean()+3*Im.std()
  iraf.set(stdimage='imt4096')
  iraf.display(FitsFileName,1,z1=1,z2=MaxDisp,zrange='no')
  if os.path.exists(DaoFindFile):
     TvMarkFile = OutDir+BaseKey+'TvMark.fits' 
     if os.path.exists(TvMarkFile): os.remove(TvMarkFile)
     iraf.tvmark(1, DaoFindFile, outimage=TvMarkFile)

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

if OverRide ==1 and SkipDAOPhot !=1 : 
   print 'Running DAOPhot'
   if os.path.exists(DaoPhotFile) : os.remove(DaoPhotFile)
   photstring = iraf.phot(mode='h',Stdout=1)
   DaoPhotTxtFile=OutDir+BaseKey+Comment+'.mag.txt'
   print 'Extracting Parameters to file : ' + DaoPhotTxtFile
   if os.path.exists(DaoPhotTxtFile) : os.remove(DaoPhotTxtFile)
   iraf.txdump(DaoPhotFile,\
   'XCENTER,YCENTER,XERR,YERR,CERROR,MSKY,FLUX,MAG,MERR','yes',\
   Stdout = DaoPhotTxtFile)
else :
   print 'Skipping DAOPhot'
   print 'Using : ' + DaoPhotFile

#Check to see what the bad fits were
if TvFlag >= 1 :
   BadPts  = iraf.pdump(DaoPhotFile,'XCENTER,YCENTER','mag==INDEF', Stdout = 1) 
   iraf.tvmark.setParam('coords','STDIN')
   iraf.tvmark(1,Stdin = BadPts, mode='h', col=204)
   GoodPts = iraf.pdump(DaoPhotFile,'XCENTER,YCENTER','mag!=INDEF', Stdout = 1)
   iraf.tvmark(1,Stdin = GoodPts, mode='h',col=205)
   if TvFlag == 2 : iraf.imexamine()

########################################################PSF
PSFStList = OutDir+BaseKey+Comment+'.pst'
PSFGrList = OutDir+BaseKey+Comment+'.psg'
PSFImHFile= OutDir+BaseKey+Comment+'.psf.fits'
if OverRide ==1 and SkipPSF !=1 :
  print 'Running PSF'
  #Set defaults, but allow use of a psf list
  iraf.psf.setParam('image',FitsFileName)
  iraf.psf.setParam('graphics', 'stdgraph')
  iraf.psf.setParam('display', 'display(image=\'$1\',frame=$2)')
  iraf.psf.setParam('verify','yes')
  iraf.psf.setParam('photfile',DaoPhotFile)
  iraf.psf.setParam('psfimage',PSFImHFile)
  iraf.psf.setParam('groupfile',PSFGrList)
  iraf.psf.setParam('opstfile',PSFStList) 

  if os.path.exists(PSFImHFile) and SkipPSF !=1 : os.remove(PSFImHFile)
  if os.path.exists(PSFGrList)  and SkipPSF !=1 : os.remove(PSFGrList)
  #Use a PSF list or allow the interactive choice
  if UsePSFList == 1 :
     iraf.psf.setParam('opstfile',"STDOUT")
     iraf.psf.setParam('pstfile',PSFStList)
     iraf.psf.setParam('interactive','no')
  else:
     if os.path.exists(PSFStList)  and SkipPSF !=1 : os.remove(PSFStList)
     iraf.psf.setParam('pstfile','')
     iraf.psf.setParam('interactive','yes')
  #Save the Parameters to the same file
  iraf.psf.saveParList(filename=OutDir+ObjectName+\
		'_'+ThisFilter+'PSFPar.par')
  PSFString = iraf.psf()

else :
  print 'Skipping PSF Fit'

########################################################NSTAR
#Set PSFRad to smaller value in order to fit simultaneously
NStarFile = OutDir+BaseKey+Comment+'.nst'
NRejFile  = OutDir+BaseKey+Comment+'.nrj'
iraf.daopars.setParam('recenter','yes')
iraf.daopars.setParam('maxgroup', 100)
iraf.daopars.setParam('psfrad',3*FWHM)
if OverRide == 1 and SkipNStar !=1:
  if os.path.exists(NStarFile) : os.remove(NStarFile)
  if os.path.exists(NRejFile)  : os.remove(NRejFile)
  iraf.nstar.setParam('image',FitsFileName)
  iraf.nstar.setParam('groupfile', PSFGrList)
  iraf.nstar.setParam('psfimage', PSFImHFile)
  iraf.nstar.setParam('nstarfile',NStarFile)
  iraf.nstar.setParam('rejfile', NRejFile)
  iraf.nstar.setParam('verify','yes')
  iraf.nstar.saveParList(filename=OutDir+ObjectName+\
                '_'+ThisFilter+'NStarPar.par')
  NStarStr=iraf.nstar()#mode='h')

########################################################SUBSTAR
SkipSubStar = 0
SubStarFile = OutDir+BaseKey+Comment+'.sub.fits'
if OverRide == 1 and SkipSubStar !=1 :
  if os.path.exists(SubStarFile) : os.remove(SubStarFile)
  iraf.substar.setParam('image',FitsFileName)
  iraf.substar.setParam('photfile',NStarFile)
  if SkipPSF2 == 1 :
    iraf.substar.setParam('exfile',  NRejFile)
  else :
    iraf.substar.setParam('exfile',  PSFStList)
  iraf.substar.setParam('psfimage',PSFImHFile)
  iraf.substar.setParam('subimage',SubStarFile)
  iraf.substar.setParam('verify','yes')
  SubStarStr = iraf.substar()#mode='h')

#Display the subtracted image if TvFlag == 1
if TvFlag >= 1 : 
  iraf.display(SubStarFile,1,z1=1,z2=30,zrange='no')
  PSFPts = iraf.pdump(PSFGrList,'XCENTER,YCENTER','mag != INDEF', Stdout = 1)
  iraf.tvmark(1,Stdin = PSFPts, mode='h', col=205)
  PSFBadPts = iraf.pdump(PSFGrList,'XCENTER,YCENTER','mag == INDEF',Stdout = 1)
  iraf.tvmark(1,Stdin = PSFBadPts, mode='h', col=204)
  if TvFlag == 2 : iraf.imexamine()

########################################################PSF 2ND TIME
iraf.daopars.setParam('psfrad',4*FWHM)
if OverRide ==1 and SkipPSF2 !=1 :
  PSFStList = OutDir+BaseKey+'.pst'
  PSFGrList = OutDir+BaseKey+'.psg'
  PSFImHFile= OutDir+BaseKey+'.psf.fits'
  print 'Running Second PSF'
  iraf.psf.setParam('image',SubStarFile)
  if os.path.exists(PSFImHFile) and SkipPSF !=1 : os.remove(PSFImHFile)
  if os.path.exists(PSFGrList)  and SkipPSF !=1 : os.remove(PSFGrList)
  #Use a PSF list or allow the interactive choice
  if UsePSFList == 1 :
     iraf.psf.setParam('opstfile',"STDOUT")
     iraf.psf.setParam('pstfile',PSFStList)
     iraf.psf.setParam('interactive','no')
  else:
     if os.path.exists(PSFStList)  and SkipPSF !=1 : os.remove(PSFStList)
     iraf.psf.setParam('pstfile','')
     iraf.psf.setParam('interactive','yes')
  #Save the Parameters to the same file
  iraf.psf.saveParList(filename=OutDir+ObjectName+\
                '_'+ThisFilter+'PSFPar.par')
  PSFString = iraf.psf()

########################################################ALLSTAR
AllStarFile = OutDir+BaseKey+Comment+'.als'
RejFile     = OutDir+BaseKey+Comment+'.arj'
AllSubFile  = OutDir+BaseKey+Comment+'All.sub.fits'
iraf.daopars.setParam('psfrad',4*FWHM)

if OverRide == 1 and SkipAllStar !=1 :
  if os.path.exists(AllStarFile) : os.remove(AllStarFile)
  if os.path.exists(RejFile) : os.remove(RejFile)
  if os.path.exists(AllSubFile) : os.remove(AllSubFile)
  iraf.allstar.setParam('image',FitsFileName)
  iraf.allstar.setParam('photfile',DaoPhotFile) 
  iraf.allstar.setParam('psfimage',PSFImHFile)
  iraf.allstar.setParam('allstarfile',AllStarFile)
  iraf.allstar.setParam('rejfile',RejFile)
  iraf.allstar.setParam('subimage',AllSubFile)
  AllStarStr = iraf.allstar()

#Open up the subtracted Image and show what the psfs look like
SubIm = pyfits.getdata(AllSubFile)
PSFIm = Im - SubIm

#Make the file nae
PSFImName = OutDir+BaseKey+Comment+'PSFIm.fits'
if os.path.exists(PSFImName) : os.remove(PSFImName)
hdu = pyfits.PrimaryHDU(PSFIm)
hdu.writeto(PSFImName)

if TvFlag >= 1 :
   iraf.display(AllSubFile,1,z1=-4,z2=30,zrange='no')