#! /usr/bin/env python
#
# Paul Martini (OSU) 
#
#  oalign.py [omsFile] maskImage [fieldImage]
#
# if [fieldImage] not specified, perform alignment through maskImage. 
# [omsFile] should not specified in longslit mode 
#
# Steps: 
#  1. identify alignment boxes and mark them on the mask image
#     - interactively adjust positions 
#  2. predict x,y values for alignment stars and plot on image
#     - interactively adjust positions
#  3. calculate translations and rotation
#
#  18 Feb 2011: implemented recenter of box coords 
#  20 Feb 2011: added bells and whistles, such as plotting of residuals
#  21 Feb 2011: streamlined code further, such as command line flags 
#  22 Feb 2011: added geomap for calculation of transformation 
#  18 Mar 2011: improved performance when aligning with stars in boxes already
#  2011-03-20 [0.3.1]: tested with MODS data, rotates around correct center, 
#       further improved performance when bright stars are in alignment boxes
#  2011-04-15: various cuts at new prompts, making error messages similar to
#              other MODS scripts, and notes for possible chnages marked
#              with ?? [rwp]
#  2011-05-12: modified from mosalign.py 
#  2011-05-13: de-sophisticated to run on hiltner, agung 
#  2011-05-16: added longslit option 
#  2011-09-12: fixed biasc on line 411, started to implement direct 
#	       rotator communication with Ray Gonzalez's code 
#
#---------------------------------------------------------------------------

from string import *
from time import sleep
# from subprocess import call
import socket
import sys
import os 
import math 
from sys import argv, exit
from numpy import mean, std, zeros, median 
# from numpy.linalg import lstsq
from pyraf import iraf 
import pyfits 
iraf.images.tv() 
iraf.images.immatch() 
iraf.set(stdimage='imt4096')

versNum = "0.9.0"
versDate = "2011-05-16"

# Rotation center (initial guess) 
xrotc = 2000
yrotc = 2000

# Defaults for various flags
verbose=1		# Verbose mode
shiftonly=0		# Only compute xy shift (no rotation) 
boxdelete=1		# Re-selects alignment star boxes if set to 1
debug=0			# Unadvertised debug mode (lots of output) 	
version=0		# Print version information 
longslit=0		# Long slit acquisition mode [default is MOS] 

############################
#### Define various routines
############################

def usage(): 
  print "\nMulti-Slit Usage: \n" 
  print "	%s -flags omsFile maskImage [acqImage]" % (argv[0])
  print "\nWhere: omsFile is an OSMOS mask design file" 
  print "       maskImage is an image through the mask" 
  print "       acqImage is an image of the target field (no mask)" 
  print "       Note: Do not specify fieldImage if maskImage has stars" 
  print "             well inside the alignment boxes." 
  print "\nLongslit Usage: (specify -l) \n" 
  print "	%s -flags slitImage acqImage" % (argv[0])
  print "\nWhere: slitImage is an image of the longslit" 
  print "       acqImage is an image of the target field (no slit)" 
  print "\nOptions: -v print verbose output" 
  print "         -s compute the x,y shift only" 
  print "         -b use previously-selected alignment boxes"
  print "         -V print version information" 
  print "         -l longslit acquisition (multi-slit is default)\n" 
#
# old sense of -b
#  print "         -b force re-selection of alignment boxes"
# 

def parseinput():
  flags = []
  files = []
  # check for any command-line arguments and input files 
  for i in range(1,len(argv)): 
    if find(argv[i], "-") == 0: 
      flags.append(argv[i].strip("-"))
    else: 
      files.append(argv[i]) 
  # check that the input files exist
  for i in range(1,len(files)): 
    if os.path.isfile(files[i]) == 0:
      print "\n** ERROR: "+files[i]+" does not exist." 
      exit(1) 
  return files, flags

def parseoms(file, verbose=0): 
  x = []
  y = []
  dx = []
  dy = []
  F=open(file, "r")
  M=F.readlines()[::]
  F.close()
  nrefslits = 0 
  if verbose: 
    print "Found %d lines in %s" % (len(M), file) 
  for i in range(len(M)): 
    # get the number of reference slits and regular slits
    if find(M[i], "INS.RSLIT.NUMBER") >= 0: 
      nrefslits = int(M[i].split()[1]) 
      # first note reference slit should be 100+nrefslits+1
      firstnonref = "INS.TARG%3d" % (101+nrefslits) 
      if verbose: 
        print "First non-reference slits should be: ", firstnonref
    if find(M[i], "INS.SLIT.NUMBER") >= 0: 
      nslits = int(M[i].split()[1])
    if find(M[i], "XMM") >= 0: 
      x.append(atof(M[i].split()[1])) 
    if find(M[i], "WIDMM") >= 0: 
      dx.append(atof(M[i].split()[1])) 
    if find(M[i], "YMM") >= 0: 
      y.append(atof(M[i].split()[1])) 
    if find(M[i], "LENMM") >= 0: 
      dy.append(atof(M[i].split()[1])) 
  if verbose: 
    print "Found %d reference boxes and %d regular slits" % (nrefslits, nslits)
  return x[nrefslits::], y[nrefslits::], dx[nrefslits::], dy[nrefslits::] 

