def ICDPLA(alignment_matrix, path_trace_matrix, superior, inferior, row_offset, transposed=False):
    # stopping condition -> if max value is 0, there are no more alignments..or the matrix is empty, so just return nothing 
    if len(alignment_matrix) == 0 or (np.max(alignment_matrix) == 0):
        return 0, []
    # find max in Smith-Waterman matrix
    max_index = pl.unravel_index(alignment_matrix.argmax(), alignment_matrix.shape)
    # calculate distance over alignment path
    local_distance = 0.0
    # the start index will be where the alignment begins (necessary to remove submatrix)
    start_x = row_offset
    start_y = 0
    # these indices trace path backwards
    x_index = max_index[0]+row_offset
    y_index = max_index[1]
    while (x_index != -1 and x_index >= row_offset):
        start_x = x_index
        start_y = y_index
        local_distance += (np.linalg.norm(superior[x_index][:]-inferior[y_index][:])) / (4*np.sqrt(2))
        if transposed:
            [y_index, x_index] = path_trace_matrix[x_index][y_index]
        else:
            [x_index, y_index] = path_trace_matrix[x_index][y_index]
    # remove appropriate rows from sequence1 and split into two matrices to be involved in the same process
    alignment_top_submatrix = alignment_matrix[:(start_x-row_offset),:]
    alignment_bottom_submatrix = alignment_matrix[max_index[0]+1:, :]
    [distance_top, alignments_top] = ICDPLA(alignment_top_submatrix, path_trace_matrix, superior, inferior, row_offset, transposed)
    [distance_bottom, alignments_bottom] = ICDPLA(alignment_bottom_submatrix, path_trace_matrix, superior, inferior, row_offset+max_index[0], transposed)
    total_distance = distance_top+distance_bottom+local_distance
    alignments = [[[start_x, start_y], [max_index[0]+row_offset, max_index[1]]]]
    alignments.extend(alignments_top)
    alignments.extend(alignments_bottom)
    return total_distance, alignments

 def contour(self,tmin,tmax,t0min=-1,t0max=-1,steps=100.,amax=3000.,showopt=True,statusreport=False):
     if t0min<0:
         t0min=tmin
     if t0max<0:
         t0max=tmax    
     nbar=self.nbar
     sideband=self.sideband
     delta = self.delta
     nu = self.nu
     omega = self.omega
     Z=list(np.zeros((steps,steps)))
     z=list(np.zeros(steps))
     xdelta = (tmax-tmin)/steps
     ydelta = (t0max-t0min)/steps
     x = np.arange(tmin, tmax, xdelta)
     y = np.arange(t0min, t0max, ydelta)
     X, Y = np.meshgrid(x, y)
     dyn=Sideband(nbar=nbar,sideband=sideband,delta = delta,nu = nu,omega = omega,amax=amax)
     print 'generating list for analytical contour plot'
     half=False
     index=0
     for i in y:
         try:
             Z[index]=dyn.localsignal(i,x+i)
             z[index]=np.sum(Z[index])/steps
         except IndexError:
             print 'Warning: index ran out of range, removing last element'
             t=list(X)
             t.pop()
             X=np.array(t)
             t=list(Y)
             t.pop()
             Y=np.array(t)
         if i>(t0max-t0min)/2.+t0min and not half:
             print '    50% done'
             half=True
         if statusreport:
             print 'step '+str(index+1)+' of '+str(len(y))
         index=index+1
     print 'done'
     fig2=pyplot.figure()
     pyplot.title('Analytical local signal')
     pyplot.xlabel('t+t0')
     pyplot.ylabel('t0')
     pyplot.contourf(X, Y, Z)
     if showopt:
         m=pylab.unravel_index(np.array(Z).argmax(), np.array(Z).shape)
         pyplot.annotate('nbar = {:.2f}'.format(nbar)+'. Optimal t0 in plotted range is {:.6f}'.format(X[m])+'. Highest contrast is {:.2f}'.format(np.array(Z).max()),xy=(0.,-0.115), xycoords='axes fraction')
         pyplot.axvline(x=X[m],ls=':',color='k')
         pyplot.axhline(y=Y[m],ls=':',color='k')
     
     fig3=pyplot.figure()
     pyplot.title('Time-averaged local signal')
     pyplot.xlabel('t0')
     pyplot.plot(Y,z)        

 def __init__(self, parent,ident):
     self.parent = parent
     dataset, directory, index = ident
     dataX, dataY = self.getData(dataset, directory, index)
     Xmax = dataX.max()
     Xmin = dataX.min()
     center = dataX[pylab.unravel_index(np.array(dataY).argmax(),np.array(dataY).shape)]
     FWHM = (Xmax-Xmin)/6.0
     height = dataY.max()
     self.curveName = 'Lorentzian'
     self.parameterNames = ['FWHM', 'Center','Height', 'Offset']
     self.parameterValues = [FWHM,center,height,0.0]
     self.parameterFit = [True,True,True,True]

    def setPiTimeBoxes(self,which):
        # UPDATE Pi-TIME SPIN BOXES
        if self.curveName=='Rabi Flop':
            dataX, dataY = self.fitRabiflop.getData(self.dataset, self.directory, self.index)
            params = self.getParameter(which)
            detailedX = np.linspace(dataX.min(),dataX.max(),1000)
            dataY = self.fitRabiflop.fitFunc(detailedX, params)
            m=pylab.unravel_index(np.array(dataY).argmax(), np.array(dataY).shape)
            piTime=detailedX[m]
            self.TwoPiTimeBox.setValue(2.0*piTime)
            self.PiTimeBox.setValue(piTime)
            self.PiOverTwoTimeBox.setValue(piTime/2.0)

        if self.curveName in ['Cosine','Ramsey Fringes']:
            f=self.parent.savedAnalysisParameters[self.dataset, self.directory, self.index, self.curveName][which]['Frequency']
            self.TwoPiTimeBox.setValue(1.0/f*10**6)
            self.PiTimeBox.setValue(0.5/f*10**6)
            self.PiOverTwoTimeBox.setValue(0.25/f*10**6)

