def estrampcpx(inData, estData=None):
    if estData is not None:
        if np.any(np.iscomplex(estData)):
            estData=estData.conj()        
            inData=inData*estData
        else:
            estData=np.exp(-1.j*estData)
            inData=inData*estData
    indx=np.log10(abs(np.fft.fftshift(np.fft.fft2(np.angle(inData))))).argmax() #to display fft you may want to take log10 of the value here
    subi=np.unravel_index(indx, inData.shape)
    X,Y=np.meshgrid(np.r_[0:inData.shape[1]], np.r_[0:inData.shape[0]])
    fxmax=inData.shape[1]/2.
    fymax=inData.shape[0]/2.
    surface=np.pi*( (subi[0]-fymax)*Y/fymax + (subi[1]-fxmax)*X/fxmax )
    #print [ subi, fxmax, fymax ]
    outData=inData*np.exp(-1.j*surface)
    outVal=[(subi[0]-fymax)/fymax , (subi[1]-fxmax)/fxmax]
    #test if we get the direction right.
    indx2=np.log10(abs(np.fft.fftshift(np.fft.fft2(np.angle(outData))))).argmax() 
    subi2=np.unravel_index(indx2, inData.shape)
    dist1=np.sqrt( (subi[0]-fymax)**2. + (subi[1]-fxmax)**2.)
    dist2=np.sqrt( (subi2[0]-fymax)**2. + (subi2[1]-fxmax)**2.)
    #print [dist1, dist2]
    if dist1<dist2:
        #we got the direction wrong
        outVal=[-(subi[0]-fymax)/fymax , -(subi[1]-fxmax)/fxmax]
    #outData=np.exp(-1.j*surface)
    return outVal

def find_correct_sign(z1,z2,z_approx):
    '''
    Create new vector from z1, z2 choosing elements with sign matching z_approx
    
    This is used when you have to make a root choice on a complex number.
    and you know the approximate value of the root. 
    
    .. math:: 
        
        z1,z2 = \\pm \\sqrt(z^2)
        

    Parameters
    ------------
    z1 : array-like
        root 1
    z2 : array-like
        root 2
    z_approx : array-like
        approximate answer of z
    
    Returns 
    ----------
    z3 : npy.array
        array built from z1 and z2 by 
        z1 where sign(z1) == sign(z_approx), z2 else
    
    '''
    return npy.where(
    npy.sign(npy.angle(z1)) == npy.sign(npy.angle(z_approx)),z1, z2)    

def main():
    sample_rate = 128
    N = 128
    freq = 10.5
    k = np.pi * 2 * freq / sample_rate
    x = np.sin(np.arange(N, dtype='float64') * k)
    #x = np.pad(x, (0, N * 3), mode='constant')
    xfreq = np.arange(N/2) * (float(sample_rate) / len(x))

    fv = np.fft.fft(x)
    amp = np.abs(fv[:N/2])
    phase = np.angle(fv[:N/2])
    print amp
    print phase
    mid = fv[N/2-3:N/2+4]
    print amp[0], np.abs(mid), mid
    plt.figure()
    plt.subplot(2, 1, 1)
    plt.plot(xfreq, amp)
    plt.subplot(2, 1, 2)
    plt.plot(xfreq, phase)

    N = len(x)
    fv[:N/2] *= np.exp(np.random.uniform(0, np.pi*2, N/2) * 1j)
    fv[N/2:] = np.conjugate(fv[:N/2])
    recon = np.fft.ifft(fv)
    plt.figure()
    plt.subplot(2, 1, 1)
    plt.plot(xfreq, np.angle(fv[:len(xfreq)]))
    plt.subplot(2, 1, 2)
    plt.plot(np.real(recon))
    plt.show()

