def flow_rate_hist(sheets):
    ant_rates = []
    weights = []
    for sheet in sheets:
        ants, seconds, weight = flow_rate(sheet)
        ant_rate = seconds / ants
        #ant_rate = ants / seconds
        ant_rates.append(ant_rate)
        weights.append(float(weight))
        #weights.append(seconds)

    weights = pylab.array(weights)
    weights /= sum(weights)

    #print "ants per second"
    print "seconds per ant"
    mu = pylab.mean(ant_rates)
    print "mean", pylab.mean(ant_rates)
    wmean = pylab.average(ant_rates, weights=weights)
    print "weighted mean", wmean
    print "median", pylab.median(ant_rates)
    print "std", pylab.std(ant_rates, ddof=1)
    ant_rates = pylab.array(ant_rates)
    werror = (ant_rates - mu) * weights
    print "weighted std", ((sum(werror ** 2))) ** 0.5
    print "weighted std 2", (pylab.average((ant_rates - mu)**2, weights=weights)) ** 0.5
    pylab.figure()
    pylab.hist(ant_rates)
    pylab.savefig('ant_flow_rates.pdf', format='pdf')
    pylab.close()

def makeplot(filename):
    T0 = 2452525.374416
    P = 0.154525
    
    X = pl.load(filename)
    x = X[:,0]
    y = X[:,1]
    print x[0] # check for HJD faults
    
    #orbital phase
    p = (x-T0)/P
    
    pl.figure(figsize=(6,4))
    pl.subplots_adjust(hspace=0.47,left=0.16)
    
    pl.subplot(211)
    pl.scatter(p,y,marker='o',s=0.1,color='k')
    pl.ylim(-0.06,0.06)
    pl.xlim(pl.average(p)-1.25,pl.average(p)+1.25)
    pl.ylabel('Intensity')
    pl.xlabel('Orbital Phase')
    
    pl.subplot(212)
    f,a = ast.signal.dft(x,y,0,4000,1)
    pl.plot(f,a,'k')
    pl.ylabel('Amplitude')
    pl.xlabel('Frequency (c/d)')
    #pl.ylim(yl[0],yl[1])
    
    #pl.vlines(3636,0.002,0.0025,color='k',linestyle='solid')
    #pl.vlines(829,0.002,0.0025,color='k',linestyle='solid')
    #pl.text(3500,0.00255,'DNO',fontsize=11)
    #pl.text(700,0.00255,'lpDNO',fontsize=11)
    pl.ylim(0.0,0.004)
    pl.savefig('%spng'%filename[:-3])

def AT(s_ij=0, s_ik=0, s_jk=0, S=0):
	'''
	Calculates: SEFD = (2k/S) * (s_ij*s_ik)/(s_jk)
	'''
	kb = 1.38e3 # Boltzmann's Constant in Jy m^2 K^-1
	s_ij = pl.average(s_ij)
	s_ik = pl.average(s_ik)
	s_jk = pl.average(s_jk)
	return (2*kb/S)*(s_ij*s_ik)/(s_jk-s_ij*s_ik)

def visualize ():
    sample_rate, snd = load_sample(".\\hh-closed\\dh9.WAV")
    print snd.dtype
    data = normalize(snd)
    print data.shape
    n = data.shape[0]
    length = float(n)
    print length / sample_rate, "s"
    timeArray = arange(0, length, 1)
    timeArray = timeArray / sample_rate
    timeArray = timeArray * 1000  #scale to milliseconds
    ion()
    if False:
        plot(timeArray, data, color='k')
        ylabel('Amplitude')
        xlabel('Time (ms)')
        raw_input("press enter")
        exit()
    p = fft(data) # take the fourier transform
    nUniquePts = ceil((n+1)/2.0)
    print nUniquePts
    p = p[0:nUniquePts]
    p = abs(p)
    p = p / float(n) # scale by the number of points so that
                 # the magnitude does not depend on the length
                 # of the signal or on its sampling frequency
    p = p**2  # square it to get the power

    # multiply by two (see technical document for details)
    # odd nfft excludes Nyquist point
    if n % 2 > 0: # we've got odd number of points fft
        p[1:len(p)] = p[1:len(p)] * 2
    else:
        p[1:len(p) -1] = p[1:len(p) - 1] * 2 # we've got even number of points fft

    print p
    freqArray = arange(0, nUniquePts, 1.0) * (sample_rate / n);
    plot(freqArray/1000, 10*log10(p), color='k')
    xlabel('Frequency (kHz)')
    ylabel('Power (dB)')
    raw_input("press enter")

    m = average(freqArray, weights = p)
    v = average((freqArray - m)**2, weights= p)
    r = sqrt(mean(data**2))
    s = var(data**2)
    print "mean freq", m #TODO: IMPORTANT: this is currently the mean *power*, not the mean freq.  What we want is mean freq weighted by power
    print "var freq", v
    print "rms", r
    print "squared variance", s