flop_fit_y_axis = evo.state_evolution_fluc(flop_x_axis, nb(), f_Rabi(), delta(),delta_fluc())

#red_chi2 = chi_square(flop_y_axis[fitting_region], flop_fit_y_axis[fitting_region], flop_errors[fitting_region], True,len(fit_params))
figure = pyplot.figure()

i=0
for par in fit_params:
    print 'P[{}] = {} +- {}'.format(i,par(),np.sqrt(cov[i][i]))
    i+=1
    
enb = np.sqrt(cov[0][0])
ef_Rabi = np.sqrt(cov[1][1])

#pyplot.plot(flop_x_axis*10**6,flop_fit_y_axis,'r-')
m=pylab.unravel_index(np.array(flop_fit_y_axis).argmax(), np.array(flop_fit_y_axis).shape)
#print 'Flop maximum at {:.2f} us'.format(flop_x_axis[m]*10**6)+' -> Expected optimal t0 at {:.2f} us'.format(flop_x_axis[m]/2.0*10**6)
#print 'Actual t0 = {}'.format(t0)
#print '2pi time {}'.format(flop_x_axis[m]*f_Rabi()*2.0)

#pyplot.plot(flop_x_axis*10**6,flop_y_axis, 'ro')
#pyplot.plot(deph_x_axis*10**6,deph_y_axis, 'bs')
pyplot.xlabel(r'Subsequent evolution time $\frac{\Omega t}{2\pi}$',fontsize=size*22)
pyplot.ylim((0,ymax))
pyplot.ylabel('Local Hilbert-Schmidt Distance',fontsize=size*22)
#pyplot.legend()

subseq_evolution=np.where(flop_x_axis>=t0)
nicer_resolution = np.linspace(t0,flop_x_axis.max(),1000)
deph_fit_y_axis = evo.deph_evolution_fluc(nicer_resolution, t0,nb(),f_Rabi(),delta(),delta_fluc())
flop_fit_y_axis = evo.state_evolution_fluc(nicer_resolution, nb(), f_Rabi(), delta(),delta_fluc())

