def plot_track_props(tracks, nx, ny, len_cutoff=20):
    pl.ioff()
    wdist = wraparound_dist(nx, ny)
    val_fig = pl.figure()
    area_fig = pl.figure()
    psn_fig = pl.figure()
    delta_vals = []
    delta_dists = []
    for tr in tracks:
        if len(tr) < len_cutoff:
            continue
        idxs, regs = zip(*tr)
        delta_vals.extend([abs(regs[idx].val - regs[idx + 1].val) for idx in range(len(regs) - 1)])
        dists = [wdist(regs[i].loc, regs[i + 1].loc) for i in range(len(regs) - 1)]
        delta_dists.extend([abs(dists[idx] - dists[idx + 1]) for idx in range(len(dists) - 1)])
        pl.figure(val_fig.number)
        pl.plot(idxs, [reg.val for reg in regs], "s-", hold=True)
        pl.figure(area_fig.number)
        pl.semilogy(idxs, [reg.area for reg in regs], "s-", hold=True)
        pl.figure(psn_fig.number)
        pl.plot(idxs[:-1], dists, "s-", hold=True)
    pl.figure(val_fig.number)
    pl.savefig("val_v_time.pdf")
    pl.figure(area_fig.number)
    pl.savefig("area_v_time.pdf")
    pl.figure(psn_fig.number)
    pl.savefig("psn_v_time.pdf")
    pl.figure()
    pl.hist(delta_vals, bins=pl.sqrt(len(delta_vals)))
    pl.savefig("delta_vals.pdf")
    pl.figure()
    pl.hist(delta_dists, bins=pl.sqrt(len(delta_dists)))
    pl.savefig("delta_dists.pdf")
    pl.close("all")

def simple_hierarchical_data(n):
    """ Generate data based on the simple one-way hierarchical model
    given in section 3.1.1::

        y[i,j] | alpha[j], sigma^2 ~ N(alpha[j], sigma^2) i = 1, ..., n_j, j = 1, ..., J;
        alpha[j] | mu, tau^2 ~ N(mu, tau^2) j = 1, ..., J.

        sigma^2 ~ Inv-Chi^2(5, 20)
        mu ~ N(5, 5^2)
        tau^2 ~ Inv-Chi^2(2, 10)

    Parameters
    ----------
    n : list, len(n) = J, n[j] = num observations in group j
    """

    inv_sigma_sq = mc.rgamma(alpha=2.5, beta=50.0)
    mu = mc.rnormal(mu=5.0, tau=5.0 ** -2.0)
    inv_tau_sq = mc.rgamma(alpha=1.0, beta=10.0)

    J = len(n)
    alpha = mc.rnormal(mu=mu, tau=inv_tau_sq, size=J)
    y = [mc.rnormal(mu=alpha[j], tau=inv_sigma_sq, size=n[j]) for j in range(J)]

    mu_by_tau = mu * pl.sqrt(inv_tau_sq)
    alpha_by_sigma = alpha * pl.sqrt(inv_sigma_sq)
    alpha_bar = alpha.sum()
    alpha_bar_by_sigma = alpha_bar * pl.sqrt(inv_sigma_sq)

    return vars()

def renormalize(x_unpurt,x_before,x_purt,epsilon,N):
    # BEFORE ANYTHING: make sure particles near boundaries are shuffeled into places where where the
    # seam is not between any purturbed and fudicial trajectories.
    x_unpurt,x_purt = shuff(x_unpurt,x_purt,N)

    # The trajectory we are going to be returning is going to be the new one for the next run. lets
    # call it
    x_new = pl.copy(x_unpurt)
    # copied it because we are going to add the small amounts to it to purturb it.

    # lets find a vector pointing in the direction of the trajectories path. For this we need the
    # fiducual point at t-dt, which is given to us in the function as x_before. find the vector
    # between x_before and x_unpurt
    traj_vec = x_unpurt-x_before 
    # normalize it
    traj_vec = traj_vec/pl.sqrt(pl.dot(traj_vec,traj_vec))
    print('traj_vec magnitude (should be 1): ' + str(pl.sqrt(pl.dot(traj_vec,traj_vec))))

    # Now lets see how close the vector pointing from the fidicial to the perturbed trajectorie is
    # to orthogonal with the trajectory... should get closer to 1 as we check more because it should
    # be aligning itself with the axis of greatest expansion and that should be orthogonal.
    # First normalize the difference vector
    diff_vec = x_unpurt - x_purt
    # normalize it
    diff_vec = diff_vec/pl.sqrt(pl.dot(diff_vec,diff_vec))
    print('diff_vec magnitude (should be 1): ' + str(pl.sqrt(pl.dot(diff_vec,diff_vec))))
    print('normalized(x_unpurt-x_purt)dot(traj_vec)  (should get close to 0): '+ str(pl.dot(diff_vec,traj_vec)))

    # for now lets just return a point moved back along the difference vector. no gram shmidt or
    # anything.
    return x_new + epsilon*diff_vec