def acumPassband( nacum=5, outfile=passbandFile ):
  [mpNameList, tsysNameList] = makeMpNameLists( )
  vectorSumPbMatrix = numpy.zeros( [16,15], dtype=complex )    # vector sum of passband values
  scalarSumPbMatrix = numpy.zeros( [16,15], dtype=float )	   # scalar (magnitude) passband sum
  npts = 0
  while (npts < nacum) :
    waitForNewData( )
    visMatrix = getVisMatrix( mpNameList )
    vectorSumPbMatrix += visMatrix
    scalarSumPbMatrix += numpy.abs(visMatrix)
    npts += 1
    print "accumulating passband (%d/%d records)" % (npts,nacum)
    meritMatrix = numpy.divide( numpy.abs(vectorSumPbMatrix), scalarSumPbMatrix )   # mag(vector) / sum(mags) 
      # ... some values may be nan if data are missing
    print numpy.array_str( meritMatrix, precision=3, max_line_width=200 )
    maskedMeritMatrix = numpy.ma.array( meritMatrix, mask=(numpy.isnan(meritMatrix) ) )
      # ... mask off the nans
    merit = numpy.ma.mean( maskedMeritMatrix ) 
    print "merit = %.3f" % merit

  # for consistency with mfcal, rotate passband phases to make phase(win1) = 0 for each antenna
  print "\nbefore rotation:"
  print numpy.array_str( numpy.angle(vectorSumPbMatrix, deg=True ), precision=0, suppress_small=True, max_line_width=200 )
  for n in range(0,15) :
    Win1 = vectorSumPbMatrix[0][n]
    if numpy.abs(Win1) > 0. :
      for m in range(0,16) :
        vectorSumPbMatrix[m][n] = vectorSumPbMatrix[m][n] * numpy.conj(Win1)
          # ... don't bother normalizing, as this happens at the very end anyway
  print "\nafter rotation:"
  print numpy.array_str( numpy.angle(vectorSumPbMatrix, deg=True ), precision=0, suppress_small=True, max_line_width=200 )

  print "\nsaving normalized passband to file %s" % outfile
  numpy.save(outfile, vectorSumPbMatrix/numpy.abs(vectorSumPbMatrix) )

def change_parameterization(g, p_from, p_to):
    """

    """
    if p_from == p_to:
        return g

    elif p_from == "MAT" and p_to == "C":
        theta = np.arctan2(g[..., 1, 0], g[..., 0, 0])
        return np.exp(1j * theta)
    elif p_from == "MAT" and p_to == "ANG":
        return np.arctan2(g[..., 1, 0], g[..., 0, 0])
    elif p_from == "C" and p_to == "MAT":
        theta = np.angle(g)
        c = np.cos(theta)
        s = np.sin(theta)
        return np.array([[c, -s], [s, c]]).transpose(range(2, 2 + c.ndim) + [0, 1])
    elif p_from == "C" and p_to == "ANG":
        return np.angle(g)
    elif p_from == "ANG" and p_to == "MAT":
        c = np.cos(g)
        s = np.sin(g)
        return np.array([[c, -s], [s, c]]).transpose(range(2, 2 + c.ndim) + [0, 1])
    elif p_from == "ANG" and p_to == "C":
        return np.exp(1j * g)
    else:
        raise ValueError("Unsupported conversion:" + str(p_from) + " to " + str(p_to))