for fluctuations in fluclist:
    nbarfitlist=[]

    for nb in nbarlist:
        nbar=Parameter(nb)
        omega_R=Parameter(omega_center)    
        
        flops=[]
        n1=int(n*100*fluctuations)
        for i in range(n1):
            print 'fluctuations = {:.2f}, nb = {:.2f}, i = {}'.format(fluctuations,nb,i)
            x=random.uniform(-1,1)*fluctuations
            sb=tp.Sideband(nb, sideband=sideband_order,omega=omega_center*(1.+x),nu=2.*np.pi*trap_frequency,amax=maxn)
            sb.anaplot(0, xmax*10**-6*sb.p.omega/(2.*np.pi), 50, 0, dephasing=False, discord=False, lsig=False)
            flops.append(sb.flop)
        
        flops=np.sum(flops,axis=0)/np.float32(n1)
        sb.x=2.*np.pi*sb.x/sb.p.omega
        m=pylab.unravel_index(np.array(flops).argmax(), np.array(flops).shape)   

        fitting_region = np.where(sb.x <= 2.*sb.x[m])
        p,success = fit(f, [nbar,omega_R], y = flops[fitting_region], x = sb.x[fitting_region])
        
        nbarfitlist.append(nbar())    
        
    pyplot.plot(nbarlist,nbarfitlist,label = 'Fitted nbars with {:.1%} intensity fluctuations'.format(fluctuations))

pyplot.xlabel('True nbar')
pyplot.ylabel('Apparent (fitted) nbar on carrier')
pyplot.legend(loc=2)
pyplot.show()

    def f(t):
        evolution = flop.state_evolution_fluc(t,nbar(), f_Rabi(), delta(), delta_fluc())
        return evolution
    fit_params=[]
    if info['fit_fRabi']: fit_params.append(f_Rabi)
    if info['fit_nbar']: fit_params.append(nbar)
    if info['fit_delta']: fit_params.append(delta)
    if info['fit_delta_fluc']: fit_params.append(delta_fluc)
    p,success = fit(f, fit_params, y = prob[fitting_region], x = times[fitting_region] - offset_time)
    print 'fit for f_Rabi is ', f_Rabi()
    print 'fit for nbar is', nbar()
    if 'plot_initial_values' in info and info['plot_initial_values']:
        evolution = flop.state_evolution_fluc( fit_times,fit_init_nbar, fit_init_fRabi,fit_init_delta,fit_init_delta_fluc )
    else:
        evolution = flop.state_evolution_fluc( fit_times, nbar(),f_Rabi(), delta(),delta_fluc())
    pi_time_arg = pylab.unravel_index(np.array(evolution).argmax(),np.array(evolution).shape)
    pi_time = fit_times[pi_time_arg]
    print 'nbar = {}'.format(nbar())
    print 'Rabi Pi Time is {} us'.format((pi_time)*10**6)
    print 'Rabi Pi/2 Time is {} us'.format((pi_time)/2.0*10**6)
    print 'Rabi Frequency is {} kHz'.format(f_Rabi()*10**-3)
    print "The detuning is centered around {} kHz and spreads with a variance of {} kHz".format(delta()*10**-3,np.abs(delta_fluc())*10**-3)
    plot_fit_label = 'fit with nb = {:.2f} and f_Rabi = {:.1f} kHz'.format(nbar(),10**-3 * f_Rabi())
    plot_data_label = 'measured data, sideband = {}'.format(sideband_order)

elif info['plot_type']=='ramsey_fringe':
    def ramsey_fringe(frequency,T2,phase,contrast,offset,t):
        return contrast*np.exp(-t/T2)*(np.cos(np.pi*frequency*t+phase)**2-.5)+.5+offset
    if 'fit_init_phase' in info: fit_init_phase=info['fit_init_phase']
    else: fit_init_phase=0
    if 'fit_init_contrast' in info: fit_init_contrast=info['fit_init_contrast']

    evolution = evo.rabiflop(nb(),omega(),sideband,x)
    return evolution

fitting_region = np.where(flop_x_axis <= xmax)
p,success = fit(f, [nb,omega], y = flop_y_axis[fitting_region], x = flop_x_axis[fitting_region])

figure = pyplot.figure()

print "nbar = ",nb()
print "Rabi Frequency (Driving strength) = ", omega()*10**(-3)/(2.*np.pi)," kHz"