def start_ds9():
  ds9open = os.system("xpaaccess oalign") 
  # if ds9 display titled 'OSMOS' isn't open, open it
  if ds9open==0:
    if verbose:
      print "Opening ds9 window for oalign..."
    os.spawnlp(os.P_NOWAIT,"ds9", "ds9", "-title", "oalign")
    sleep(2)  # need a little longer sleep, 1 usually caused a race condition [rwp]
    os.spawnlp(os.P_NOWAIT,"xpaset", "xpaset", "-p","oalign","scale","mode","zscale")
  # otherwise just delete the existing region files
  else: 
    os.spawnlp(os.P_WAIT, "xpaset", "xpaset", "-p", "oalign", "frame","2") 
    os.spawnlp(os.P_WAIT, "xpaset", "xpaset", "-p", "oalign", "regions","deleteall")
    os.spawnlp(os.P_WAIT, "xpaset", "xpaset", "-p", "oalign", "frame","1") 
    os.spawnlp(os.P_WAIT, "xpaset", "xpaset", "-p", "oalign", "regions","deleteall")

def end_ds9():
  os.spawnlp(os.P_NOWAIT,"xpaset", "xpaset", "-p","oalign","exit") 

# Should have more graceful error checking here

def get_keyword(file, keyword):
  iraf.imgets(file,keyword) 
  return iraf.imgets.value

def readboxfile(file):  
  F=open(file, "r")
  M=F.readlines()
  V=[split(m) for m in M]
  F.close()
  return list(map(float,zip(*V)[0])), list(map(float,zip(*V)[1])), list(map(int,zip(*V)[2])), list(map(float,zip(*V)[3]))

# Calculate box centers if boxfile does not already exist

