##^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009
#! /usr/bin/env python

##This script will use the immatchx.ccmap utility in IRAF to 
##compute the WCS FITS Keywords based upon the pixel and sky
##coordinates of several bright stars in an image
##
##InputCooList should have the following columns
##	RA	DEC	X	Y
##
#CALLING SEQUENCE:
#run ~/python/pyraf/CorrectWCS.py '/nfs/slac/g/ki/ki04/lances/H4RG-10-007/Reduced/07Apr26/M13/M13_G_SlopeFlatFieldedCenteredDitheredMeanWCS.fits' '/nfs/slac/g/ki/ki04/lances/H4RG-10-007/Reduced/07Apr26/M13/Astrometry/M13_G_CoordinatesUBV.txt' '/nfs/slac/g/ki/ki04/lances/H4RG-10-007/Reduced/07Apr26/M13/Astrometry/M13UBVPhotometryAndProperMotionsKadla1966.txt' --TvFlag=1 
import sys,os, getopt, string
from pyraf import iraf
import pyfits
from numpy import *
from pylab import *

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

opts, extraparams = getopt.getopt(sys.argv[4:],'',keywords)
for o,p in opts:
  if o in ['--TvFlag']:
    TvFlag=int(p)

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

#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 Astrometry output
OutDir = DirName+'/Astrometry/'
if not os.path.exists(OutDir): os.mkdir(OutDir)

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

#InputCooList should have the following columns
SolutionsList = OutDir+ObjectName+'_'+ThisFilter+'_'+'Solutions.txt'
iraf.ccmap.setParam('images',FitsFileName)
iraf.ccmap.setParam('input', CoordFileName)
iraf.ccmap.setParam('database',SolutionsList)
iraf.ccmap.setParam('lngcolumn',1)
iraf.ccmap.setParam('latcolumn',2)
iraf.ccmap.setParam('xcolumn',3)
iraf.ccmap.setParam('ycolumn',4)
iraf.ccmap.setParam('lngunits','degrees')
iraf.ccmap.setParam('update','yes')
iraf.ccmap(interactive='yes')

#Now display them on an image if so desired
if TvFlag > 0:
  Im = pyfits.getdata(FitsFileName)
  MaxDisp = Im.mean()+3*Im.std()
  iraf.set(stdimage='imt4096')
  iraf.display(FitsFileName,1,z1=1,z2=MaxDisp,zrange='no')

  iraf.tproject.setParam('intable', CoordFileName)
  iraf.tproject.setParam('outtable','STDOUT')
  iraf.tproject.setParam('columns','c3,c4')
  XYCoordList = iraf.tproject(Stdout=1)
  iraf.tvmark.setParam('coords',"STDIN")
  iraf.tvmark(1,Stdin=XYCoordList,col=205)

  iraf.tproject.setParam('intable', FullCoordFileName)
  iraf.tproject.setParam('outtable','STDOUT')
  iraf.tproject.setParam('columns', 'c11,12')
  FullRACoordList = iraf.tproject(Stdout=1)
  iraf.wcsctran.setParam('input', 'STDIN')
  iraf.wcsctran.setParam('output','STDOUT')
  iraf.wcsctran.setParam('image',FitsFileName)
  iraf.wcsctran.setParam('inwcs','world')
  iraf.wcsctran.setParam('outwcs','tv')
  FullXYCoordList = iraf.wcsctran(Stdin=FullRACoordList,Stdout=1)
  iraf.tvmark(1,Stdin=FullXYCoordList,col=204)