# this calls tp to make the plots by the functions defined in TheoryPrediction.py (takes times in units of rabi frequency (driving strength)
sb=tp.Sideband(nb(), sideband=sideband,omega=omega(),nu=2.*np.pi*trap_frequency)
t0=dephasing_time
sb.anaplot(0, xmax*sb.p.omega/(2.*np.pi), 100, t0*sb.p.omega/(2.*np.pi), dephasing=True, discord=False, lsig=True)
m=pylab.unravel_index(np.array(sb.flop).argmax(), np.array(sb.flop).shape)
print 'Flop maximum at {:.2f} us'.format(sb.x[m]*10**6*2.*np.pi/sb.p.omega)+' -> Expected optimal t0 at {:.2f} us'.format(sb.x[m]*10**6*2.*np.pi/sb.p.omega/2.)
# rescale x-axis
sb.x=2.*np.pi*sb.x/sb.p.omega
pyplot.plot(sb.x*10**6,sb.flop)
pyplot.plot(sb.x*10**6,sb.deph)
#pyplot(sb.x,sb.flop)

pyplot.plot(flop_x_axis*10**6,flop_y_axis, '-o')
pyplot.plot(deph_x_axis*10**6,deph_y_axis, '-o')
pyplot.xlabel('t in us')
pyplot.ylim((0,1))
pyplot.ylabel('Population in the D-5/2 state')# + {0:.0f} kHz'.format(ymin))
#pyplot.legend()
pyplot.text(xmax*0.50*10**6,0.73, 'nbar = {:.2f}'.format(nb()))
pyplot.text(xmax*0.50*10**6,0.78, 'Rabi Frequency f = {:.2f} kHz'.format(omega()*10**(-3)/(2.*np.pi)))

 def makeplot(self,tmin,tmax,steps,nsteps=-1,dephasing=True,lsig=True,coh=True,discord=False,adiscord=False,ndiscord=False,t0=0,num=True,nmax=50.,amax=10000.,statusreport=False,plotgroundstatepop=False):
     nbar=self.nbar
     sideband=self.sideband
     delta = self.delta
     nu = self.nu
     omega = self.omega
     state = self.state 
     if discord:
         adiscord=ndiscord=True
     
     ymin=0
     ymax=1
     ysmin=0
     ysmax=.5
     fig1=pyplot.figure()
     dyn=Sideband(nbar,sideband=sideband,delta = delta,nu = nu,omega = omega,amax=amax,state=state,plotgroundstatepop=plotgroundstatepop)
     dyn.anaplot(tmin,tmax,steps,t0=t0,lsig=lsig,discord=adiscord,dephasing=dephasing)
     
     if coh and not num:
         print 'Warning: Coherences can only be plotted if numerical plot is used'
         
     if ndiscord and not num and not discord:
         print 'Warning: Numerical discord can only be plotted if numerical plot is used'
     
     if lsig and not dephasing:
         print 'Warning: Local signal is plotted without dephasing plot'
     
     numplot=0
     coplot=0
     lsplot=0
     discplot=0
     if num:
         numplot=1
         if coh:
             coplot=1
     if lsig:
         lsplot=1
     if adiscord or (ndiscord and num):
         discplot=1
     
     pyplot.subplot(1+lsplot+discplot,1+numplot+coplot,1)
     pyplot.title('Analytical, RWA, '+state)
     pyplot.plot(dyn.x,dyn.flop)
     pyplot.ylim( (ymin, ymax) )
     pyplot.xlim( (tmin, tmax) )
     m=pylab.unravel_index(np.array(dyn.flop).argmax(), np.array(dyn.flop).shape)
     print 'Flop maximum at {:.2f}'.format(dyn.x[m])+' -> Expected optimal t0 at {:.2f}'.format(dyn.x[m]/2.)
     if plotgroundstatepop:
         st='ground '
     else:
         st='excited '
     pyplot.ylabel(st+'state population')    
     s='nbar = {:.2f}'.format(nbar)+', amax = {:.1f}'.format(dyn.p.amax)+', sideband = {:.0f}'.format(sideband)
     if num:
         if nsteps==-1:
             nsteps=steps
         ndyn=numSideband(nbar,sideband=sideband,delta = delta,nu = nu,omega = omega,nmax=nmax,plotgroundstatepop=plotgroundstatepop)
         if state=='therm':
             r=ndyn.gibbsstate()
         elif state=='coh':
             r=ndyn.coherentstate()
         else:
             print 'Error, initial state not recognized.'
         s=s+', purity = {:.4f}'.format(np.trace(np.dot(r,r)))+', nmax = {0}'.format(nmax)
     if dephasing:
         pyplot.plot(dyn.x,dyn.deph)
         pyplot.axvline(x=t0,ls=':',color='k')
         s=s+', t0 = {0}'.format(t0)
     pyplot.annotate(s, xy=(0.,-0.1-1.3*(lsplot+discplot)), xycoords='axes fraction')
     if lsig:
         pyplot.subplot(1+lsplot+discplot,1+numplot+coplot,2+numplot+coplot)
         pyplot.plot(dyn.x,dyn.lsig)
         pyplot.ylim( (ysmin, ysmax) )
         pyplot.xlim( (tmin, tmax) )        
         pyplot.axvline(x=t0,ls=':',color='k')
         pyplot.ylabel('local signal')
     
     if adiscord:
         pyplot.subplot(1+lsplot+discplot,1+numplot+coplot,2+coplot+lsplot*(1+numplot+coplot)+numplot)
         pyplot.plot(dyn.x,dyn.disc)
         pyplot.ylabel('discord')
         pyplot.xlim( (tmin, tmax) )
         pyplot.axvline(x=t0,ls=':',color='k')  
                
     if num:
         ndyn.numplot(r,tmin,tmax,nsteps,t0=t0,lsig=lsig,discord=ndiscord,coh=coh,statusreport=statusreport,dephasing=dephasing)
         pyplot.subplot(1+lsplot+discplot,1+numplot+coplot,1+numplot)
         pyplot.title('Numerical, no RWA, '+state)
         pyplot.plot(ndyn.x,ndyn.flop)
         pyplot.ylim( (ymin, ymax) )
         pyplot.xlim( (tmin, tmax) )
         if dephasing:
             pyplot.axvline(x=t0,ls=':',color='k')
             pyplot.plot(ndyn.x,ndyn.deph)
         if lsig:
             pyplot.subplot(1+lsplot+discplot,1+numplot+coplot,3+numplot+coplot)
             pyplot.plot(ndyn.x,ndyn.lsig)
             pyplot.ylim( (ysmin, ysmax) ) 
             pyplot.axvline(x=t0,ls=':',color='k')
             pyplot.xlim( (tmin, tmax) )