def getboxcoords(pixscale, maskimage): 
  # get dimensions of mask image 
  naxis1 = float(get_keyword(maskimage, "naxis1") )
  naxis2 = float( get_keyword(maskimage, "naxis2") )
  binx = float( get_keyword(maskimage, "CCDXBIN") )
  biny = float( get_keyword(maskimage, "CCDYBIN") )
  gprobex = float(get_keyword(maskimage, "GPROBEX") )
  gprobey = float(get_keyword(maskimage, "GPROBEY") )
  binx = 1
  biny = 1
  # identify alignment boxes (just slits with equal width and length) 
  boxx = []
  boxy = []
  index = 1
  G = open("_oalign.reg", "w") 
  for i in range(len(x)): 
    if dx[i] == dy[i]:
      #boxx.append(0.5*naxis1+x[i]*(1.667)/(pixscale*binx)) 
      #boxy.append(0.5*naxis2-y[i]*(1.667)/(pixscale*biny))
      boxx.append(0.5*naxis1+x[i]*(11.528)/(pixscale*binx))
      boxy.append(0.5*naxis2-y[i]*(11.528)/(pixscale*biny))
      regline = "box(%6.2f,%6.2f,40,40) # text={%d} color=black\n" % (boxx[-1], boxy[-1], index)
      index=index+1
      G.write(regline) 
  G.close()
  
  # Open ds9 and draw the predicted locations of the alignment boxes 

  iraf.display(maskimage, 1)
  os.spawnlp(os.P_WAIT, "xpaset", "xpaset", "-p", "oalign", "zoom","to fit")
  # os.spawnlp(os.P_WAIT, "xpaset", "xpaset", "-p", "oalign", "zoom","4") 
  os.spawnlp(os.P_WAIT, "xpaset", "xpaset", "-p", "oalign", "regions","file","_oalign.reg")
 
  ##### Note: automate this in the future ######
  # mark the boxes to get exact x,y coords
  if os.path.isfile("_oalign.log"):
    os.remove("_oalign.log")

  print "\n** Mark Alignment Star Boxes to use with a 'x'."
  print "   and then press 'q' when finished." 
  print "   *** You must mark at least two (2) boxes ***\n" 

  iraf.imexam(maskimage, 1, maskimage, logfile="_oalign.log", keeplog="yes")
  boxx, boxy = readlogfile("_oalign.log") 
  sci = pyfits.open(maskimage) 
  scidata = sci[0].data
  nbox = len(boxx) 
  # boxxc = zeros(nbox, float)        
  # boxyc = zeros(nbox, float)        
  boxxc = []
  boxyc = []
  boxfluxc = []
  # Calculate the median flux at the input coords
  for i in range(nbox): 
    #boxfluxc.append(median(scidata[int(boxy[i])-2:int(boxy[i])+2, int(boxx[i])-2:int(boxx[i])+2]))
    #boxfluxc.append(median(scidata[int(boxy[i])-2:int(boxy[i])+2, int(boxx[i])-2:int(boxx[i])+2], axis=None))
    # fix for old pyfits, etc. on hiltner and agung -- 
    boxfluxc.append(median(scidata[int(boxy[i])-2:int(boxy[i])+2, int(boxx[i])-2:int(boxx[i])+2], axis=None))
    if verbose: 
      print "%d Median at xc,yc is %4.2f" % (i, boxfluxc[i]) 
  if debug: 
    print boxfluxc
    print min(boxfluxc)
  maxboxflux = min(boxfluxc)
  if verbose:
    print "Min for all boxes is %4.2f" % (maxboxflux) 

  # update the centering of the boxes
  
  for i in range(nbox):
    if debug: 
      print boxx[i], boxy[i]
    xc, yc, fluxc = centerbox(scidata, boxx[i], boxy[i], maxboxflux)
    boxxc.append(xc)
    boxyc.append(yc)
    boxfluxc.append(fluxc)
    # print xc, yc, boxxc[i], boxyc[i]
  boxid = [int(i) for i in range(1,len(boxx))]
  
  # Write the alignment box locations to a file
  F = open(boxfile, "w") 
  G = open(boxregfile, "w") 
  for i in range(len(boxx)): 
    # line = "%6.2f %6.2f %3d\n" % (boxx[i], boxy[i], i+1)
    line = "%6.2f %6.2f %3d %6.2f\n" % (boxxc[i], boxyc[i], i+1, boxfluxc[i])
    #regline = "box(%6.2f,%6.2f,40,40) # text={%d} color=black\n" % (boxxc[i], boxyc[i], i+1)
    regline = "box(%6.2f,%6.2f,40,40) # text={%d} color=white\n" % (boxxc[i], boxyc[i], i+1)
    F.write(line) 
    G.write(regline) 
  F.close()
  G.close()
  return boxxc, boxyc, boxid, boxfluxc 

def getslitcoords(pixscale, maskimage): 
  # get dimensions of image 
  naxis1 = float(get_keyword(maskimage, "naxis1") )
  naxis2 = float( get_keyword(maskimage, "naxis2") )
  gprobex = float(get_keyword(maskimage, "GPROBEX") )
  gprobey = float(get_keyword(maskimage, "GPROBEY") )
  binx = float( get_keyword(maskimage, "CCDXBIN") )
  biny = float( get_keyword(maskimage, "CCDYBIN") )
  binx = 1
  biny = 1
  print "\n** Mark where you would like the object on the slit with 'x'."
  print "   and then press 'q' when finished." 

  iraf.imexam(maskimage, 1, maskimage, logfile="_oalign.log", keeplog="yes")
  boxx, boxy = readlogfile("_oalign.log") 
  boxx = boxx[-1]
  boxy = boxy[-1]
  nbox = 1
  sci = pyfits.open(maskimage) 
  scidata = sci[0].data
  boxxc = []
  boxyc = []
  boxfluxc = []
  # Calculate the median flux at the input coords
  if longslit==0: 
    boxfluxc = median(scidata[int(boxy[i])-2:int(boxy[i])+2, int(boxx[i])-2:int(boxx[i])+2], axis=None)
    if verbose: 
      print "Median at xc,yc is %4.2f" % (boxfluxc[i]) 
  else: 
    boxfluxc = median(scidata[int(boxy)-2:int(boxy)+2, int(boxx)-2:int(boxx)+2], axis=None)
    if verbose: 
      print "Median at xc,yc is %4.2f" % (boxfluxc) 
  print boxfluxc
  maxboxflux = boxfluxc
  if verbose:
    print "Min for all boxes is %4.2f" % (maxboxflux) 

  # update the centering of the boxes
  
  if verbose:
    print boxx, boxy
  xc, yc, fluxc = centerslit(scidata, boxx, boxy, maxboxflux)
  boxxc = xc
  boxyc = yc
  boxfluxc = fluxc
  boxid = 1
  
  # Write the alignment box locations to a file
  F = open(boxfile, "w") 
  G = open(boxregfile, "w") 
  line = "%6.2f %6.2f %3d %6.2f\n" % (boxxc, boxyc, 1, boxfluxc)
  regline = "box(%6.2f,%6.2f,10,1000) # text={%d} color=white\n" % (boxxc, boxyc, 1)
  F.write(line) 
  G.write(regline) 
  F.close()
  G.close()
  return boxxc, boxyc, boxid, boxfluxc 