def show_grey_channels(I):
    K = average(I, axis=2)
    for i in range(3):
        J = zeros_like(I)
        J[:, :, i] = K
        figure(i+10)
        imshow(J)

def hist_values(parameter, group, strategy, decay_type, label, y_limit=None):
    values = group[parameter]
    values = list(values)
    
    binsize = 0.05
    if 'zoom' in label:
        binsize = 0.01
    
    if y_limit == None:
        y_limit = max(values)
    cutoff = 1
    #weights = np.ones_like(values)/float(len(values))
    
    weights = group['num_lines'] / sum(group['num_lines'])
    weights = np.array(weights)
    
    mu = pylab.average(values, weights=weights)
    sigma2 = pylab.var(values)
    
    pylab.figure()
    pylab.hist(values, weights=weights, bins=np.arange(0, cutoff + binsize, binsize))
    title_items = []
    title_items.append('%s maximum likelihood values %s %s %s' % (parameter, strategy, decay_type, label))
    title_items.append('mean of estimates = %f' % mu)
    title_items.append('variance of estimates = %f' % sigma2)
    title_str = '\n'.join(title_items)
    #pylab.title(parameter + ' maximum likelihood values ' + str(strategy) + ' ' + str(outname))
    #pylab.title(title_str)
    print title_str
    pylab.xlabel('%s mle' % parameter, fontsize=20)
    pylab.ylabel('weighted proportion', fontsize=20)
    pylab.xlim((0, 1))
    pylab.ylim((0, y_limit))
    pylab.savefig('repair_ml_hist_%s_%s_%s_%s.pdf' % (parameter, strategy, decay_type, label), format='pdf')
    pylab.close()

 def movingaverage(x,L):
     ma = pl.zeros(len(x),dtype='Float64')
     # must take the lead-up zone into account (prob slow)
     for i in range(0,L):
         ma[i] = pl.average(x[0:i+1])
 
     for i in range(L,len(x)):
         ma[i] = ma[i-1] + 1.0/L*(x[i]-x[i-L])
         
     return ma

def img2ascii(filename, map_array=None):
    a = imread(filename)

    print "Converting ..."
    # useful only when reading .jpg files.
    # PIL is used for jpegs; converting PIL image to numpy array messes up. 
    # a = a[::-1, :] 

    # convert image to grayscale.
    if len(a.shape) > 2:
        a = 0.21 * a[:,:,0] + 0.71 * a[:,:,1] + 0.07 * a[:,:,2]
    a_r, a_c = a.shape[:2]
    a_max = float(a.max())

    blk_siz = 1 #size of block

    if map_array == None:
        # just linearly map gray level to characters.
        # giving lowest gray level to space character.
        out_file = open(filename + 'lin' + str(blk_siz) + '.txt', 'w')
        print "File %s opened" %out_file.name
        for i in range(0, a_r, blk_siz*2):
            for j in range(0, a_c, blk_siz):
                b = a[i:i+2*blk_siz, j:j+blk_siz]
                b_char = chr(32+int((1-average(b))*94))
                out_file.write(b_char)
            out_file.write("\n")
        out_file.close()
    else:
        # map based on visual density of characters.
        out_file = open(filename + 'arr' + str(blk_siz) + '.txt', 'w')
        print "File %s opened" %out_file.name
        for i in range(0, a_r, blk_siz*2):
            for j in range(0, a_c, blk_siz):
                b = a[i:i+2*blk_siz, j:j+blk_siz]
                b_mean = int(average(b)/a_max*(len(map_array)-1))
                b_char = chr(map_array[b_mean])
                out_file.write(b_char)
            out_file.write("\n")
        out_file.close()
        
    print "%s Converted! \nWritten to %s" %(filename, out_file.name)