def makeGnuFig(filename):
    resultmat = np.zeros([len(squid.xaxis), 9])
    S11 = elem.S11
    ydat = measdata.ydat
    resultmat[:, 0] = squid.xaxis
    resultmat[:, 1] = ydat.real
    resultmat[:, 2] = S11.real
    resultmat[:, 3] = ydat.imag
    resultmat[:, 4] = S11.imag
    resultmat[:, 5] = abs(ydat)
    resultmat[:, 6] = abs(S11)
    resultmat[:, 7] = np.unwrap(np.angle(S11), discont=pi)
    resultmat[:, 8] = np.unwrap(np.angle(ydat), discont=pi)
    np.savetxt(filename, resultmat, delimiter='\t')
    # Plot in Gnuplot
    g1 = gp.Gnuplot(persist=1, debug=1)
    g1("plot '" + str(filename) + "' u 1:2 w l t 'Meas.real'")
    g1("replot '" + str(filename) + "' u 1:3 w l t 'Fit.real'")
    g1("replot '" + str(filename) + "' u 1:4 w l t 'Meas.imag'")
    g1("replot '" + str(filename) + "' u 1:5 w l t 'Fit.imag'")
    g2 = gp.Gnuplot(persist=1, debug=1)
    g2("plot '" + str(filename) + "' u 1:6 w l t 'Meas.mag'")
    g2("replot '" + str(filename) + "' u 1:7 w l t 'Fit.mag'")
    g2("replot '" + str(filename) + "' u 1:8 w l t 'Meas.phase'")
    g2("replot '" + str(filename) + "' u 1:9 w l t 'Fit.phase'")
    return

    def updateFigs(ZarcFitWindow):

        if ZarcFitWindow.radioButtonSerial.isChecked():
            Z = ZarcFitWindow.zarc.Zseries(ZarcFitWindow.frequency)  
        elif ZarcFitWindow.radioButtonParallel.isChecked():
            Z = ZarcFitWindow.zarc.Zparallel(ZarcFitWindow.frequency)  
        else:
            Exception("Not implemented!! choose either series or parallel")

        ZarcFitWindow.lineCole.set_data(Z.real, Z.imag)
        ZarcFitWindow.lineColeobs.set_data(ZarcFitWindow.obs.real, ZarcFitWindow.obs.imag)
        ZarcFitWindow.axCole.draw_artist(ZarcFitWindow.axCole.patch)        
        ZarcFitWindow.axCole.draw_artist(ZarcFitWindow.lineCole)
        ZarcFitWindow.axCole.draw_artist(ZarcFitWindow.lineColeobs)


        ZarcFitWindow.lineBodeMagn.set_ydata(abs(Z))
        ZarcFitWindow.lineBodeMagnobs.set_ydata(abs(ZarcFitWindow.obs))
        ZarcFitWindow.axBodeMagn.draw_artist(ZarcFitWindow.axBodeMagn.patch)
        ZarcFitWindow.axBodeMagn.draw_artist(ZarcFitWindow.lineBodeMagn)
        ZarcFitWindow.axBodeMagn.draw_artist(ZarcFitWindow.lineBodeMagnobs)        

        ZarcFitWindow.lineBodePhase.set_ydata(abs(np.angle(Z, deg=True)))
        ZarcFitWindow.lineBodePhaseobs.set_ydata(abs(np.angle(ZarcFitWindow.obs, deg=True)))
        ZarcFitWindow.axBodePhase.draw_artist(ZarcFitWindow.axBodePhase.patch)
        ZarcFitWindow.axBodePhase.draw_artist(ZarcFitWindow.lineBodePhase)        
        ZarcFitWindow.axBodePhase.draw_artist(ZarcFitWindow.lineBodePhaseobs)  

        ZarcFitWindow.figCole.canvas.update()

 def plot_line_on_signal(self,   color ='m',
                                                 sampling_rate = None,
                                                 **kargs):
     if 'ax'in kargs:
         ax = kargs['ax']
     
     if sampling_rate is None:
         x = self.time_line
         v = self.value_line
         y = np.cos(np.angle(v))*np.abs(v)
     else :
         if self.time_line.size>1:
             old_dt = (self.time_line[1]-self.time_line[0]).rescale('s').magnitude
             x = np.arange(self.time_start, self.time_stop+old_dt, 1./sampling_rate.rescale('Hz').magnitude)
             #~ l = int((self.time_stop-self.time_start)*sampling_rate.rescale('Hz').magnitude)
             #~ x = self.time_start + np.arange(l) / sampling_rate.rescale('Hz').magnitude
             
         else:
             x=self.time_line
         v = self.value_line
         y = np.cos(np.angle(v))*np.abs(v)
         
         # Before resampling, in order to avoid slow down due the use of ifft in scipy.resample
         # y is padded with 0 proportionnally to the distance from x.size to the next 2**N 
         # QUESTION: does it lead to some strange edge effects???
         N=np.ceil(np.log2(x.size))
         yy=np.r_[y,np.zeros(np.floor(y.size*(2**N-x.size)/x.size))]
         yy = scipy.signal.resample( yy, 2**N)
         y = yy[:x.size]
     
     l = ax.plot(x,y  , linewidth = 1, color=color)
     return l

    def write_cube(self, index, fname, repeat=None, real=True):
        """Dump specified Wannier function to a cube file"""
        from ase.io.cube import write_cube

        # Default size of plotting cell is the one corresponding to k-points.
        if repeat is None:
            repeat = self.kptgrid
        atoms = self.calc.get_atoms() * repeat
        func = self.get_function(index, repeat)

        # Handle separation of complex wave into real parts
        if real:
            if self.Nk == 1:
                func *= np.exp(-1.0j * np.angle(func.max()))
                if 0:
                    assert max(abs(func.imag).flat) < 1e-4
                func = func.real
            else:
                func = abs(func)
        else:
            phase_fname = fname.split(".")
            phase_fname.insert(1, "phase")
            phase_fname = ".".join(phase_fname)
            write_cube(phase_fname, atoms, data=np.angle(func))
            func = abs(func)

        write_cube(fname, atoms, data=func)