def getdisplaypars(flux): 
  zz1 = 0. 
  if longslit==0: 
    zz2 = max(flux)+3.*std(flux)  
  else: 
    zz2 = flux+500. 
  return zz1, zz2

def getstarcoords(file):
  if longslit == 0: 
    print "\n** Mark the center of each alignment star with an 'a'" 
    print "   and then press 'q' when finished." 
    print "   Stars must be selected in the *SAME ORDER* as the boxes.\n" 
  else: 
    print "\n** Mark the target with an 'a' and then press 'q'\n" 
  sleep(2)
  os.spawnlp(os.P_WAIT, "xpaset", "xpaset", "-p", "oalign", "regions","file", boxregfile) 
  # set use_display=no to force imexam to not redisplay file
  iraf.imexam(file, 2, file, logfile="_oalign.log", keeplog="yes", autoredraw="no", use_display="no")
  starx, stary = readlogfile("_oalign.log") 
  return starx, stary

def readlogfile(file): 
  # read in the reference star locations from the log file
  F = open(file, "r")
  M=F.readlines()
  F.close()
  x = []
  y = []
  # 
  for i in range(len(M)): 
    if find(M[i], "#") <0: 
      x.append(float(M[i].split()[0])) 
      y.append(float(M[i].split()[1])) 

  # check that at least two stars are present - abort if not
  # ?? What other cleanup we should be doing here on abort ??
  
  if len(x)<2 and longslit==0: 
    print "\n** ERROR: Need at least two points to solve for"
    print "   translation and rotation - %s aborting.\n" % (argv[0])
    end_ds9()
    exit(1)
  return x, y

def centerbox(data, x, y, maxflux, N=50, dw=9,verbose=0):
  # N is the search tolerance for box edges in pixels (could get from MMS)
  # dw is the half width of the stripe used for the gradient search
  xw=int(x)  # int of the clicked box center 
  yw=int(y) 
  x1=xw-N    # int of the start of the range 
  x2=xw+N    # int of the end of the range
  y1=yw-N
  y2=yw+N
  dx = []    # for the first deriv in x 
  dy = []
  if dw == 0:
    for i in range(2*N):
      dx.append( min(data[y-dw:y+dw, x1+i],maxflux)-min(data[y-dw:y+dw, x1+i+1],maxflux) )
      dy.append( min(data[y1+i, x-dw:x+dw],maxflux)-min(data[y1+i+1, x-dw:x+dw],maxflux) )
      if verbose:
        print x1+i, min(mean(data[y-dw:y+dw, x1+i]),maxflux), dx[i], y1+i, min(mean(data[y1+i, x-dw:x+dw]),maxflux), dy[i]
  else:
    for i in range(2*N):
      #dx.append( mean(data[y-dw:y+dw, x1+i])-mean(data[y-dw:y+dw, x1+i+1]) )
      #dy.append( mean(data[y1+i, x-dw:x+dw])-mean(data[y1+i+1, x-dw:x+dw]) )
      #dx.append( min(median(data[y-dw:y+dw, x1+i]),maxflux)-min(median(data[y-dw:y+dw, x1+i+1]),maxflux) )
      #dy.append( min(median(data[y1+i, x-dw:x+dw]),maxflux)-min(median(data[y1+i+1, x-dw:x+dw]),maxflux) )
      # fix for old pyfits, etc. on hiltner, agung -- 
      dx.append( min([median(data[yw-dw:yw+dw, x1+i]),maxflux])-min([median(data[yw-dw:yw+dw, x1+i+1]),maxflux]) )
      dy.append( min([median(data[y1+i, xw-dw:xw+dw]),maxflux])-min([median(data[y1+i+1, xw-dw:xw+dw]),maxflux]) )
      if verbose:
        #print x1+i, data[y, x1+i], dx[i], y1+i, data[y+i,x], dy[i]
        # print x1+i, min(median(data[y-dw:y+dw, x1+i]),maxflux), dx[i], y1+i, min(median(data[y1+i, x-dw:x+dw]),maxflux), dy[i]
        print x1+i, min([median(median(data[y-dw:y+dw, x1+i])),maxflux]), dx[i], y1+i, min([median(median(data[y1+i, x-dw:x+dw])),maxflux]), dy[i]
  x1=dx.index(min(dx))
  x2=dx.index(max(dx))
  if verbose:
    print x1, x2, x-N+x1, x-N+x2
  xc=x-N+x1+0.5*(x2-x1)+1.
  y1=dy.index(min(dy))
  y2=dy.index(max(dy))
  if verbose:
    print y1, y2, y-N+y1, y-N+y2
  yc=y-N+y1+0.5*(y2-y1)+1.
  # print x, y, xc, yc
  # fluxc = median(data[int(yc)-5:int(yc)+5, int(xc)-5:int(xc)+5]) 
  fluxc = median(median(data[int(yc)-5:int(yc)+5, int(xc)-5:int(xc)+5])) 
  # estimate the bias level for that box
  biasc = median(median(data[int(yc)-N:int(yc)+N, int(xc)-N:int(xc)+N])) 
  if verbose: 
    print "centerbox(): ", xc, yc, fluxc, biasc 
  fluxc = fluxc-biasc
  return xc, yc, fluxc 