def sigclip(im,nsig):
    # returns min and max values of image inside nsig sigmas
    temp = im.ravel()
    sd = pl.std(temp)
    m = pl.average(temp)
    gt = temp > m-nsig*sd
    lt = temp < m+nsig*sd
    temp = temp[gt*lt]
    mini = min(temp)
    maxi = max(temp)
    
    return mini,maxi

    def plot(self):
        """generate the plot formatting"""
        if self.data == None:
            print "Must load and parse data first"
            sys.exit()
            
        for k,v in self.data.iteritems():
            for type, data in v.iteritems():
                pylab.clf()
                height = int(self.height)
                width = int(self.width)
                pylab.figure()
                ax = pylab.gca()
                ax.set_xlabel('<--- Width = %s wells --->' % str(width))
                ax.set_ylabel('<--- Height = %s wells --->' % str(height))
                ax.set_yticks([0,height/10])
                ax.set_xticks([0,width/10])
                ax.set_yticklabels([0,height])
                ax.set_xticklabels([0,width])
                ax.autoscale_view()
                pylab.jet()
            #color = self.makeColorMap()
            #remove zeros for calculation of average
                flattened = []
                for i in data:
                    for j in i:
                        flattened.append(j)
                flattened = filter(lambda x: x != 0.0, flattened)
                Aver=pylab.average(flattened)
                name = type.replace(" ", "_")
                fave = ("%.2f") % Aver
                pylab.title(k.strip().split(" ")[-1] + " Heat Map (Average = "+fave+"%)")
                ticks = None
                vmax = None
                if type == "Region DR":
                    ticks = [0.0,0.2,0.4,0.6,0.8,1.0]
                    vmax = 1.0
                else:
                    ticks = [0.0,0.4,0.8,1.2,1.6,2.0]
                    vmax = 2.0
                    
                pylab.imshow(data, vmin=0, vmax=vmax, origin='lower')
                pylab.colorbar(format='%.2f %%',ticks=ticks)
                pylab.vmin = 0.0
                pylab.vmax = 2.0
            #pylab.colorbar()     

                if self.savePath is None:
                    save = "%s_heat_map_%s.png" % (name,k)
                else:
                    save = path.join(self.savePath,"%s_heat_map_%s.png" % (name,k))
                pylab.savefig(save)
                pylab.clf()

def process_window(sample_rate, data):
    # print "processing window"
    # print data.dtype
    # print data.shape
    n = data.shape[0]
    length = float(n)
    # print length / sample_rate, "s"
    p = fft(data) # take the fourier transform
    nUniquePts = ceil((n+1)/2.0)
    p = p[0:nUniquePts]
    p = abs(p)
    p = p / float(n) # scale by the number of points so that
                 # the magnitude does not depend on the length
                 # of the signal or on its sampling frequency
    p = p**2  # square it to get the power

    # multiply by two (see technical document for details)
    # odd nfft excludes Nyquist point
    if n % 2 > 0: # we've got odd number of points fft
        p[1:len(p)] = p[1:len(p)] * 2
    else:
        p[1:len(p) -1] = p[1:len(p) - 1] * 2 # we've got even number of points fft
    freqArray = arange(0, nUniquePts, 1.0) * (sample_rate / n);

    if sum(p) == 0:
        raise Silence
    m = average(freqArray, weights = p)
    v = sqrt(average((freqArray - m)**2, weights= p))
    r = sqrt(mean(data**2))
    s = var(data**2)
    print "mean freq", m #TODO: IMPORTANT: this is currently the mean *power*, not the mean freq.  What we want is mean freq weighted by power
    # print freqArray
    # print (freqArray - m)
    # print p
    print "var freq", v
    print "rms", r
    print "squared variance", s
    return [m, v, r, s]

 def writeMetricsFile(self, filename):
     '''Writes the lib_cafie.txt file'''
     if self.data == None:
         print "Must load and parse data first"
         sys.exit()
     cafie_out = open(filename,'w')
     for k,v in self.data.iteritems():
         for type, data in v.iteritems():
             flattened = []
             for i in data:
                 for j in i:
                     flattened.append(j)
             flattened = filter(lambda x: x != 0.0, flattened)
             Aver=pylab.average(flattened)
             name = type.replace(" ", "_")
             if 'LIB' in k:
                 if len(flattened)==0:
                     cafie_out.write('%s = %s\n' % (name,0.0))
                     Aver = 0
                 else:
                     cafie_out.write('%s = %s\n' % (name,Aver))
                     Aver=pylab.average(flattened)
     cafie_out.close()

 def calcTDData(self,tdDatas):
     #tdDatas is a a 3d array of measurements, along with their uncertainties
     #meantdData is the weighted sum of the different measurements
     #meantdData,sumofweights=py.average(tdDatas[:,:,1:3],axis=0,weights=1.0/tdDatas[:,:,3:]**2,returned=True)
     meantdData=py.average(tdDatas[:,:,1:3],axis=0)
     #use error propagation formula
     noise=py.sqrt(py.mean(self.getAllPrecNoise()[0]**2))
     if tdDatas.shape[0]==1:
         rep = py.zeros((len(tdDatas[0,:,0]),2))
     else:
         rep = py.std(tdDatas[:,:,1:3],axis=0, ddof=1)/py.sqrt(self.numberOfDataSets)
     unc = py.sqrt(rep**2+noise**2)
     #unc=py.sqrt(1.0/sumofweights)
     #time axis are all equal
     return py.column_stack((tdDatas[0][:,0],meantdData,unc))       