#.........这里部分代码省略.........

def get_nbar(flop_directory,blue_file,red_file,fit_until=U.WithUnit(1000.0,'us'),show=False):
    print 'obtaining nbar from peak ratio of red and blue flop ...',
    #parameters and initial guesses for fit
    sideband = 1.0
    amax=2000.0
    f_Rabi_init = U.WithUnit(150.0,'kHz')
    nb_init = 0.1
    delta_init = U.WithUnit(1000.0,'Hz')
    fit_range_min=U.WithUnit(0.0,'us')
    fit_range_max=fit_until
    delta_fluc_init=U.WithUnit(100.0,'Hz')
    
    #actual script starts here
    class Parameter:
        def __init__(self, value):
                self.value = value
    
        def set(self, value):
                self.value = value
    
        def __call__(self):
                return self.value
            
    def fit(function, parameters, y, x = None):
        def f(params):
            i = 0
            for p in parameters:
                p.set(params[i])
                i += 1
            return y - function(x)
    
        if x is None: x = np.arange(y.shape[0])
        p = [param() for param in parameters]
        return optimize.leastsq(f, p)
    
    #get access to servers
    cxn = labrad.connect('192.168.169.197', password = 'lab')
    dv = cxn.data_vault
    
    #get trap frequency
    dv.cd(flop_directory)
    dv.cd(blue_file)
    dv.open(1)
    sideband_selection = dv.get_parameter('RabiFlopping.sideband_selection')
    sb = np.array(sideband_selection)
    trap_frequencies = ['TrapFrequencies.radial_frequency_1','TrapFrequencies.radial_frequency_2','TrapFrequencies.axial_frequency','TrapFrequencies.rf_drive_frequency']
    trap_frequency = dv.get_parameter(str(np.array(trap_frequencies)[sb.nonzero()][0]))            
    
    #SET PARAMETERS
    nb = Parameter(nb_init)
    f_Rabi = Parameter(f_Rabi_init['Hz'])
    delta = Parameter(delta_init['Hz'])
    delta_fluc=Parameter(delta_fluc_init['Hz'])
    #which to fit?
    fit_params = [nb,f_Rabi,delta,delta_fluc]
    
    # take Rabi flops
    data = dv.get().asarray
    blue_flop_y_axis = data[:,1]
    blue_flop_x_axis = data[:,0]*10**(-6)
    dv.cd(1)
    
    dv.cd(red_file)
    dv.open(1)
    data = dv.get().asarray
    red_flop_y_axis = data[:,1]
    red_flop_x_axis = data[:,0]*10**(-6)
    
    #fit Rabi Flops to theory
    blue_evo=tp.time_evolution(trap_frequency, sideband,nmax = amax)
    def blue_f(x):
        evolution = blue_evo.state_evolution_fluc(x,nb(),f_Rabi(),delta(),delta_fluc())
        return evolution
    
    red_evo=tp.time_evolution(trap_frequency, -sideband,nmax = amax)
    def red_f(x):
        evolution = red_evo.state_evolution_fluc(x,nb(),f_Rabi(),delta(),delta_fluc())
        return evolution
    
    #FIT BLUE
    
    fitting_region = np.where((blue_flop_x_axis >= fit_range_min['s'])&(blue_flop_x_axis <= fit_range_max['s']))
    fit(blue_f, fit_params, y = blue_flop_y_axis[fitting_region], x = blue_flop_x_axis[fitting_region])
    
    blue_nicer_resolution = np.linspace(0,blue_flop_x_axis.max(),1000)
    blue_flop_fit_y_axis = blue_evo.state_evolution_fluc(blue_nicer_resolution, nb(), f_Rabi(), delta(),delta_fluc())
    m=pylab.unravel_index(np.array(blue_flop_fit_y_axis).argmax(), np.array(blue_flop_fit_y_axis).shape)
    
    blue_max = np.array(blue_flop_fit_y_axis).max()
    
    fit_params = [nb,delta,delta_fluc]
    #FIT RED
    fitting_region = np.where((red_flop_x_axis >= fit_range_min['s'])&(red_flop_x_axis <= fit_range_max['s']))
    fit(red_f, fit_params, y = red_flop_y_axis[fitting_region], x = red_flop_x_axis[fitting_region])
    
    red_nicer_resolution = np.linspace(0,red_flop_x_axis.max(),1000)
    red_flop_fit_y_axis = red_evo.state_evolution_fluc(red_nicer_resolution, nb(), f_Rabi(), delta(),delta_fluc())
    
    red_max = red_flop_fit_y_axis[m]
    