def haversine (latlong1, latlong2, r):

    deltaLatlong = latlong1 - latlong2
    
    dLat = deltaLatlong[0]
    dLon = deltaLatlong[1]

    lat1 = latlong1[0]
    lat2 = latlong2[0]

    a = (sin (dLat/2) * sin (dLat/2) +
         sin (dLon/2) * sin (dLon/2) * cos (lat1) * cos (lat2))
    c = 2 * arctan2 (sqrt (a), sqrt (1-a))
    d = r * c

    # initial bearing
    y = sin (dLon) * cos (lat2)
    x = (cos (lat1)*sin (lat2) -
         sin (lat1)*cos (lat2)*cos (dLon))
    b1 = arctan2 (y, x);

    # final bearing
    dLon = -dLon
    dLat = -dLat
    tmp = lat1
    lat1 = lat2
    lat2 = tmp
    y = sin (dLon) * cos (lat2)
    x = (cos (lat1) * sin (lat2) - 
         sin (lat1) * cos (lat2) * cos (dLon))
    b2 = arctan2 (y, x)
    b2 = mod ((b2 + pi), 2*pi)

    return (d, b1, b2)

    def sample(self, model, evidence):
        g = evidence['g']
        h = evidence['h']
        C = evidence['C']
        z = evidence['z']
        shot_id = evidence['shot_id']
        noise_proportion = evidence['noise_proportion']
        observation_var_g = evidence['observation_var_g']
        observation_var_h = evidence['observation_var_h']

        canopy_cover = model.known_params['canopy_cover']
        z_min = model.known_params['z_min']
        z_max = model.known_params['z_max']

        prior_p = model.hyper_params['T']['p']

        N = len(z)
        T = zeros(N)
        noise_rv = stats.uniform(z_min, z_max - z_min)
        min_index = min(z.index)
        for i in shot_id.index:
            l = zeros(3)
            index = i-min_index
            shot_index = shot_id[i]-min(shot_id)
            l[0] = noise_proportion*noise_rv.pdf(z[i])
            g_norm = stats.norm(g[shot_index], sqrt(observation_var_g))
            C_i = canopy_cover[C[shot_index]]
            l[1] = (1-noise_proportion)*(1-C_i)*g_norm.pdf(z[i])
            h_norm = stats.norm(h[shot_index] + g[shot_index], sqrt(observation_var_h))
            if z[i] > g[shot_index]+3:
                l[2] = (1-noise_proportion)*(C_i)*h_norm.pdf(z[i])
            p = l/sum(l)
            T[index] = Categorical(p).rvs()

        return T

    def rank_by_distance_bhatt(self, qkeys, ikeys, rkeys, dists):
        """
        ::

            Reduce timbre-channel distances to ranks list by ground-truth key indices
            Bhattacharyya distance on timbre-channel probabilities and Kullback distances
        """
        # timbre-channel search using pre-computed distances
        ranks_list = []
        t_keys, t_lens = self.get_adb_lists(0) 
        rdists=pylab.ones(len(t_keys))*float('inf')
        qk = self._get_probs_tc(qkeys)
        for i in range(len(ikeys[0])): # number of include keys
            ikey=[]
            dk = pylab.zeros(self.timbre_channels)
            for t_chan in range(self.timbre_channels): # timbre channels
                ikey.append(ikeys[t_chan][i])
                try: 
                    # find dist of key i for query
                    i_idx = rkeys[t_chan].index( ikey[t_chan] ) # dataset include-key match
                    # the reduced distance function in include_keys order
                    # distance is Bhattacharyya distance on probs and dists
                    dk[t_chan] = dists[t_chan][i_idx]
                except:
                    print "Key not found in result list: ", ikey, "for query:", qkeys[t_chan]
                    raise error.BregmanError()
            rk = self._get_probs_tc(ikey)
            a_idx = t_keys.index( ikey[0] ) # audiodb include-key index
            rdists[a_idx] = distance.bhatt(pylab.sqrt(pylab.absolute(dk)), pylab.sqrt(pylab.absolute(qk*rk)))
        #search for the index of the relevant keys
        rdists = pylab.absolute(rdists)
        sort_idx = pylab.argsort(rdists)   # Sort fields into database order
        for r in self.ground_truth: # relevant keys
            ranks_list.append(pylab.where(sort_idx==r)[0][0]) # Rank of the relevant key
        return ranks_list, rdists

    def intersection(self,ray) :
        cv = self.placement.location+self.placement.orientation*self.radcurv
        dv = ray.p0 - self.placement.orientation*self.radcurv - self.placement.location        
        a = 1
        b = 2*pl.linalg.dot(ray.d,dv)
        c = pl.linalg.dot(dv,dv)-self.radcurv**2
        
        qs  = b**2-4*a*c
        if qs == 0 :
            lam = -b/(2*a)
        elif qs < 0 :
            lam = None
        else :
            lamp = (-b+pl.sqrt(b**2-4*a*c))/(2*a)
            lamn = (-b-pl.sqrt(b**2-4*a*c))/(2*a)
            pd   = pl.linalg.norm(ray.propagate(lamp)-ray.p0)
            nd   = pl.linalg.norm(ray.propagate(lamn)-ray.p0)