def amptime(uv, baseline="0_1", pol="xx", applycal=False):
	'''
	Plots Amp vs Time for a single baseline. 
	'''
	fig = pl.figure()
	ax = fig.add_subplot(111)
	aipy.scripting.uv_selector(uv, baseline, pol)
	for preamble, data, flags in uv.all(raw=True):
		uvw, t, (i, j) = preamble
		ax.plot(t, pl.average(pl.absolute(data)), 'ks', mec='None', alpha=0.2, ms=5)
	hfmt = dates.DateFormatter('%m/%d %H:%M')
	ax.xaxis.set_major_locator(dates.HourLocator())
	ax.xaxis.set_major_formatter(hfmt)
	ax.set_ylim(bottom = 0)
	pl.xticks(rotation='vertical')	

def calculateScores(arr, HEIGHT, WIDTH):
    rowlen,collen = arr.shape
    scores = pylab.zeros(((rowlen/INCREMENT),(collen/INCREMENT)))
    score = []
    for row in range(rowlen/INCREMENT):
        for column in range(collen/INCREMENT):            
            keypassed,size = getAreaScore(row,column,arr, HEIGHT, WIDTH)
            scores[row,column] = round(float(keypassed)/float(size)*100,2)
            if keypassed > 2:
                score.append(round(float(keypassed)/float(size)*100,2))         
            
            #scores[0,0] = 0
            #scores[HEIGHT/INCREMENT -1,WIDTH/INCREMENT -1] = 100
    print scores
    
    flattened = []
    for i in score:
        flattened.append(i)
    flattened = filter(lambda x: x != 0.0, flattened)    
    average=pylab.average(flattened)
    
    return score, scores, average

def plot_effect_num_cashiers_on_cust_wait_time(customersPerMinute = 10, numCashiersToTestUpTo = 12):
    assert customersPerMinute > 0
    assert numCashiersToTestUpTo > 0
    assert type(customersPerMinute) == type(numCashiersToTestUpTo) == int
    store = Store(customersPerMinute)
    worstCase = []
    averageCase = []
    rangeOfNumCashiers = range(1, numCashiersToTestUpTo + 1)
    for i in rangeOfNumCashiers:
        store.run_simulation(i)
        timeOnLineData = [x.timeOnLine / 60. for x in store.completedCustomers]
        averageCase.append(pylab.average(timeOnLineData))
        worstCase.append(max(timeOnLineData))
        store.reset_store()
    
    pylab.plot(rangeOfNumCashiers, worstCase, label='Longest Time on Line') 
    pylab.plot(rangeOfNumCashiers, averageCase, label = 'Average Time on Line')
    
    pylab.title('Effect of Adding Additional Cashiers \n on Customer Wait Time')
    pylab.xlabel('Number of Cashiers')
    pylab.ylabel('Customer Wait Time in Minutes \n (if store receives {} customers per minute)'.format(store.customersPerMinute))
    pylab.legend()
    pylab.xticks(rangeOfNumCashiers)  
    pylab.show()