def plot_parameters(gid, data):
    print("Plotting transformed wavepacket parameters of group '"+str(gid)+"'")

    # Grid of mother and first spawned packet
    grid_m = data[0][0]
    grid_s = data[1][0]

    # Parameters of mother and first spawned packet
    P, Q, S, p, q = data[0][1]
    B, A, S, b, a = data[1][1]

    X = P*abs(Q)/Q

    # Various interesting figures

    fig = figure()
    ax = fig.gca()

    ax.plot(grid_m, real(X), "*", label=r"$\Re \frac{P |Q|}{Q}$")
    ax.plot(grid_s, real(B), "o", label=r"$\Re B$")

    ax.legend()
    ax.grid(True)
    fig.savefig("test_spawned_PI_realparts_group"+str(gid)+GD.output_format)



    fig = figure()
    ax = fig.gca()

    ax.plot(grid_m, imag(X), "*", label=r"$\Im \frac{P |Q|}{Q}$")
    ax.plot(grid_s, imag(B), "o", label=r"$\Im B$")

    ax.legend()
    ax.grid(True)
    fig.savefig("test_spawned_PI_imagparts_group"+str(gid)+GD.output_format)



    fig = figure()
    ax = fig.gca()

    ax.plot(real(X), imag(X), "-*", label=r"traject $\frac{P |Q|}{Q}$")
    ax.plot(real(B), imag(B), "-*", label=r"traject $B$")

    ax.legend()
    ax.grid(True)
    fig.savefig("test_spawned_PI_complex_trajectories_group"+str(gid)+GD.output_format)



    fig = figure()
    ax = fig.gca()

    ax.plot(grid_m, angle(X), label=r"$\arg \frac{P |Q|}{Q}$")
    ax.plot(grid_s, angle(B), label=r"$\arg B$")

    ax.legend()
    ax.grid(True)
    fig.savefig("test_spawned_PI_angles_group"+str(gid)+GD.output_format)

def stretch_lattice(x, y, fov, stretch, angle):
    if stretch < 1:
        raise RunTimeError("Error. Stretch values must be >1.  {} entered.".format(stretch))
    print "Stretching by a factor of {} at an angle of {} degrees to the x-axis about the image center".format(
        stretch, angle
    )
    x_new = []
    y_new = []
    x_shifted, y_shifted = x - fov / 2.0, y - fov / 2.0
    x = (
        x_shifted
        + np.sqrt(x_shifted ** 2 + y_shifted ** 2)
        * np.cos(np.angle(x_shifted + y_shifted * 1j) - np.radians(angle))
        * np.cos(np.radians(angle))
        * (stretch - 1)
        + fov / 2.0
    )
    y = (
        y_shifted
        + np.sqrt(x_shifted ** 2 + y_shifted ** 2)
        * np.cos(np.angle(x_shifted + y_shifted * 1j) - np.radians(angle))
        * np.sin(np.radians(angle))
        * (stretch - 1)
        + fov / 2.0
    )
    for i in range(len(x)):
        if x[i] < fov and x[i] >= 0 and y[i] < fov and y[i] >= 0:
            x_new.append(x[i])
            y_new.append(y[i])
    print "{} particles discarded".format(len(x) - len(x_new))
    return np.array(x_new), np.array(y_new)