#            lam = min(lamp,lamn)
            
            if self.radcurv > 0 :
                lam = min(lamp,lamn)
            elif self.radcurv < 0 :
                lam = max(lamp,lamn)
            
            # assign intersection
        ray.p1 = ray.propagate(lam)

def displayData(X):
    print "Visualizing"
    m, n = X.shape
    width = round(sqrt(n))
    height = width
    display_rows = int(floor(sqrt(m)))
    display_cols = int(ceil(m/display_rows))

    print "Cell width:", width
    print "Cell height:", height    
    print "Display rows:", display_rows
    print "Display columns:", display_cols
        
    display = zeros((display_rows*height,display_cols*width))

    # Iterate through the training sets, reshape each one and populate
    # the display matrix with the letter matrixes.    
    for xrow in range(0, m):
        rowindex = divide(xrow, display_cols)
        columnindex = remainder(xrow, display_cols)
        rowstart = int(rowindex*height)
        rowend = int((rowindex+1)*height)
        colstart = int(columnindex*width)
        colend = int((columnindex+1)*width)
        display[rowstart:rowend, colstart:colend] = X[xrow,:].reshape(height,width).transpose()
         
    imshow(display, cmap=get_cmap('binary'), interpolation='none')
    
    # Show plot without blocking
    draw()    

 def __calculate__(self):
     global USE_IDENTITY_LINE
     sd1 = (self.signal_plus - self.signal_minus) / pl.sqrt(2)
     if USE_IDENTITY_LINE:
         return pl.sqrt(pl.sum((sd1 ** 2)) / len(self.signal_plus))
     else:
         return pl.sqrt(pl.var(sd1))

    def calculateFDunc(self):
        #Calculates the uncertainty of the FFT according to:
        #   - J. M. Fornies-Marquina, J. Letosa, M. Garcia-Garcia, J. M. Artacho, "Error Propagation for the transformation of time domain into frequency domain", IEEE Trans. Magn, Vol. 33, No. 2, March 1997, pp. 1456-1459
        #return asarray _tdData
        #Assumes tha the amplitude of each time sample is statistically independent from the amplitude of the other time
        #samples

        # Calculates uncertainty of the real and imaginary part of the FFT and ther covariance
        unc_E_real = []
        unc_E_imag = []
        cov = []
        for f in self.getfreqs():
            unc_E_real.append(py.sum((py.cos(2*py.pi*f*self._tdData.getTimes())*self._tdData.getUncEX())**2))
            unc_E_imag.append(py.sum((py.sin(2*py.pi*f*self._tdData.getTimes())*self._tdData.getUncEX())**2))
            cov.append(-0.5*sum(py.sin(4*py.pi*f*self._tdData.getTimes())*self._tdData.getUncEX()**2))
        
        unc_E_real = py.sqrt(py.asarray(unc_E_real))
        unc_E_imag = py.sqrt(py.asarray(unc_E_imag))
        cov = py.asarray(cov)
        
        # Calculates the uncertainty of the modulus and phase of the FFT
        unc_E_abs = py.sqrt((self.getFReal()**2*unc_E_real**2+self.getFImag()**2*unc_E_imag**2+2*self.getFReal()*self.getFImag()*cov)/self.getFAbs()**2)
        unc_E_ph = py.sqrt((self.getFImag()**2*unc_E_real**2+self.getFReal()**2*unc_E_imag**2-2*self.getFReal()*self.getFImag()*cov)/self.getFAbs()**4)
        
        t=py.column_stack((self.getfreqs(),unc_E_real,unc_E_imag,unc_E_abs,unc_E_ph))
        return self.getcroppedData(t)  