def Froudenumber(flmlname):
  print "\n********** Calculating the Froude number\n"
  # warn user about assumptions
  print "Froude number calculations makes three assumptions: \n i) domain height = 0.1m \n ii) mid point domain is at x = 0.4 \n iii) initial temperature difference is 1.0 degC"
  domainheight = 0.1
  domainmid = 0.4
  rho_zero, T_zero, alpha, g = le_tools.Getconstantsfromflml(flmlname)
  gprime = rho_zero*alpha*g*1.0 # this has assumed the initial temperature difference is 1.0 degC

  # get list of vtus
  filelist = le_tools.GetFiles('./')
  logs = ['diagnostics/logs/time.log','diagnostics/logs/X_ns.log','diagnostics/logs/X_fs.log']
  try:
    # if have extracted information already just use that
    os.stat('diagnostics/logs/time.log')
    os.stat('diagnostics/logs/X_ns.log')
    os.stat('diagnostics/logs/X_fs.log')
    time = le_tools.ReadLog('diagnostics/logs/time.log')
    X_ns = [x-domainmid for x in le_tools.ReadLog('diagnostics/logs/X_ns.log')]
    X_fs = [domainmid-x for x in le_tools.ReadLog('diagnostics/logs/X_fs.log')]
  except OSError:
    # otherwise get X_ns and X_fs and t from vtus
    time, X_ns, X_fs = le_tools.GetXandt(filelist)
    f_time = open('./diagnostics/logs/time.log','w')
    for t in time: f_time.write(str(t)+'\n')
    f_time.close()
    f_X_ns = open('./diagnostics/logs/X_ns.log','w')
    for X in X_ns: f_X_ns.write(str(X)+'\n')
    f_X_ns.close()
    f_X_fs = open('./diagnostics/logs/X_fs.log','w')
    for X in X_fs: f_X_fs.write(str(X)+'\n')
    f_X_fs.close()

    # shift so bot X_ns and X_fs are 
    # distance of front from 
    #initial position (mid point of domain)
    X_ns = [x-domainmid for x in X_ns]
    X_fs = [domainmid-x for x in X_fs]

  # Calculate U_ns and U_fs from X_ns, X_fs and t
  U_ns = le_tools.GetU(time, X_ns)
  U_fs = le_tools.GetU(time, X_fs)
  U_average = [[],[]]

  # If possible average 
  # (if fronts have not travelled far enough then will not average)
  start_val, end_val, average_flag_ns = le_tools.GetAverageRange(X_ns, 0.2, domainheight)
  if average_flag_ns == True: U_average[0].append(pylab.average(U_ns[start_val:end_val]))
  
  start_val, end_val, average_flag_fs = le_tools.GetAverageRange(X_fs, 0.25, domainheight)
  if average_flag_fs == True: U_average[1].append(pylab.average(U_fs[start_val:end_val]))
  
  # plot
  fs = 18
  pylab.figure(num=1, figsize = (16.5, 11.5))
  pylab.suptitle('Front speed', fontsize = fs)

  pylab.subplot(221)
  pylab.plot(time,X_ns, color = 'k')
  pylab.axis([0,45,0,0.4])
  pylab.grid('on')
  pylab.xlabel('$t$ (s)', fontsize = fs)
  pylab.ylabel('$X$ (m)', fontsize = fs)
  pylab.title('no-slip', fontsize = fs)
    
  pylab.subplot(222)
  pylab.plot([x/domainheight for x in X_ns],[U/math.sqrt(gprime*domainheight) for U in U_ns], color = 'k')
  pylab.axis([0,4,0,0.6])
  pylab.grid('on')
  pylab.axhline(0.406, color = 'k')
  pylab.axhline(0.432, color = 'k')
  pylab.text(3.95,0.396,'Hartel 2000',bbox=dict(facecolor='white', edgecolor='black'), va = 'top', ha = 'right')
  pylab.text(3.95,0.442,'Simpson 1979',bbox=dict(facecolor='white', edgecolor='black'), ha = 'right')
  pylab.xlabel('$X/H$', fontsize = fs)
  pylab.ylabel('$Fr$', fontsize = fs)
  pylab.title('no-slip', fontsize = fs)
  if average_flag_ns == True:
    pylab.axvline(2.0, color = 'k')
    pylab.axvline(3.0, color = 'k')
    pylab.text(0.05, 0.01, 'Average Fr = '+'{0:.2f}'.format(U_average[0][0]/math.sqrt(gprime*domainheight))+'\nvertical lines indicate the range \nover which the average is taken', bbox=dict(facecolor='white', edgecolor='black'))
  
  pylab.subplot(223)
  pylab.plot(time,X_fs, color = 'k')
  pylab.axis([0,45,0,0.4])
  pylab.grid('on')
  pylab.xlabel('$t$ (s)', fontsize = fs)
  pylab.ylabel('$X$ (m)', fontsize = fs)
  pylab.title('free-slip', fontsize = fs)
    
  pylab.subplot(224)
  pylab.plot([x/domainheight for x in X_fs],[U/math.sqrt(gprime*domainheight) for U in U_fs], color = 'k')
  pylab.axis([0,4,0,0.6])
  pylab.grid('on')
  pylab.axhline(0.477, color = 'k')
  pylab.text(3.95,0.467,'Hartel 2000', va = 'top',bbox=dict(facecolor='white', edgecolor='black'), ha = 'right')
  pylab.xlabel('$X/H$', fontsize = fs)
  pylab.ylabel('$Fr$', fontsize = fs)
  pylab.title('free-slip', fontsize = fs)
  if average_flag_fs == True:
    pylab.text(0.05, 0.01, 'Average Fr  = '+'{0:.2f}'.format(U_average[1][0]/math.sqrt(gprime*domainheight))+'\nvertical lines indicate the range \nover which the average is taken', bbox=dict(facecolor='white', edgecolor='black'))