def update(attrname,old,new):
    ori=orientation.value*(numpy.pi/180.0)
    oriWidth=orientationWidth.value*(numpy.pi/180.0)
    oband=[ori, ori-oriWidth]
    #oband=[ori, numpy.pi/8]
    t1 = numpy.angle(numpy.exp(complex(0,1)*(theta-oband[0]))) # center it
    t1 = numpy.exp(-(t1**2/(2*oband[1]*oband[1]))) # and put a gaussian around it
    # make the 2nd lobe on the other side of k-space
    t2 = numpy.angle(numpy.exp(complex(0,1)*(theta-oband[0]-numpy.pi)))
    t2 = numpy.exp(-(t2**2/(2*oband[1]*oband[1])))

    thisSF=spatialFreq.value
    sfWidth=spatialFreqWidth.value
    fband=[thisSF,thisSF-sfWidth]
    SF = numpy.exp(-((r-fband[0])**2/(2.0*fband[1]*fband[1])))  

    customFilter = SF*(t1 + t2)
    source2 = ColumnDataSource(data={'image': [customFilter]})
    p2.image(image="image", x=[0], y=[0], dw=[10], dh=[10],source=source2,palette="Spectral11")
    filteredFFT = imageFFT*customFilter
    filteredImage = numpy.fft.ifft2(numpy.fft.fftshift(filteredFFT))
    filteredImageReal=numpy.real(filteredImage)
    source3 = ColumnDataSource(data={'image': [filteredImageReal]})
    p3.image(image="image", x=[0], y=[0], dw=[10], dh=[10],source=source3)
    thisTitle="Ori: %2.1f [%2.1f, %2.1f], SF: %2.1f, %2.1f"%(orientation.value,oband[0]*180/numpy.pi,
                oband[1]*180/numpy.pi,thisSF, sfWidth)
    p3.title=thisTitle

    def write(self, MA=False):

        with open(self.filename, 'w') as g:
            g.write( "! TOUCHSTONE file generated by Touchstone Python Converter\n" )
            g.write( "! Date and time: %s\n"%time.strftime("%b %d %Y %H:%M:%S") )
            if MA:
                g.write( "# GHZ S MA R 50\n" )
                for row in zip(self.raw['freq'], self.raw['s11'], self.raw['s12'], self.raw['s21'], self.raw['s22']):
                    freq, s11, s12, s21, s22 = row
                    freq *= 1e-9
                    ms11, as11 = abs(s11), angle(s11,deg=True)
                    ms12, as12 = abs(s12), angle(s12,deg=True)
                    ms21, as21 = abs(s21), angle(s21,deg=True)
                    ms22, as22 = abs(s22), angle(s22,deg=True)
                    pattern = "{:<15.6e}{:>15.6e}{:>15.6e}{:>15.6e}{:>15.6e}{:>15.6e}{:>15.6e}{:>15.6e}{:>15.6e}\n"
                    g.write( pattern.format(freq, ms11, as11, ms21, as21, ms12, as12, ms22, as22) )
            else:
                g.write( "# GHZ S RI R 50\n" )
                for row in zip(self.raw['freq'], self.raw['s11'], self.raw['s12'], self.raw['s21'], self.raw['s22']):
                    freq, s11, s12, s21, s22 = row
                    freq *= 1e-9
                    pattern = "{:<15.6e}{:>15.6e}{:>15.6e}{:>15.6e}{:>15.6e}{:>15.6e}{:>15.6e}{:>15.6e}{:>15.6e}\n"
                    g.write( pattern.format(freq, s11.real, s11.imag, s21.real, s21.imag, s12.real, s12.imag, s22.real, s22.imag) )
           
            g.write( "! noise parameters\n" )
            for row in zip(self.raw['freq'], self.raw['nfmin'], self.raw['sopt'], self.raw['rn']):
                freq, nfmin, sopt, rn = row
                freq *= 1e-9
                msopt, asopt = abs(sopt), angle(sopt,deg=True)
                rn *= 1./50.
                pattern = "{:<15.6e}{:>15.6e}{:>15.6e}{:>15.6e}{:>15.6e}\n"
                #print freq, nfmin, msopt, asopt, rn
                g.write( pattern.format(freq, nfmin, msopt, asopt, rn) )