# Solve for x, y, and rotation 

def centerslit(data, x, y, maxflux, N=50, dw=9,verbose=0):
  # N is the search tolerance for box edges in pixels (could get from MMS)
  # dw is the half width of the stripe used for the gradient search
  xw=int(x) 
  yw=int(y) 
  x1=xw-N
  x2=xw+N
  y1=yw-N
  y2=yw+N
  dx = []
  dy = []
  if dw == 0:
    for i in range(2*N):
      dx.append( min(data[y-dw:y+dw, x1+i],maxflux)-min(data[y-dw:y+dw, x1+i+1],maxflux) )
      if verbose:
        print x1+i, min(mean(data[y-dw:y+dw, x1+i]),maxflux), dx[i], y1+i, min(mean(data[y1+i, x-dw:x+dw]),maxflux)
  else:
    for i in range(2*N):
      dx.append( min([median(data[yw-dw:yw+dw, x1+i]),maxflux])-min([median(data[yw-dw:yw+dw, x1+i+1]),maxflux]) )
      if verbose:
        print x1+i, min([median(median(data[y-dw:y+dw, x1+i])),maxflux]), dx[i], y1+i, min([median(median(data[y1+i, x-dw:x+dw])),maxflux])
  x1=dx.index(min(dx))
  x2=dx.index(max(dx))
  if verbose:
    print x1, x2, x-N+x1, x-N+x2
  xc=x-N+x1+0.5*(x2-x1)+1.
  yc=y
  # print x, y, xc, yc
  # fluxc = median(data[int(yc)-5:int(yc)+5, int(xc)-5:int(xc)+5]) 
  fluxc = median(median(data[int(yc)-5:int(yc)+5, int(xc)-5:int(xc)+5])) 
  if verbose: 
    print "centerbox(): ", xc, yc, fluxc 
  return xc, yc, fluxc 

def calcgeotrans(sx, sy, bx, by, fittype, verbose): 
  # create input file
  geomapinput = "_oalign.in" 
  geomapoutput = "_oalign.out" 
  if os.path.isfile(geomapoutput): 
    os.remove(geomapoutput) 
  fmt = "%8.2f %8.2f %8.2f %8.2f \n" 
  F = open(geomapinput, "w") 
  for i in range(len(sx)): 
    # Adjust data so center of rotation is the origin
    # x(ref), y(ref), x(in), y(in) 
    gline = fmt % (sx[i]-xrotc, sy[i]-yrotc, bx[i]-xrotc, by[i]-yrotc)  
    F.write(gline) 
  F.close() 
  # run geomap
  iraf.geomap(geomapinput, "_mods", "INDEF","INDEF","INDEF","INDEF", results=geomapoutput, fitgeometry=fittype, verbose="no", interac="no")
  dx, dy, theta, nx, ny = readgeofile(geomapoutput, fittype) 
  # draw the new positions on the display
  drawstars(nx, ny, color="magenta") 
  fmt = "%d %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %6.2f %6.2f" 
  if verbose: 
    for i in range(len(nx)): 
      # nx, ny (new star positions) are x(fit) and y(fit) 
      print fmt % (i, bx[i], by[i], sx[i], sy[i], nx[i], ny[i], nx[i]-bx[i], ny[i]-by[i]) 
    print "# Geomap fit yields dx = %4.2f dy = %4.2f theta = %6.4f" % (dx, dy, theta) 
  # this is not necessary anymore [20 Mar 2011] 
  # nx,ny,res = applytrans(sx, sy, bx, by, dx, dy, theta) 
  res = []
  for i in range(len(nx)): 
    res.append(math.sqrt((nx[i] - bx[i])**2 + (ny[i] - by[i])**2))
  return dx, dy, theta, nx, ny, res

