##^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009
#! /usr/bin/env python
#
#
#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     = ""- String describing the conditions used for the psf fit
# 
#Example Calling Sequence:
#run ~/python/H4RGAnly/color_mag_dg '/nfs/slac/g/ki/ki09/lances/Reduced/07Apr26/M13/M13_G_SlopeFlatFieldedCenteredDitheredMean.fits' '/nfs/slac/g/ki/ki09/lances/Reduced/07Apr26/M13/M13_I_SlopeFlatFieldedCenteredDitheredMean.fits' --TvFlag=1 
import sys,os,getopt
from   pyraf import iraf 
import pyfits
from   pylab import *
import numpy
from   numpy import *
import numdisplay
import matplotlib.pylab as mplot

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

#Now parse the command line keyword arguments 
keywords = ['TvFlag=', 'SkipDAOFind=', 'SkipDAOPhot=', 'SkipPSF=', \
	    'SkipNStar=', 'SkipSubStar=','SkipAllStar=','OverRide=', \
            'SkipPSF2=', 'OverWrite=', 'Comment=' ]
OverWrite   = 0
OverRide    = 0
TvFlag      = 0
SkipDAOFind = 0
SkipDAOPhot = 0
SkipPSF     = 0
SkipPSF2    = 0
SkipNStar   = 0
SkipSubStar = 0
SkipAllStar = 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 ['--Comment']:
    Comment   = p
#Extract various keys from the full path
GBaseName     = os.path.basename(GFitsFileName)
DirName      = os.path.dirname(GFitsFileName)
GBaseParts    = GBaseName.split('_')
GObjectName   = GBaseParts[0]
GThisFilter   = GBaseParts[1] 
GBaseParts    = GBaseName.split('.')
GBaseKey      = GBaseParts[0]