def istft(spec, hop):
    """
    :param spec: The result of an stft, spectra at different times
    :param hop: number of smaples to move forward between windows
    :return: time domain signal which is related to the original by a constant coefficient.

    fftlen inferred from shape of input vector.
    """
    fftlen = 2 * (spec.shape[1] - 1)
    xs = np.empty((spec.shape[0], fftlen))
    w = np.hanning(fftlen)

    out = np.empty(spec.shape[0] * hop + fftlen,
                   dtype=complex)  # Need a quick calc to find the length of the output array
    phase = np.zeros(spec.shape[1])
    adj = np.zeros(spec.shape[1], dtype=complex)

    for t in range(spec.shape[0]):
        if t > 0:
            phase += np.angle(spec[t]) - np.angle(spec[t - 1])
            phase %= 2 * np.pi
            adj.real, adj.imag = np.cos(phase), np.sin(phase)

        # xs[t] = np.fft.ifft(Xs[t])
        xs[t] = np.fft.irfft(np.abs(spec[t]) * adj)
        out[t * hop: t * hop + fftlen] += xs[t] * w

    return out

def plot_vis_3(d1,d2,d3):
    fig = plt.figure()
    p1=fig.add_subplot(3,2,1)
    p1.set_title("data_amp")
    i1=p1.imshow(np.abs(d1), interpolation='nearest',aspect='auto')
    fig.colorbar(i1)
    p1.xaxis.set_ticklabels([])
    p2=fig.add_subplot(3,2,2)
    p2.set_title("data_phs")
    i2=p2.imshow(np.angle(d1),interpolation='nearest',aspect='auto')
    fig.colorbar(i2)
    p2.xaxis.set_ticklabels([])
    p3=fig.add_subplot(3,2,3)
    p3.set_title("model_amp")
    i3=p3.imshow(np.abs(d2),interpolation='nearest',aspect='auto')
    fig.colorbar(i3)
    p3.xaxis.set_ticklabels([])
    p4=fig.add_subplot(3,2,4)
    p4.set_title("model_phs")
    i4=p4.imshow(np.angle(d2),interpolation='nearest',aspect='auto')
    fig.colorbar(i4)
    p4.xaxis.set_ticklabels([])
    p5=fig.add_subplot(3,2,5)
    p5.set_title("sol*model_amp")
    i5=p5.imshow(np.abs(d3),interpolation='nearest',aspect='auto')
    fig.colorbar(i5)
    p6=fig.add_subplot(3,2,6)
    p6.set_title("sol*model_phs")
    i6=p6.imshow(np.angle(d3),interpolation='nearest',aspect='auto')
    fig.colorbar(i6)
    plt.show()

def get_closure_phase(infile='L401323_SB349_uv.dppp.MS',\
                 triangle = ['TS001','DE601HBA','DE605HBA']):
    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']='mytask'
    myf['taskname']='mytask'
    tb=myf['tb']
    oroot = infile.split('uv')[0]
    for lfile in np.sort(glob.glob(oroot+'*ms')):
        os.system('ms2uvfits in='+lfile+' out='+lfile.replace('ms','fits')+' writesyscal=F')
        if lfile == infile:
            continue
        tb.open(lfile+'/ANTENNA')
        names = tb.getcol('NAME')
        trnum = []
        for itr in range(3):
            trnum.append(np.argwhere(names==triangle[itr])[0][0])
        tb.close()
        trnum.sort()
        tb.open(lfile)
        ant1 = tb.getcol('ANTENNA1')
        ant2 = tb.getcol('ANTENNA2')
        data = tb.getcol('DATA')
        ph12 = +np.angle(data[0,0,(ant1==trnum[0])&(ant2==trnum[1])])
        ph23 = +np.angle(data[0,0,(ant1==trnum[1])&(ant2==trnum[2])])
        ph31 = -np.angle(data[0,0,(ant1==trnum[0])&(ant2==trnum[2])])
        clph = ph12+ph23+ph31
        np.putmask(clph,clph>np.pi,clph-2.*np.pi)
        np.putmask(clph,clph<-np.pi,clph+2.*np.pi)
#        np.savetxt(lfile.replace('ms','txt'),np.unwrap(clph))
        np.savetxt(lfile.replace('ms','txt'),clph)