def render_network(A):
    [L, M] = shape(A)
    sz = int(sqrt(L))
    buf = 1
    A = asarray(A)

    if floor(sqrt(M)) ** 2 != M:
        m = int(sqrt(M / 2))
        n = M / m
    else:
        m = int(sqrt(M))
        n = m

    array = -ones([buf + m * (sz + buf), buf + n * (sz + buf)], "d")

    k = 0
    for i in range(m):
        for j in range(n):
            clim = max(abs(A[:, k]))
            x_offset = buf + i * (sz + buf)
            y_offset = buf + j * (sz + buf)
            array[x_offset : x_offset + sz, y_offset : y_offset + sz] = reshape(A[:, k], [sz, sz]) / clim

            k += 1
    return array

def haversine(location1, location2=None):  # calculates great circle distance
    __doc__ = """Returns the great circle distance of the given
    coordinates.
    
    INPUT:  location1 = ((lat1, lon1), ..., n(lat1, lon1))
           *location2 = ((lat2, lon2), ..., n(lat2, lon2))
           *if location2 is not given a square matrix of distances
             for location1 will be put out
    OUTPUT: distance in km
            (dist1  ...  ndist
              :            : 
             ndist1 ...  ndist)
            shape will depend on the input
    METHOD: a = sin(dLat / 2) * sin(dLat / 2) + 
                sin(dLon / 2) * sin(dLon / 2) * 
                cos(lat1) * cos(lat2)
            c = 2 * arctan2(sqrt(a), sqrt(1 - a))
            d = R * c
            
            where R is the earth's radius (6371 km)
            and d is the distance in km"""
    
    from itertools import product, combinations
    from pylab import   deg2rad, sin, cos, arctan2, \
                        meshgrid, sqrt, array, arange
    
    if location2: 
        location1 = array(location1, ndmin=2)
        location2 = array(location2, ndmin=2)
    elif location2 is None:
        location1 = array(location1, ndmin=2)
        location2 = location1.copy()
    
    # get all combinations using indicies
    ind1 = arange(location1.shape[0])
    ind2 = arange(location2.shape[0])
    ind  = array(list(product(ind1, ind2)))
    
    # using combination inds to get lats and lons
    lat1, lon1 = location1[ind[:,0]].T
    lat2, lon2 = location2[ind[:,1]].T
    
    # setting up variables for haversine
    R = 6371.
    dLat = deg2rad(lat2 - lat1)
    dLon = deg2rad(lon2 - lon1)
    lat1 = deg2rad(lat1)
    lat2 = deg2rad(lat2)
    
    # haversine formula
    a = sin(dLat / 2) * sin(dLat / 2) + \
        sin(dLon / 2) * sin(dLon / 2) * \
        cos(lat1) * cos(lat2)
    c = 2 * arctan2(sqrt(a), sqrt(1 - a))
    d = R * c
    
    # reshape accodring to the input
    D = d.reshape(location1.shape[0], location2.shape[0])
    
    return D

def ICeuklid_to_ICcircle(IC):
    """
    converts from IC in euklidean space to IC in circle parameters (rotational invariant).
    The formats are:
    IC_euklid: [x, y, z, vx, vy, vz]
    IC_circle: [y, vy, |v|, |l|, phiv], where |v| is the magnitude of CoM velocity, |l| 
        is the distance from leg1 (assumed to be at [0,0,0]) to CoM, and phiv the angle
        of the velocity in horizontal plane wrt x-axis
    *NOTE* for re-conversion, the leg position is additionally required
    
    :args:
        IC (6x float): the initial conditions in euklidean space

    :returns:
        IC (5x float): the initial conditions in circular coordinates
    
    """
    x,y,z,vx,vy,vz = IC
    v = sqrt(vx**2 + vy**2 + vz**2)
    l = sqrt(x**2 + y**2 + z**2)
    #phiv = arctan2(vz, vx)
    #phiv = arctan2(-vz, vx)
    phiv = -arctan2(-vz, vx)
    #phix = arctan2(-z, -x)
    phix = arctan2(z, -x)
    # warnings.warn('TODO: fix phi_x (add)')
    # print "phix:", phix * 180 / pi
    return [y, vy, v, l, phiv + phix]

