##^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009
#! /usr/bin/env python
#
##This script will run through a given fits file trying various values of 
##C_threshold using DAOFIND
#
#Example Calling Sequence:
#run ~/python/H4RGAnly/daofind_loop_fwhmpsf '/nfs/slac/g/ki/ki09/lances/Reduced/07Apr26/M13/M13_G_SlopeFlatFieldedCenteredDitheredMean.fits' 1
import sys,os
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]
TvFlag	     = sys.argv[2]

#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/keywords.py'):
  execfile('/afs/slac/u/ki/lances/python/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()

#Paramters to use
FWHM = 11.8
print datapars[0]
#CENTERPARS
iraf.centerpars.setParam('cbox',4.5*FWHM)
iraf.centerpars.setParam('cthreshold','0')
iraf.centerpars.saveParList(filename=OutDir+'datcentes.par')

#FINDPARS
iraf.findpars.setParam('threshold',2.5)
iraf.findpars.setParam('sharphi',.6)

#DATAPARS
#Good description at 
#http://stsdas.stsci.edu/cgi-bin/gethelp.cgi?datapars.hlp
iraf.datapars.setParam('fwhmpsf',FWHM)  #FWHM of biggest star
iraf.datapars.setParam('sigma','0.31') #Std. Dev of sky pixels
iraf.datapars.setParam('datamin','1')   #Min good data value ~Sky-3Sigma
iraf.datapars.setParam('datamax','600') #Max good data value
iraf.datapars.setParam('epadu','1')
iraf.datapars.setParam('readnoise','0.3')
iraf.datapars.saveParList(filename=OutDir+'datdataps.par')

#DAOPARS
iraf.daopars.setParam('psfrad',4*FWHM)
iraf.daopars.setParam('fitrad',1.1*FWHM)
iraf.daopars.saveParList(filename=OutDir+'datdaos.par')

#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.saveParList(filename=OutDir+'daofind.par') #Save values

#If the output file exists, overwrite it
DaoFindFile = OutDir+BaseKey+'.coo'
if os.path.isfile(DaoFindFile):
   os.remove(DaoFindFile);
iraf.daofind.setParam('output',DaoFindFile)  #Output File

#Run DAOFIND once
daostring=iraf.daofind(mode='h',Stdout=1)    

#Mark the centroids on the image if TvFlag = 1
if int(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)

'''
imax=20
num_points=zeros(imax,Int32)
threshold_arr=zeros(imax,Int32)

for i in range(imax):
	fwhm=float(i+10)
	iraf.datapars.setParam('fwhmpsf',str(fwhm))
	daostring=iraf.daofind(mode='h',Stdout=1)		#Run hidden
	lendaostring=len(daostring)
	print lendaostring
	num_points[i]=lendaostring
	threshold_arr[i]=fwhm

	if i==0: first_string=daostring[1].split()[1]
	second_string=daostring[1]
	third_string=daostring[2]
	splitdao1=daostring[1].split()[1]	
	#header=splitdao1.split()		#[0] is empty [1] contains header
	splitdao2=daostring[3].split()		#[2] is empty [3]-end contain data
#	print first_string
#	print second_string
#	print third_string
	print splitdao1
	print splitdao2[0]
	print splitdao2[1]

plot(threshold_arr,num_points)
title('For image: '+first_string)
xlabel('Threshold (Sigma)')
ylabel('Number of Points Detected')
savefig('~/Test'+Month+'_'+Day+'_fwhm'+
'_'+first_string+'_Num_spots.ps')
show()
'''