#.........这里部分代码省略.........

#FIT BLUE

fitting_region = np.where((blue_flop_x_axis >= fit_range_min['s'])&(blue_flop_x_axis <= fit_range_max['s']))
print 'Fitting blue...'
p,success = fit(blue_f, fit_params, y = blue_flop_y_axis[fitting_region], x = blue_flop_x_axis[fitting_region])
print 'Fitting DONE.'

print "nbar = {}".format(nb())
print "Rabi Frequency = {} kHz".format(f_Rabi()*10**(-3))
print "The detuning is ({:.2f} +- {:.2f}) kHz".format(delta()*10**-3,np.abs(delta_fluc())*10**-3)

blue_nicer_resolution = np.linspace(0,blue_flop_x_axis.max(),1000)
blue_flop_fit_y_axis = blue_evo.state_evolution_fluc(blue_nicer_resolution, nb(), f_Rabi(), delta(),delta_fluc())

m=pylab.unravel_index(np.array(blue_flop_fit_y_axis).argmax(), np.array(blue_flop_fit_y_axis).shape)
print 'blue sideband highest peak value = {} at {}'.format(np.array(blue_flop_fit_y_axis).max(),10**6*blue_nicer_resolution[m])
blue_max = np.array(blue_flop_fit_y_axis).max()

blue_fit_nbar = nb()

fit_params = [nb,delta,delta_fluc]
#FIT RED
fitting_region = np.where((red_flop_x_axis >= fit_range_min['s'])&(red_flop_x_axis <= fit_range_max['s']))
print 'Fitting red...'
p,success = fit(red_f, fit_params, y = red_flop_y_axis[fitting_region], x = red_flop_x_axis[fitting_region])
print 'Fitting DONE.'

print "red nbar = {}".format(nb())
print "Rabi Frequency = {} kHz".format(f_Rabi()*10**(-3))
print "The detuning is ({:.2f} +- {:.2f}) kHz".format(delta()*10**-3,np.abs(delta_fluc())*10**-3)