def f(x, t):
    # for now masses just = 1.0

    # the 4.0 only works for 2D
    N = len(x) / 4
    xdot = pl.array([])

    for i in range(N):
        temp = 0.0
        for j in range(N):
            if i == j:
                continue
            temp += -(x[2 * N + i] - x[2 * N + j]) / (
                pl.sqrt((x[2 * N + i] - x[2 * N + j]) ** 2 + (x[3 * N + i] - x[3 * N + j]) ** 2) ** 3
            )
        xdot = pl.append(xdot, temp)
    for i in range(N):
        temp = 0.0
        for j in range(N):
            if i == j:
                continue
            temp += -(x[3 * N + i] - x[3 * N + j]) / (
                pl.sqrt((x[2 * N + i] - x[2 * N + j]) ** 2 + (x[3 * N + i] - x[3 * N + j]) ** 2) ** 3
            )
        xdot = pl.append(xdot, temp)
    for i in range(N):
        xdot = pl.append(xdot, x[i])
    for i in range(N):
        xdot = pl.append(xdot, x[N + i])

    print("len xdot is: " + str(len(xdot)))
    return xdot

def data_to_ch(data):
    ch = {}
    for ch_ind in range(1, 97):
        ch[ch_ind] = {}
        ch[ch_ind]["bl"] = data[ch_ind]["blanks"]
        ch[ch_ind]["bl_mu"] = pl.mean(ch[ch_ind]["bl"])
        ch[ch_ind]["bl_sem"] = pl.std(ch[ch_ind]["bl"]) / pl.sqrt(len(ch[ch_ind]["bl"]))
        for ind in sorted(data[ch_ind].keys()):
            if ind != "blanks":
                k = ind[0]
                if k not in ch[ch_ind]:
                    ch[ch_ind][k] = {}
                    ch[ch_ind][k]["fr"] = []
                    ch[ch_ind][k]["fr_mu"] = []
                    ch[ch_ind][k]["fr_sem"] = []
                    ch[ch_ind][k]["pos_y"] = []
                    ch[ch_ind][k]["dprime"] = []
                ch[ch_ind][k]["fr"].append(data[ch_ind][ind]["on"])
                ch[ch_ind][k]["fr_mu"].append(pl.mean(data[ch_ind][ind]["on"]))
                ch[ch_ind][k]["fr_sem"].append(pl.std(data[ch_ind][ind]["on"]) / pl.sqrt(len(data[1][ind]["on"])))
                ch[ch_ind][k]["pos_y"].append(ind[2])
                # print ch[ch_ind][k]['pos_y']
                # print pl.std(data[ch_ind][ind]['on'])
                ch[ch_ind][k]["dprime"].append(
                    (pl.mean(data[ch_ind][ind]["on"]) - ch[ch_ind]["bl_mu"])
                    / ((pl.std(ch[ch_ind]["bl"]) + pl.std(data[ch_ind][ind]["on"])) / 2)
                )
                # print ch[ch_ind]['OSImage_5']['pos_y']
    return ch

def ICcircle_to_ICeuklid(IC):
    """
    converts from IC in cirle parameters to IC in euklidean space (rotational invariant).
    The formats are:
    IC_euklid: [x, y, z, vx, vy, vz]
    IC_circle: [y, vy, |v|, |l|, phiv], where |v| is the magnitude of CoM velocity, |l| 
        is the distance from leg1 (assumed to be at [0,0,0]) to CoM, and phiv the angle
        of the velocity in horizontal plane wrt x-axis
    *NOTE* for re-conversion, the leg position is additionally required, assumed to be [0,0,0]
    Further, it is assumed that the axis foot-CoM points in x-axis
    
    :args:
        IC (5x float): the initial conditions in circular coordinates

    :returns:
        IC (6x float): the initial conditions in euklidean space
    
    """
    y, vy, v, l, phiv = IC
    z = 0
    xsq = l**2 - y**2
    if xsq < 0:
        raise RuntimeError('Error in initial conditions: y > l!')
    x = -sqrt(xsq)
    vhsq = v**2 - vy**2
    if vhsq < 0:
        raise RuntimeError('Error in initial conditions: |vy| > |v|!')
    v_horiz = sqrt(vhsq)
    vx = v_horiz * cos(phiv)
    #vz = v_horiz * sin(phiv)
    vz = v_horiz * sin(phiv)
    return [x, y, z, vx, vy, vz]