def readgeofile(file, fittype): 
  F = open(file, "r") 
  M=F.readlines()[::]
  F.close()
  nx = []
  ny = []
  theta = 0.  
  for i in range(len(M)): 
    if find(M[i], "X and Y shift:") >= 0: 
      dx,dy = float(M[i].split()[5]), float(M[i].split()[6])  
    if find(M[i], "X and Y axis rotation:") >= 0: 
      theta = float(M[i].split()[6])
    if find(M[i][0], "#") <0 and len(M[i])>10: 
      # read the geomap x(fit) and y(fit) 
      nx.append(float(M[i].split()[4])) 
      ny.append(float(M[i].split()[5])) 
  for i in range(len(nx)): 
    nx[i] = nx[i] + xrotc
    ny[i] = ny[i] + yrotc
  if theta>180.: 
    theta = 360.-theta  
  return dx, dy, theta, nx, ny 
  
def applytrans(sx, sy, bx, by, dx, dy, theta): 
  thetar = theta*math.pi/180. 
  outx = []
  outy = []
  ds = []
  for i in range(len(bx)): 
    #outx.append(math.cos(thetar)*(sx[i]-dx) - math.sin(thetar)*(sy[i]-dy) ) 
    #outy.append(math.sin(thetar)*(sx[i]-dx) + math.cos(thetar)*(sy[i]-dy) ) 
    outx.append(math.cos(thetar)*sx[i] - math.sin(thetar)*sy[i] + dx)
    outy.append(math.sin(thetar)*sx[i] + math.cos(thetar)*sy[i] + dy)
    ds.append(math.sqrt((outx[i] - bx[i])**2 + (outy[i] - by[i])**2))
  return outx, outy, ds

def drawstars(x,y,color="magenta"): 
  # Create a region file from x,y coords and draw stars on the display 
  F = open("_oalignxy.reg", "w") 
  for i in range(len(x)): 
    regline = "point(%6.2f,%6.2f) # point=x width=3 color=%s\n" % (x[i], y[i], color)
    F.write(regline)
  F.close()
  os.spawnlp(os.P_WAIT, "xpaset", "xpaset", "-p", "oalign", "regions","file","_oalignxy.reg")
  # clean up 
  os.remove("_oalignxy.reg") 

def rotator(deg): 
  # rotator module from Ray Gonzalez
  if math.fabs(deg) > 1.: 
    print "Error: Rotator move of %.1f degrees exceeds the 1 deg safety limit" % (deg) 
    print "Please use the hand paddle in the dome" 
    return False
  else: 
    angle = deg/0.7222
    HOST = '140.252.83.100' # The remote host
    PORT = 8003             # The same port as used by the server
    command = 'ANGLE=' + repr(angle) + '\r' # server address
    s = socket.socket(socket.AF_INET, socket.SOCK_STREAM) # setup socket
    s.connect((HOST, PORT)) # connect to Rotator comtrol server socket
    s.send(command)         # send the command.
    data = s.recv(24)     # retrive data if any
    s.close()                 # Close the connection
    print 'oalign.py: Sent/Received Rotator Angle', repr(data)
    return True

############################
#### Script starts here ####
############################

# Parse the input files and various flags

files, flags = parseinput()
if len(flags)>0: 
  for i in range(len(flags)): 
    if find(flags[i], "d")>=0: 
      debug=1
      print "Debug mode set" 
    if find(flags[i], "v")>=0: 
      verbose=1
      if debug: 
        print "Verbose mode set" 
    if find(flags[i], "s")>=0: 
      shiftonly=1
      if debug: 
        print "Shift-only mode set" 
    if find(flags[i], "V")>=0: 
      version=1
      if debug: 
        print "Print version mode set" 
    if find(flags[i], "l")>=0: 
      longslit=1
      if debug: 
        print "Long slit acquisition mode set" 
    if find(flags[i], "b")>=0: 
      boxdelete=0
      if debug: 
        print "Using previous alignment star boxes"
