#! /usr/bin/env python
# read a liris image corresponding to any of the longslit
#
# command line parameters
# 1: name of image
# 2: file containing reference positions
# 3: verbose (give detailed information about star and region statistics)
#
# Last version:
#    Jul-08 - Corrected algorithm to take stars as count excess over 
#             quartiles (J.Acosta) 

import pyfits
import sys
import scipy
#import numdisplay
#from pylab import *


def printHelp():

   print " ----------------------------------------------------------"
   print " Check mask centering in imaging mode when mask is inserted "
   print " Syntax: recenter_mos imagen.fit mask_center.mask [verbose]"
   print "              "
   print " ----------------------------------------------------------"
   return

def extract_regions(im,xref,yref,npos,rad,npix_ref,npix_star,verb):
   
   xscent = []
   yscent = []
   wcent = []
   xcent = []
   ycent = []

  ## image to obtain hole center is larger to have background and use 
  ## image with not subtracted sky
  ## image to obtain star center is smaller to contain only sky and 
  ## sky must be subtracted
   for ireg in range(0,npos):
     ## subtract -1 to take into account zero start in Py
     if verb == 1: print 'Reference position ',ireg+1
     xr1 = int(xref[ireg] - npix_ref+0.5) - 1 
     xr2 = int(xref[ireg] + npix_ref+0.5) - 1
     yr1 = int(yref[ireg] - npix_ref+0.5) - 1 
     yr2 = int(yref[ireg] + npix_ref+0.5) - 1
     imref = im[yr1:yr2,xr1:xr2]
     xs1 = int(xref[ireg] - npix_star+0.5) - 1
     xs2 = int(xref[ireg] + npix_star+0.5) 
     ys1 = int(yref[ireg] - npix_star+0.5) - 1
     ys2 = int(yref[ireg] + npix_star+0.5) 
     if verb == 1: print '  region limits for star center ', xs1,xs2,ys1,ys2
     if verb == 1: print '  region limits for hole center ', xr1,xr2,yr1,yr2
     imstar = im[ys1:ys2,xs1:xs2]
     sky = scipy.median(scipy.median(imstar))
     ind3qua = scipy.where(imstar > sky)
     sky3qua = scipy.median(imstar[ind3qua]) - sky
     ind1qua = scipy.where(imstar < sky)
     sky1qua = sky - scipy.median(imstar[ind1qua])
     skysg = (sky3qua+sky1qua)/2.*1.2 
     if verb == 1: print '  Sky value and sigma', sky,skysg
     ## get center of mass for reference hole
     imref_sky = imref < float((sky*1.1))
     imref_clip = imref >= float((sky*1.1))
     imref = imref_clip * float(sky*1.1) + imref_sky * imref 
     nx = xr2 - xr1 
     ny = yr2 - yr1 
     #if (ireg == 9999): 
     #  numdisplay.display(imref,z1=0,z2=3500)
     shy = scipy.add.reduce(imref)
     shx = scipy.add.reduce(imref,axis=1)
     #print shy
     #print shx

     # xcenter 
     xind = scipy.arange(nx) + xr1 + 1
     shyx = shy * xind
     #print shy
     xcent.append(float(scipy.add.reduce(shyx)) / float(scipy.add.reduce(shy)))

     # ycenter
     yind = scipy.arange(ny) + yr1 + 1
     shxy = shx * yind
     ycent.append(float(scipy.add.reduce(shxy) / scipy.add.reduce(shx)))
     if verb == 1: print '  Expected Center of the hole [x,y]:', xref[ireg],yref[ireg]
     if verb == 1: print '  Measured Center of the hole [x,y]:', xcent[ireg],ycent[ireg]

     ## get center of mass for star
     nx = xs2 - xs1
     ny = ys2 - ys1
     # cut values below sky 
     imstar_clip = imstar < float(sky*1.)
     imstar_val = imstar >= float(sky*1.)
     imstar =  imstar_clip* float(sky*1.) + imstar_val * imstar
     imstar = imstar  - sky
     imstar_wcent = imstar >= float(skysg*3.)
     wcent.append(float(scipy.add.reduce(scipy.add.reduce(imstar_wcent))))
     if verb == 1: print '  Number of pixels above 3*sigma of sky ',wcent[ireg],' over a total of ',nx*ny
     #if (ireg == 0): 
     #  numdisplay.display(imstar,z1=-500,z2=5200)
     shy = scipy.add.reduce(imstar)
     shx = scipy.add.reduce(imstar,axis=1)

     # xcenter 
     xind = scipy.arange(nx) + xs1 + 1
     shyx = shy * xind
     xscent.append(float(scipy.add.reduce(shyx)) / float(scipy.add.reduce(shy)))

     # ycenter
     yind = scipy.arange(ny) + ys1 + 1
     shxy = shx * yind
     yscent.append(float(scipy.add.reduce(shxy)) / float(scipy.add.reduce(shx)))
     if verb == 1: print '  Measured star centroid ', xscent[ireg],yscent[ireg]
     if verb == 1: print ' '
     if verb == 1: print ' '

   wtot = float(scipy.add.reduce(wcent)) 
   if wtot > 0:
      xoff = []
      yoff = []
      for ireg in range(0,npos):
       print 'Reference Position ',ireg+1,' weight ',wcent[ireg]/wtot
       print '  Offset - x,y [pixels]',-(xscent[ireg] - xcent[ireg]),yscent[ireg] - ycent[ireg]
       #print 'Off-y [pixels]',yscent[ireg] - ycent[ireg]
       xoff.append(float((xscent[ireg] - xcent[ireg])*wcent[ireg]))
       yoff.append(float((yscent[ireg] - ycent[ireg])*wcent[ireg]))

      print '  '
      print 'Average Star-Reference Offset - x,y [arcsec]:'
      print '              ', -float(scipy.add.reduce(xoff))/wtot*0.25, float(scipy.add.reduce(yoff))/wtot*0.25   
   else: 
      printError_TooFaint()

   
def printError_TooFaint():

   print " ERROR: Offset cannot be determined.  "
   print "        Centering targets are too faint "
   print "              "
   return
  
# Start main 
if __name__ == '__main__':
 _args = sys.argv[1:]
 _nargs = len(_args)

if _nargs < 2:
  printHelp()
else:     
  filein = _args[0]
  im = pyfits.getdata(filein) 
  filecent = _args[1]
  fcent = open(filecent,'r')
  npos = 0
  iref = []
  xref = []
  yref = []
  for line in fcent:
    data = line.split()
    if data[0]=='REF': 
      iref.append(data[0])
      xref.append(float(data[1]))
      yref.append(float(data[2]))
      npos = npos+1

  fcent.close()  
  
  # check verbose
  if _nargs == 3: 
    verb  = 1 
  else:
    verb = 0  

  ## extract regions centered at refpos
  rad = 3./0.25
  npix_ref = int(rad)
  npix_star = 7

  extract_regions(im,xref,yref,npos,rad,npix_ref,npix_star,verb)