#.........这里部分代码省略.........

#print '\nDone...\n'



# bin spectra together in 0.1 phase bins
im3 = []
klist = []
k = 0
temp = pl.zeros(len(pf.getdata(ff[0])[xx]),dtype=float)

for i in range(0,100):
    for j in range(len(ff)):
        if phase[j] >= i*0.01 and phase[j] < (i+1)*0.01:
             temp += pf.getdata(ff[i])[xx]-pl.median(pf.getdata(ff[i])[xx])#pf.getdata(ff[argsort[j]])
             k+=1
    ave = pl.average(temp)
    im3.append(temp)
    temp = pl.zeros(len(pf.getdata(ff[0])[xx]),dtype=float)
    klist.append(k)
    k = 0

#print sum(klist)
#print len(klist)




pl.figure()
pl.gray()
pl.imshow(im3, interpolation='nearest', aspect='auto',cmap=pl.cm.gray_r,extent=(6500,6625,1,0))
#ax = pl.axes()

sa = []

for f in files:
    print 'Reading %s ' % f
    name, eph = string.split(f)
    T0 = ephemeris[eph]
    P = 0.154525
    
    X = pl.load(name)
    x = (X[:,2] - T0)/P
    xx.append(x - int(x[0]))
    aa.append(X[:,0])
    
    # let phase at 0.8 -> 1.0
    tpp = X[:,1]
    tpp -= pl.average(tpp[0])
    #tpp += 1.0
    
    pp.append(tpp)
    sa.append(X[:,3])
    sp.append(X[:,4])
    

# now sort observations in terms of orbital phase
xx = pl.array([i for i in pl.flatten(xx)])
pp = pl.array([i for i in pl.flatten(pp)])
aa = pl.array([i for i in pl.flatten(aa)])
sa = pl.array([i for i in pl.flatten(sa)])
sp = pl.array([i for i in pl.flatten(sp)])

arg = xx.argsort()

 
 
 # write average subtracted spectrum to new fits file
 #pf.writeto('avesub%s'%i,data=data,header=head)
 
 start = head['CRVAL1']
 step = head['CDELT1']
 length = head['NAXIS1']
 x = start + pl.arange(0,length)*step
 
 # hydrogen alpha
 dl = v/c*6563.0
 w1 = x > 6563 - dl
 w2 = x < 6563 + dl
 
 imHa.append(data[w1*w2]-pl.average(data[w1*w2]))
 #imHa.append((data[w1*w2]))
 #imHa.append((data[w1*w2]-pl.average(data[w1*w2])-(ave[w1*w2]-pl.average(ave[w1*w2]))))
 
 dl = v/c*4860.0
 w1 = x > 4860 - dl
 w2 = x < 4860 + dl
 #data = pf.getdata(i)
 imHb.append(data[w1*w2]-pl.average(data[w1*w2]))
 #imHb.append((data[w1*w2]-pl.average(data[w1*w2])-(ave[w1*w2]-pl.average(ave[w1*w2]))))
 #imHb.append(data[w1*w2])
 
 dl = v/c*4686
 w1 = x > 4686 - dl
 w2 = x < 4686 + dl
 #data = pf.getdata(i)