#
# old sense of -b flag
#      boxdelete=1
#      if debug: 
#        print "Box reselection forced" 
#

if version: 
  print "\n%s v%s [%s]\n" % (argv[0], versNum, versDate) 
  exit(0)
  
# Check that the appropriate number of input files where specificed 

if len(files) < 2 or len(files) > 3:
  usage()
  exit(1) 

if longslit==0: 
  omsfile = files[0]
  maskimage = files[1]
  if len(files)==3: 
    starimage = files[2]
else: 
  maskimage = files[0]	# slit image in longslit mode 
  starimage = files[1]	# acq image in longslit mode 

start_ds9()
iraf.unlearn("imexamine") 

object = get_keyword(maskimage, "OBJECT")

# ?? hide behind verbose ??

if verbose:
  if longslit==0: 
    print "Maskimage %s title is: %s" % (maskimage, object) 
  else: 
    print "Slitimage %s title is: %s" % (maskimage, object) 

# The FITS header card SECPIX contains the unbinned pixel scale in "/pixel

# Note: this assumes the mask image and field image have the same binning
# pixscale = float( get_keyword(maskimage, "SECPIX") )
# turned off because R4K does not have SECPIX keyword [29 Oct 2011 | PM]
pixscale = 0.273

# Parse the oms file to identify the preliminary box coordinates
if longslit==0: 
  x,y,dx,dy = parseoms(omsfile) 
# Establish the name for the box files 
starti = 0
if maskimage.rfind("/")>=0: 
  starti = maskimage.rfind("/")+1 
boxfile = strip(maskimage[starti::], ".fits")+".box"
boxregfile = strip(maskimage[starti::], ".fits")+".reg"
if verbose: 
  if longslit==0: 
    print "Alignment box coordinates file: ", boxfile
  else: 
    print "Slit coordinates file: ", boxfile
# Delete the box file if boxdelete=True and it exists
if boxdelete: 
  if os.path.isfile(boxfile): 
    if verbose: 
      # "Flag -b is set, so deleting %s" % (boxfile)
      "Deleting old box file %s" % (boxfile) 
      os.remove(boxfile)
  if os.path.isfile(boxregfile): 
      os.remove(boxregfile)

### Measure the box positions or get them from an existing file -- 

# If boxfile exists, do not repeat marking of box

if longslit==0: 
  if os.path.isfile(boxfile) and os.path.isfile(boxregfile):
    if verbose:
      print "\n** Box coordinate files %s and %s found." % (boxfile, boxregfile) 
    print "\n** Proceeding directly to marking alignment stars" 
    boxxc, boxyc, boxid, boxfluxc = readboxfile(boxfile)
  # otherwise generate the boxfile interactively 
  else:
    if verbose:
      print "Box coordinate files %s and/or %s do not exist" % (boxfile, boxregfile)
    # make sure both are gone 
    if os.path.isfile(boxfile): 
      os.remove(boxfile)
    if os.path.isfile(boxregfile): 
      os.remove(boxregfile)
    boxxc, boxyc, boxid, boxfluxc = getboxcoords(pixscale, maskimage) 
else: 
  if os.path.isfile(boxfile) and os.path.isfile(boxregfile):
    if verbose:
      print "\n** Box coordinate files %s and %s found." % (boxfile, boxregfile) 
    print "\n** Proceeding directly to marking alignment stars" 
    boxxc, boxyc, boxid, boxfluxc = readboxfile(boxfile)
  # otherwise generate the boxfile interactively 
  else:
    if verbose:
      print "Box coordinate files %s and/or %s do not exist" % (boxfile, boxregfile)
    if os.path.isfile(boxfile): 
      os.remove(boxfile)
    if os.path.isfile(boxregfile): 
      os.remove(boxregfile)
    boxxc, boxyc, boxid, boxfluxc = getslitcoords(pixscale, maskimage) 


### Display the box coordinates in both frames 

os.spawnlp(os.P_WAIT, "xpaset", "xpaset", "-p", "oalign", "regions","deleteall")
os.spawnlp(os.P_WAIT, "xpaset", "xpaset", "-p", "oalign", "regions","file",boxregfile) 
os.spawnlp(os.P_WAIT, "xpaset", "xpaset", "-p", "oalign", "frame","2") 

#if len(files) == 2 and mask and longslit==0: 
#  boxz1, boxz2 = getdisplaypars(boxfluxc) 
#  print "\n** Doing alignment with the mask image." 
#  iraf.display(maskimage, 2, zscale="no", zrange="no", z1=boxz1, z2=boxz2)
#else:
#  iraf.display(starimage, 2)