IBaseName     = os.path.basename(IFitsFileName)
DirName      = os.path.dirname(IFitsFileName)
IBaseParts    = IBaseName.split('_')
IObjectName   = IBaseParts[0]
IThisFilter   = IBaseParts[1]
IBaseParts    = IBaseName.split('.')
IBaseKey      = IBaseParts[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
GIm = pyfits.getdata(GFitsFileName)
IIm = pyfits.getdata(GFitsFileName)

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

#KEY PARAMETERS Paramters to use
FWHM     = 10.5		  #Full Width at Half Max 
GAIN     = 5.0 		  #e-/ADU
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+GObjectName+'_'+GThisFilter+'PhotPar.par'
iraf.photpars.setParList(ParList = PhotParsFileName)

#CENTERPARS #############################################################
CenterParsFileName = OutDir+GObjectName+'_'+GThisFilter+'CenterPar.par'
iraf.centerpars.setParList(ParList=CenterParsFileName)

#SKYPARS ################################################################
FitSkyParsFileName = OutDir+GObjectName+'_'+GThisFilter+'FitSkyPar.par'
iraf.fitskypars.setParList(ParList=FitSkyParsFileName)

#FINDPARS ##############################################################
FindParsFileName =  OutDir+GObjectName+'_'+GThisFilter+'FindPar.par'
iraf.findpars.setParList(ParList = FindParsFileName)

#DATAPARS - http://stsdas.stsci.edu/cgi-bin/gethelp.cgi?datapars.hlp
DataParsFileName =  OutDir+GObjectName+'_'+GThisFilter+'DataPar.par'
iraf.datapars.setParList(ParList=DataParsFileName)

#DAOPARS ###############################################################
DAOParsFileName =  OutDir+GObjectName+'_'+GThisFilter+'DAOPar.par'
iraf.daopars.setParList(ParList = DAOParsFileName)

#ALL OF THE FILES THAT WERE GENERATED BY DAOPHOT.PY
GDaoFindFile = OutDir+GBaseKey+'.coo'
GDaoPhotFile = OutDir+GBaseKey+'.mag'
GFStList     = OutDir+GBaseKey+'.pst'
GPSFGrList   = OutDir+GBaseKey+'.psg'
GPSFImHFile  = OutDir+GBaseKey+'.psf.fits'
GNStarFile   = OutDir+GBaseKey+'.nst'
GNRejFile    = OutDir+GBaseKey+'.nrj'
GAllStarFile = OutDir+GBaseKey+'.als'
GRejFile     = OutDir+GBaseKey+'.arj'
GAllSubFile  = OutDir+GBaseKey+'All.sub.fits'
GSubIm = pyfits.getdata(GAllSubFile)

#ALL OF THE FILES THAT WERE GENERATED BY DAOPHOT.PY
IDaoFindFile = OutDir+IBaseKey+'.coo'
IDaoPhotFile = OutDir+IBaseKey+'.mag'
IFStList     = OutDir+IBaseKey+'.pst'
IPSFGrList   = OutDir+IBaseKey+'.psg'
IPSFImHFile  = OutDir+IBaseKey+'.psf.fits'
INStarFile   = OutDir+IBaseKey+'.nst'
INRejFile    = OutDir+IBaseKey+'.nrj'
IAllStarFile = OutDir+IBaseKey+'.als'
IRejFile     = OutDir+IBaseKey+'.arj'
IAllSubFile  = OutDir+IBaseKey+'All.sub.fits'
ISubIm = pyfits.getdata(IAllSubFile)

#Get the magnitudes, magnitude errors, CHI and SHARP from the .als file
#G Image Data ########################
GAllStarComps = iraf.pdump(GAllStarFile,\
                'xcenter, ycenter, mag, merr, chi, sharpness, msky', \
                'mag != INDEF', Stdout = 1)
GAllStarLen = len(GAllStarComps)
print 'There are : '+str(GAllStarLen) + ' stars in G'
GXCen   = zeros(GAllStarLen, dtype=float64)
GYCen   = zeros(GAllStarLen, dtype=float64)
GMags   = zeros(GAllStarLen, dtype=float64)
GMErrs  = zeros(GAllStarLen, dtype=float64) 
GChi    = zeros(GAllStarLen, dtype=float64)
GSharp  = zeros(GAllStarLen, dtype=float64)
GMSky   = zeros(GAllStarLen, dtype=float64)

for i, line in zip(range(GAllStarLen),GAllStarComps) : 
    Pars = line.split()
    GXCen[i]  = float(Pars[0])
    GYCen[i]  = float(Pars[1])
    GMags[i]  = float(Pars[2])
    GMErrs[i] = float(Pars[3])
    GChi[i]   = float(Pars[4])
    GSharp[i] = float(Pars[5])
    GMSky[i]  = float(Pars[6])

#I Image Data #####################################################
IAllStarComps = iraf.pdump(IAllStarFile, \
		'xcenter, ycenter, mag, merr, chi, sharpness, msky', \
                'mag != INDEF', Stdout = 1)
IAllStarLen = len(IAllStarComps)
print 'There are: '+str(IAllStarLen)+ ' stars in I'
IXCen   = zeros(IAllStarLen, dtype=float64)
IYCen   = zeros(IAllStarLen, dtype=float64)
IMags   = zeros(IAllStarLen, dtype=float64)
IMErrs  = zeros(IAllStarLen, dtype=float64)
IChi    = zeros(IAllStarLen, dtype=float64)
ISharp  = zeros(IAllStarLen, dtype=float64)
IMSky   = zeros(IAllStarLen, dtype=float64)
 
#Attempt to match the magnitudes
GDupMags  = [0]
IDupMags  = [0]
GIDupMags = [0]

#Get the offsets between the files
iraf.imcentroid.setParam('input',GFitsFileName)
iraf.imcentroid.setParam('reference',IFitsFileName)
iraf.imcentroid.setParam('coords',IDaoPhotFile)
GIXOffset = 10	#I X Coord - GIXOffset = I Y Coord
GIYOffset = 0   #I Y Coord - GIYOffset = G Y Coord
MatchNum  = 0
for i, line in zip(range(IAllStarLen),IAllStarComps) :
    Pars = line.split()
    IXCen[i]  = float(Pars[0])
    IYCen[i]  = float(Pars[1])
    IMags[i]  = float(Pars[2])
    IMErrs[i] = float(Pars[3])
    IChi[i]   = float(Pars[4])
    ISharp[i] = float(Pars[5])
    IMSky[i]  = float(Pars[6])

    #Try to match the array
    XCenMatch = (where(abs((IXCen[i]-GIXOffset)-GXCen) < 4))[0] 
    YCenMatch = (where(abs((IYCen[i]-GIYOffset)-GYCen[XCenMatch]) < 4))[0]
    if len(YCenMatch) == 1 :
       MatchInd = XCenMatch[YCenMatch]
       print 'Y Index'+str(MatchInd) 
       print 'Found Match : '+str(MatchNum) + ' at star number ' +str(i)
       print 'XCen G: ' + str(GXCen[MatchInd]) + 'XCen I: '+str(IXCen[i])  
       print 'YCen G: ' + str(GYCen[MatchInd]) + 'YCen I: '+str(IYCen[i])
       MatchNum +=1
       GDupMags=append(GDupMags,GMags[MatchInd])
       IDupMags=append(IDupMags,IMags[i])
       GIDupMags=append(GIDupMags,GMags[MatchInd]-IMags[i])

#Plot the magnitudes vs. error for G
mplot.clf()
mplot.semilogy(GMags,GMErrs,'b+')
mplot.xlabel('Magnitude')
mplot.ylabel('Magnitude Error')
mplot.title('Magnitudes for '+GThisFilter+' Band')
mplot.savefig(OutDir+'GMagVsMagErr.png')

#Plot the magnitudes vs. error for I
mplot.cla()
mplot.semilogy(IMags,IMErrs,'r+')
mplot.xlabel('Magnitude')
mplot.ylabel('Magnitude Error')
mplot.title('Magnitudes for '+IThisFilter+' Band')
mplot.savefig(OutDir+'IMagVsMagErr.png')

#Plot the magnitudes vs. error for I
mplot.cla()
mplot.plot(IChi,IMags,'r+')
mplot.xlabel('Chi (Goodness of fit)')
mplot.ylabel('Magnitude')
mplot.title('Magnitudes for '+IThisFilter+' Band')
mplot.savefig(OutDir+'IMagVsChi.png')

mplot.cla()
mplot.plot(GChi,GMags,'b+')
mplot.xlabel('Chi (Goodness of fit)')
mplot.ylabel('Magnitude')
mplot.title('Magnitudes for '+GThisFilter+' Band')
mplot.savefig(OutDir+'GMagVsChi.png')

#Plot the median sky vs. magnitude
mplot.cla()
mplot.plot(GMSky ,GMags,'b+')
mplot.xlabel('Median Sky')
mplot.ylabel('Magnitude')
mplot.title('Magnitudes for '+GThisFilter+' Band')
mplot.savefig(OutDir+'GMagVsSky.png')

mplot.cla()
mplot.plot(IMSky ,IMags,'r+')
mplot.xlabel('Median Sky')
mplot.ylabel('Magnitude')
mplot.title('Magnitudes for '+GThisFilter+' Band')
mplot.savefig(OutDir+'IMagVsSky.png')

#Plot the Color Magnitude Diagram - 0th element is zero
GIDupMags = GIDupMags[1:]
IDupMags  = IDupMags[1:]
GDupMags  = GDupMags[1:]
mplot.clf()
mplot.axis([GIDupMags.min(), GIDupMags.max(),\
		       IDupMags.max(), IDupMags.min()])
mplot.plot(GIDupMags, IDupMags, 'b+')
mplot.title('Color Magnitude Diagram for M13 : No Cuts Implemented')
mplot.xlabel('G-I')
mplot.ylabel('Magnitude in I')
mplot.figtext(.15,.85, str(GAllStarLen) + '    stars in G')
mplot.figtext(.15,.81, str(IAllStarLen) + '  stars in I')
mplot.figtext(.15,.77, str(MatchNum)    + '    matched stars')
mplot.savefig(OutDir+'M13ColorMagDiag.png')

#Print some values
print 'There are : '+str(GAllStarLen) + ' stars in G'
print 'There are : '+str(IAllStarLen) + ' stars in I'
print 'There are : '+str(MatchNum) + ' stars matched in G and I'