def measureDistance(lat1, lon1, lat2, lon2):
    R = 6383.137 # Radius of earth at Chajnantor aprox. in KM
    dLat = (lat2 - lat1) * np.pi / 180.
    dLon = (lon2 - lon1) * np.pi / 180.
    a = pl.sin(dLat/2.) * pl.sin(dLat/2.) + pl.cos(lat1 * np.pi / 180.) * pl.cos(lat2 * np.pi / 180.) * pl.sin(dLon/2.) * pl.sin(dLon/2.)
    c = 2. * atan2(pl.sqrt(a), pl.sqrt(1-a))
    d = R * c
    return d * 1000. # meters

def exhaustVelocity( energyDensity=None, pressureOut=None, pressureIn=const_atmosphericPressure*25, temperature=3500, kappa=1.666666, molarMass=2 ):
	if (energyDensity!=None):	# from kinetic energy
		return pylab.sqrt( 2*energyDensity )
	else:	 # lavalNozzle adiabatic expansion http://en.wikipedia.org/wiki/De_Laval_nozzle
		gamma = (kappa)/(kappa-1)
		ExpansionTerm = 1.0
		if (( pressureOut != None ) ):
				ExpansionTerm = (1.0 - (pressureOut/pressureIn)**(1/gamma) )
		return pylab.sqrt(2000.0*const_universalGas*temperature*gamma*ExpansionTerm/molarMass)			

 def pkj1pk(self, r, z=0):
     """
       integral of the power spectrum*k over j1
     """
     return (
         M.sqrt(M.pi / 2.0)
         / r
         * self.besselInt.besselInt(lambda k: self.delta(k / r, z) / M.sqrt(k), 1.5, self.besselN, self.besselh)
     )

def xyamb(xytab,qu,xyout=''):

    mytb=taskinit.tbtool()

    if not isinstance(qu,tuple):
        raise Exception,'qu must be a tuple: (Q,U)'

    if xyout=='':
        xyout=xytab
    if xyout!=xytab:
        os.system('cp -r '+xytab+' '+xyout)

    QUexp=complex(qu[0],qu[1])
    print 'Expected QU = ',qu   # , '  (',pl.angle(QUexp)*180/pi,')'

    mytb.open(xyout,nomodify=False)

    QU=mytb.getkeyword('QU')['QU']
    P=pl.sqrt(QU[0,:]**2+QU[1,:]**2)

    nspw=P.shape[0]
    for ispw in range(nspw):
        st=mytb.query('SPECTRAL_WINDOW_ID=='+str(ispw))
        if (st.nrows()>0):
            q=QU[0,ispw]
            u=QU[1,ispw]
            qufound=complex(q,u)
            c=st.getcol('CPARAM')
            fl=st.getcol('FLAG')
            xyph0=pl.angle(pl.mean(c[0,:,:][pl.logical_not(fl[0,:,:])]),True)
            print 'Spw = '+str(ispw)+': Found QU = '+str(QU[:,ispw])  # +'   ('+str(pl.angle(qufound)*180/pi)+')'
            #if ( (abs(q)>0.0 and abs(qu[0])>0.0 and (q/qu[0])<0.0) or
            #     (abs(u)>0.0 and abs(qu[1])>0.0 and (u/qu[1])<0.0) ):
            if ( pl.absolute(pl.angle(qufound/QUexp)*180/pi)>90.0 ):
                c[0,:,:]*=-1.0
                xyph1=pl.angle(pl.mean(c[0,:,:][pl.logical_not(fl[0,:,:])]),True)
                st.putcol('CPARAM',c)
                QU[:,ispw]*=-1
                print '   ...CONVERTING X-Y phase from '+str(xyph0)+' to '+str(xyph1)+' deg'
            else:
                print '      ...KEEPING X-Y phase '+str(xyph0)+' deg'
            st.close()
    QUr={}
    QUr['QU']=QU
    mytb.putkeyword('QU',QUr)
    mytb.close()
    QUm=pl.mean(QU[:,P>0],1)
    QUe=pl.std(QU[:,P>0],1)
    Pm=pl.sqrt(QUm[0]**2+QUm[1]**2)
    Xm=0.5*atan2(QUm[1],QUm[0])*180/pi

    print 'Ambiguity resolved (spw mean): Q=',QUm[0],'U=',QUm[1],'(rms=',QUe[0],QUe[1],')','P=',Pm,'X=',Xm

    stokes=[1.0,QUm[0],QUm[1],0.0]
    print 'Returning the following Stokes vector: '+str(stokes)
    
    return stokes