import sys
import os
import casac
import string
import numpy
import inspect

from parameter_check import *

# Get visibilities
# Sort visibilities by baseline length.
# For each bin of nvisperbin visibilities,
#    Calculate the phase dispersion after subtracting the expected phase from a
#    source in direction.
# Save those dispersions.
#
def coherence(vis=None, deltal=None, deltam=None, nvisperbin=50):
	"""Analyze phase noise of an ALMA dataset:

	Keyword arguments:
	vis  -- Name of measurement set.  Example: vis='mysim.ms'
	deltal -- East-west component of (peak position - phase tracking
	          center), in arcsec.
	deltam -- North-south component of (peak position - phase tracking
	          center), in arcsec.
	nvisperbin -- How many visibilities to include in each bin for
	              determining phase dispersions.

	"""
        a=inspect.stack()
        stacklevel=0
        for k in range(len(a)):
          if (string.find(a[k][1], 'ipython console') > 0):
                stacklevel=k
	myf=sys._getframe(stacklevel).f_globals
	myf['__last_task']='coherence'
	myf['taskname']='coherence'
#	###fill unfilled parameters with defaults
#	myf['update_params'](func=myf['taskname'], printtext=False)
        tb=myf['tb']
#	default=myf['default']
			
#	vis=myf['vis']
#        deltal=myf['deltal']
#        deltam=myf['deltam']
#        nvisperbin=myf['nvisperbin']
	###
	#
        #Add type/menu/range error checking here
#	if type(mask)==str: mask=[mask]
        arg_names  = ['vis', 'deltal', 'deltam', 'nvisperbin']
        arg_values = [vis,    deltal,   deltam,   nvisperbin]
        arg_types  = [str,    float,    float,    int]
	#parameter_printvalues(arg_names,arg_values,arg_types)
        try:
                parameter_checktype(arg_names, arg_values, arg_types)
#		parameter_checkmenu('mode',mode,['mfs','channel','velocity'])
        except TypeError, e:
                print "coherence -- TypeError: ", e
		return
        except ValueError, e:
                print "coherence -- OptionError: ", e
                return
        ###

	print "Made it past the argument checking."

	try:
#	    pdb.set_trace()
# Get visibilities
            msfile = vis
            tb.open(msfile, nomodify=True)
	    
	    ### shape of npol, nchan, nrows
	    dat = numpy.array(tb.getcol('DATA'))

	    uvw = numpy.array(tb.getcol('UVW'))  ### uvw is in meters in the ms
	    nvis = len(uvw[0,:])

	    print "nvis =", nvis

	    # Remove phase due to distance between the source and the phase
	    # tracking center.  This will probably only work with 1
	    # polarization and 1 channel.
	    RArad  = numpy.pi * deltal / (180.0 * 3600.0)
	    decrad = numpy.pi * deltam / (180.0 * 3600.0)
#	    dat *= numpy.exp(complex(0.0, 2.0 * numpy.pi * 
#				     (RArad * uvw[0,:] + decrad * uvw[1,:])))
#
#	    print "dat has been deshifted."
	    
	    phase = numpy.arctan2(dat.imag, dat.real)
	    
	    # Sort visibilities by baseline length.
	    uvdist = 0.001 * numpy.sqrt(uvw[0,:]**2+uvw[1,:]**2)
	    sortedindices = numpy.argsort(uvdist)

	    # For each bin of nvisperbin visibilities,
	    nbins = int(nvis / nvisperbin)
	    phasedisp = numpy.zeros(nbins)  # Just to declare them as NumArrays
	    bll       = numpy.zeros(nbins)  # (which default to float).
	    for i in xrange(nbins):
		    binindices = sortedindices[nvisperbin * i : nvisperbin * (i
	    + 1) + 1]
		    phasedisp[i] = numpy.std(phase.take(binindices))
		    bll[i]       = numpy.mean(uvdist.take(binindices))
		    
	    # Print those dispersions.
	    print "baseline (km)\tphase dispersion (radians)"
	    for i in xrange(nbins):
		    print "%.2f\t%.2f" % (bll[i], phasedisp[i])

	except TypeError, e:
                print "coherence -- TypeError: ", e
		return
        except ValueError, e:
                print "coherence -- OptionError: ", e
                return
	except Exception, instance:
		print '***Error***',instance
		return

#        saveinputs=myf['saveinputs']
#        saveinputs('coherence','coherence.last')

###End of coherence