def update_plots(widget, baselines):
    try:
        (corr_time, left, right, current, 
         total, lags, visibility, phase_fit, m, c) = correlator.get_correlation()
        baseline = left, right
        logger.debug('received baseline %s' % repr(baseline))
    except NoCorrelations:
        widget.after(1, update_plots, widget, baselines)
        return # it never comes to this
    if baseline not in baselines.keys():
        corr.axes.grid()
        #corr.axes.set_xlabel('Lag', size='large')
        #corr.axes.set_ylabel('Correlation Function', size='large')
        phase_line = corr.plot(f, angle(visibility), '%so' % colors[current%len(colors)], linewidth=1, label=repr(baseline))[0]
        fit_line = corr.plot(f, real(phase_fit), '%s-' % colors[current%len(colors)], linewidth=1, label=None)[0]
        baselines[baseline] = phase_line, fit_line
    else:
        corr.axes.legend()
        phase_line, fit_line = baselines[baseline]
        corr.update_line(phase_line, f, angle(visibility))
        corr.update_line(fit_line, f, real(phase_fit))
    if current == total-1:
        widget.update()
        logger.info('update in')
    widget.after_idle(update_plots, widget, baselines)

	def getPhaseNoise(self, psi0, psi1):
		"""
		Warning: this function considers spin distribution ellipse to be horizontal,
		which is not always so.
		"""

		ensembles = psi0.shape[0]
		get = self._env.fromDevice
		reduce = self._reduce
		creduce = self._creduce

		i = self._getInteraction(psi0, psi1)

		i = get(creduce(i, ensembles)) # Complex numbers {S_xj + iS_yj, j = 1..N}
		phi = numpy.angle(i) # normalizing

		# Center of the distribution can be shifted to pi or -pi,
		# making mean() return incorrect values.
		# The following approximate method will allow us to shift the center to zero
		# It will work only if the maximum of the distribution is clearly
		# distinguished; otherwise it can give anything as a result

		Pperp = numpy.exp(1j * phi) # transforming Pperp to distribution on the unit circle
		Pmean = Pperp.mean() # Center of masses is supposed to be close to the center of distribution

		# Normalizing the direction to the center of masses
		# Now angle(Pmean) ~ proper mean of Pperp
		Pmean /= numpy.abs(Pmean)

		# Shifting the distribution
		Pcent = Pperp * Pmean.conj()
		phi_centered = numpy.angle(Pcent)

		return phi_centered.std()

def subsample(args):
    """
    Rebin / Congrid variant
    """
    arr, shape, mode = args  # unpack arguments

    if mode == 'phase' and np.iscomplexobj(arr):
        arr = np.angle(arr)
    if np.iscomplexobj(arr):
        arr = np.abs(arr)

    if arr.shape == shape:
        return arr

    oshap = arr.shape

    for d in range(arr.ndim):
        n1 = shape[d]
        n2 = oshap[d]
        if n1 < n2:
            s = list(arr.shape)
            s.insert(d + 1, n2 // n1)
            s[d] = n1
            if mode == 'phase':
                arr = np.angle(np.exp(1j * arr.reshape(s)).mean(d + 1))
            elif mode == 'lables':
                arr = np.take(arr.reshape(s), 0, d + 1)
            else:
                arr = arr.reshape(s).mean(d + 1)
        else:
            arr = arr.repeat(n1 // n2, axis=d)
    return arr

def drawDataCallback(baseline):
    matplotlib.pyplot.clf()
    acc_n,interleave_a,interleave_b = get_data(baseline)

    matplotlib.pyplot.subplot(211)
    if ifch == True:
        matplotlib.pyplot.semilogy(numpy.abs(interleave_a))
        matplotlib.pyplot.xlim(0,1024)
    else:
        matplotlib.pyplot.semilogy(xaxis,numpy.abs(interleave_a))
    matplotlib.pyplot.grid()
    matplotlib.pyplot.title('Integration number %i \n%s'%(acc_n,baseline))
    matplotlib.pyplot.ylabel('Power (arbitrary units)')

    matplotlib.pyplot.subplot(212)
    if ifch == True:
        matplotlib.pyplot.plot(numpy.unwrap(numpy.angle(interleave_b)))
        matplotlib.pyplot.xlim(0,1024)
        matplotlib.pyplot.xlabel('FFT Channel')
    else:
        matplotlib.pyplot.plot(xaxis,(numpy.angle(interleave_b)))
        matplotlib.pyplot.xlabel('FFT Frequency')
    matplotlib.pyplot.ylabel('Phase')
    matplotlib.pyplot.ylim(-numpy.pi,numpy.pi)
    matplotlib.pyplot.grid()
    
    matplotlib.pyplot.draw()
    fig.canvas.manager.window.after(100, drawDataCallback,baseline)