iraf.display(starimage, 2)

if os.path.isfile("_oalign.log"):
  os.remove("_oalign.log") 

### Measure the star positions

if len(files) == 2: 
  starx, stary = getstarcoords(maskimage) 
else: 
  starx, stary = getstarcoords(starimage) 

# Check that we have the right number of each

if longslit == 0: 
  if len(starx) != len(boxxc): 
    print "\n** ERROR: number of alignment stars and boxes are not equal..."
    print "   %s aborting.\n" % (argv[0])
    end_ds9()
    # ?? any other cleanup we should do on abort?
    exit()
  if len(starx)<2: 
    print "\n** ERROR: At least two alignment stars required to solve"
    print "   for mask translation and rotation - %s aborting.\n" % (argv[0])
    end_ds9()
    # ?? any other cleanup we should do on abort?
    exit()

nobj = len(starx) 

####################################
### Calculate the transformation ###
####################################

# Calculate translation-only or translation+rotation 

if longslit==0: 
  if shiftonly: 
    # calculate shift only 
    if verbose:
      print "Geomap fit" 
    dx,dy,theta,nstarx,nstary,res = calcgeotrans(starx, stary, boxxc, boxyc, "shift", verbose)
    theta = 0. 
  else: 	
    # calculate shift and rotation
    if verbose: 
      print "Geomap fit"  
    dx,dy,theta,nstarx,nstary,res = calcgeotrans(starx, stary, boxxc, boxyc, "rotate", verbose)
else: 
  if verbose: 
    print starx, boxxc, stary, boxyc
  if type(boxxc) is float: 
    boxxc = [boxxc] 
    boxyc = [boxyc] 
  dx=boxxc[0]-starx[0]
  dy=boxyc[0]-stary[0]
  nstarx = starx
  nstary = stary
  theta=0. 
  res=[0.]

# Issue a warning about large residuals 

mn = mean(res)
std = std(res) 
fmt = "%d %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %6.2f %6.2f %6.2f" 
if mn>20.: 
  print "\n** WARNING: Mean offset relative to solution is %4.2f pixels" % (mn) 
  print "ID BOX_X BOX_Y Star_X Star_Y New_Star_X New_Star_Y Delta_X Delta_Y Delta_S"
  for i in range(len(res)): 
    print fmt % (i, boxxc[i], boxyc[i], starx[i], stary[i], nstarx[i], nstary[i], nstarx[i] - boxxc[i], nstary[i]-boxyc[i], res[i])
if verbose: 
  print "Mean offset relative to solution is %4.2f pixels" % (mn) 
  print "ID BOX_X BOX_Y Star_X Star_Y New_Star_X New_Star_Y Delta_X Delta_Y Delta_S"
  for i in range(len(res)): 
    print fmt % (i, boxxc[i], boxyc[i], starx[i], stary[i], nstarx[i], nstary[i], nstarx[i] - boxxc[i], nstary[i]-boxyc[i], res[i])

# Print the output (units are arcseconds and degrees)

theta = -1*theta 	# sign flip
dxa = (dx*pixscale)
dya = (dy*pixscale)
gscale = 3.4821
gdx = -1.*(dy*gscale)  # note x and y reversed 	[-1 added 29 Oct 2011 by PM] 
gdy = (dx*gscale) 

if longslit==0: 
  print "\nComputed Mask Alignment Offset:"
else: 
  print "\nComputed Slit Alignment Offset:"
print "   dx: %7.3f arcsec\n   dy: %7.3f arcsec\n  dPA: %7.4f degrees" % (dxa, dya, theta)
print "\nThese are the guider offsets to apply:"
print "   gdx: %7.4f steps\n   gdy: %7.4f steps\n" % (gdx, gdy)
if longslit==0: 
  print "\nAnd enter this rotation offset in the Rotator GUI:"
  print "   dPA: %7.4f degrees\n" % (theta)  
else: 
  print "** Note that the gdy offset is the critical one\n" 

# Clean up (uncomment more below when in "production")

if os.path.isfile("_oalign.log"):
  os.remove("_oalign.log")
if os.path.isfile("_oalign.reg"):
  os.remove("_oalign.reg")
# if os.path.isfile("_oalign.out"):
#   os.remove("_oalign.out")
# if os.path.isfile("_oalign.in"):
#   os.remove("_oalign.in")
  
# end_ds9()

print "\n%s done." % (argv[0])