def tHigh(data,z,threshold,prior_array,phigh,pint,no_val,vert_val,fht_val,offset):
#look from the previous higher layer and start down
#then check to see if there are values below
#with  backscatter higher than phigh more than no_val times
#which may arise from a cloud, saturation or attenuation
#adjust down one half layer to reflect the top of the aerosol layer
# then if the lower layer is within offset of the upper layer, 
# set the lower layer to missing

    import numpy as np
    value=[]
    for i in range(len(data)):
        v=[]
        try:
            w=np.where(prior_array[i] == z)[0]
            zeta=list(reversed(z[:w[0]]))
            beta=list(reversed(data[i][:w[0]]))
            try:
                v=np.where((threshold + pint/2.) <= beta)
                elv = zeta[v[0][0]]-vert_val
                if elv < fht_val:
                   #ignore really low ones
                   value.append(-9999)
                else:
                	diff = prior_array[i] - elv
                	if (diff <= offset):
                        	# split the difference for both
                        	mn_val = prior_array[i] - diff/2.
		      		prior_array[i] =  mn_val
                        	value.append(mn_val)
			else:
                		value.append(elv)
            except:
                value.append(-9999)
        except:
            value.append(-9999)
    #print 'other value:', value
    return np.array(value),prior_array
def count_the_highest_gradient(bs, unit_id, dt,fht,lht,phigh,plow,pint,grad_th):
    import numpy as np
    import matplotlib.pyplot as plt
#set thresholds
    # want grad to be greater than grad_th
    # will drop to half that if don't find one

    no_bin = np.rint((phigh - plow)/pint + 1)
    a = np.linspace(plow, phigh, no_bin)
    #print a
    all_counts = []
    p = []
    total = 0
    for val in a:
        counted = []
#        for v in bs:
        for v in bs[:,fht:lht]:
# handle values below the plow threshold
            if val == plow:
                 w = np.where(v <= val)
#handle values greater than or equal to the upper threshold
            elif val == phigh:
                 w = np.where(v > val )
            else:
                 # print 'Doing everything else'
                 w = np.where((v >= val) & (v < (val + pint)))
            counted.append(len(w[0]))
        #print 'sum of counted for ' + str(val) + ' is:', np.sum(counted)
        all_counts.append(np.sum(counted))
    total = float(np.sum(all_counts))
    # p.append((np.sum(counted)/total)*100)
    p = [(val / total) * 100 for val in all_counts]
    #total_prime = float(np.sum(all_counts[1:-1]))
    #p_prime = [(val / total_prime) * 100 for val in all_counts[1:-1]]
    #compute simple gradient between two adjacent values
    grad = np.diff(p[1:-1])
    #print "a",a
    #print "p",p
    #print "grad", grad
    #print "grad_th", grad_th
    values = []
    for i in range(1, len(grad) -1):
        #requiring that a gradient must be bigger than grad_th
        #and it is a local max in gradient
		     #if (grad[i] >= grad_th and grad[i+1] < 0):
		     if (grad[i] >= grad_th ):
                         values.append(a[i+2])
    #try again with lower threshold if didn't get any good gradients
    if not values:
         grad_th = grad_th / 2. 
	 for i in range(1, len(grad)-1):
		     #if (grad[i] >= grad_th and grad[i+1] < 0):
		     if (grad[i] >= grad_th ):
                         values.append(a[i+2])
#only returning now the first one
    if not values:
        values.append(-9999)
        print "no first gradient. need to add another lower threshold"
    else:
        print 'first gradient value set',values
    return p, values[0]
def count_the_percent_gradient(bs, unit_id, dt,fht,lht,phigh,plow,pint,grad_th):
    import numpy as np
    import matplotlib.pyplot as plt

    no_bin = np.rint((phigh - plow)/pint + 1)
    a = np.linspace(plow, phigh, no_bin)
    #print a
    all_counts = []
    p = []
    total = 0
    for val in a:
        counted = []

        for i in range(len(bs)):

# handle values below the tl threshold
            if val == plow:
                 w = np.where(bs[i][fht:lht[i]] <= val)
#handle values greater than or equal to the upper threshold
            elif val == phigh:
                 w = np.where(bs[i][fht:lht[i]] >= val )
            else:
                 # print 'Doing everything else'
                 w = np.where((bs[i][fht:lht[i]] >= val) & (bs[i][fht:lht[i]] < (val + pint)))
            counted.append(len(w[0]))
        all_counts.append(np.sum(counted))
    total = float(np.sum(all_counts))
    p = [(val / total) * 100 for val in all_counts]
    total_prime = float(np.sum(all_counts[1:-1]))
    p_prime = [(val / total_prime) * 100 for val in all_counts[1:-1]]
    values = []
    grad = np.diff(p[1:-1])
    #print "a",a
    #print "p",p
    #print "grad", grad
    # make threshold more stringent after first one
    values = []
    for i in range(1, len(grad) -1):
        #requiring that a gradient must be bigger than grad_th
        #and it is a local max in gradient
        if (grad[i] >= grad_th and grad[i+1] <0 ):
             values.append(a[i+2])
    print 'gradient values', values
    return p, grad, values
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
def tLow(data,z,threshold,phigh,pint,no_val,vert_val,fht_val):
# looks from the bottom up for the first layer with
#backscatter less than the threshold
#then check to see if there are values below 
#with  backscatter higher than phigh more than no_val times
#which may arise from a cloud, saturation or attenuation
#finally adjust down one layer to reflect the top of the aerosol layer
    import numpy as np
    value=[]
    for i in range(len(data)):
       v=[]
       v = np.where(data[i] < (threshold+pint/2.))
       if len(v[0])>0:
           vv = np.where(data[i,:v[0][0]]>phigh)
           #print 'i,v,vv',i,v[0][0],len(vv[0])
           if len(vv[0])<no_val:
               elv = z[v[0][0]]-vert_val
               if elv < fht_val:
                  elv = -9999
               value.append(elv)
           else:
               value.append(-9999)
       else:
           value.append(-9999)
    return np.array(value)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
def timemean(dat, time, dt, unit_id, verbose=True):
    import numpy as np
    import scipy.stats as stats
    import matplotlib.mlab as mlab

    '''
    move data into time-oriented bins, returning the binned average of multi-dimensional data


    Parameters
    ----------

    dat: numpy array
        the data to be averaged into bins, first dimension should correspond to time length

        etc...
    '''
    if dt == 0:
        'do nothing, the function was just invariably wrapped into something'
        return dat, time
    begin = np.min(time)
    end = np.max(time)
    length = int(np.floor((end - begin) / dt) + 1)
    outT = np.zeros(length)
    verts = dat.shape[1:]
    outshape = [length] + list(dat.shape[1:])
#    if verbose:
    #print 'outshape',outshape
    outD = np.zeros(outshape)
    binlow = begin
    i = -1
    while True:
        binhigh = binlow + dt
        i += 1
        outT[i] = binlow + dt / 2
        '... do the processing'
        q = dat[(time < binhigh) & (time >= binlow)]
        #print 'q',q
        if len(q) > 0:
            res = np.nanmedian(q, axis=0)
	    iqrt = stats.iqr(q,axis=0,nan_policy='propagate')
	    #print 'iqrt',iqrt
	    #print 'idx',idx
	    if unit_id == 'LOGCS':
                idx = mlab.find(iqrt>2.0)
            	res[idx] = -8.1
            else:
                idx = mlab.find(iqrt>1.0)
            	res[idx] = -7.7
        else:
            res = np.nan
        outD[i] = res
        binlow = binhigh
        if binhigh > end:
            break
    return outD, outT
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
def timebin(dat, time, dt, unit_id, verbose=False):
    '''
    A replacemaent for the name timebin, to preserve intercompatibility.
    '''
    return timemean(dat, time, dt, unit_id, verbose)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
def mean2d(dat, dim, binsize):
    '''
    compute a chunk/bin style mean along the first axis, with dimension key of the same size

        the first index should be time, and the second height, we will average by taking chunks in time
        and averaging

        THIS MEANS THAT THERE IS NO TRANSPOSE OPERATIONS NEEDED FOR COMPUTATION
    '''
    import numpy as np
    datO = np.zeros((int(dat.shape[0] / binsize), dat.shape[1]))  # initialize
    chunk = dat[0:binsize]
    i = 0  # index
    #print 'binsize',binsize
    while True:
        try:
            datO[i] = np.median(chunk, axis=0)
            i += 1
            # take a chunk of 'profiles' in 'time' ( ||| ||| ||| = 3 chunks)
            chunk = dat[i * binsize:(i + 1) * binsize]
        except:
            break
    # we could convolve as well, but that is just not as nice! (and not much faster either)
    return datO, mean1d(dim, binsize)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
def mean1d(dat, binsize):
    """
    an equivalent method for single dimension averaging (such as getting equivalent times)
    """
    import numpy as np
    dat=np.array(dat)
    out = np.zeros(int(dat.shape[0] / binsize))  # one dimensional only!!
    chunk = dat[0:binsize]
    i = 0  # index
    while True:
        try:
            out[i] = np.median(chunk, axis=0)  # axis does not have to be specified
            i += 1
            # take a chunk of 'profiles' in 'time' ( ||| ||| ||| = 3 chunks)
            chunk = dat[i * binsize:(i + 1) * binsize]
        except:
            break
    return out
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
def runmean(field, nbins=30, ax=0, nanmean=True):
    import numpy as np
    # from scipy.stats import nanmean as scipynanmean
    '''
    Compute a running mean along the axis specified, evaluating over the
    distance 'nbins' 2D datasets only.

    '''
    if ax == 1:
        field = field.T
    i = 0
    ln = len(field)
    newdata = np.zeros(field.shape)
    if nanmean:
        try:
            meanfunc = np.nanmedian
        except:
            from scipy.stats import nanmedian as scipynanmedian
            meanfunc = scipynanmedian
    else:
        meanfunc = np.median
    while i < len(field):
        'I am assuming that int/int = int'
        i0 = i - nbins / 2
        i1 = i + nbins / 2
        if i < nbins / 2:
            i0 = 0
        if i > ln - nbins / 2:
            i1 = ln

        clump = field[i0:i1]
        newdata[i] = meanfunc(clump, axis=0)
        i += 1
    field = newdata
    del newdata
    if ax == 1:
        field = field.T
    return field
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
def runmean2d(field,dim1bin,dim2bin):
    '''
    Compute a 2-d running mean on field field with lenghts (dim1bin) bins
    in the first dimension and dim2bin in the second
    '''
    import numpy as np

    field=np.array(field)
    field = runmean(field.T,dim2bin)
    field = runmean(field.T,dim1bin)

    return field
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
def _ComputeFieldMeans(data, binsize, unit_id, inTime=False, continuous=False, vertbin=2,power=False):
    '''
    Computed the background data field given the variables from a given data variable
    (slice output is taken as the input)
    binsize has different meanings depending on the routine used. if inTime = True then it is referring to time interval
    600 would be 10 min
    but if using mean2d then it is referring to the number of array elements
    inTime=False means do the filtering using mean2d rather than inTime (which smooth vertically only)
    continuous=False means skip over this invalid option
    vertbin = 2 means smoothing over 3 layers (30 m for Vaisala)
    power=False means filtering done in log of backscatter array not on raw backscatter

    '''
    import numpy as np
    if power:
        data['BS'] = 10 ** np.array(data['BS'])
    if continuous:
        'continuous is not computed WRT time ever, so it is irrelevant'
        # FIXME - this should be redone as a c
        field, times, z = runmean2d(data['BS'], binsize, vertbin), data['DATTIM'], data['HEIGHT']
        'compute means within binned values...'
    else:
        if vertbin > 1:
            field, z = runmean(np.array(data['BS']), vertbin, ax=1), np.array(data['HEIGHT'])
            '10*vertbin vertical running mean is given...'
        else:
            field, z = np.array(data['BS']), np.array(data['HEIGHT'])


        if inTime:
            field, times = timebin(field, data['DATTIM'], binsize, unit_id,verbose=False)
        else:
            'mean2d is meant to be much faster than timebin'
            field, times = mean2d(field, data['DATTIM'], binsize)
    del data
    return field, times, z

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
def zero2hero(threshold_array):
    import numpy as np
    t_a = threshold_array

    #    ok='not ok'
    #    #get rid of 0's
    for i in range(len(t_a) - 1):
        if t_a[0] == 0:
            t_a[i] = 100
        if t_a[-1] == 0:
            t_a[-1] = 100
        if isinstance(t_a[i], int) != True:
            t_a[i] = 100
        if isinstance(t_a[i + 1], int) != True:
            t_a[i + 1] = 100
        if t_a[i] == 0 and t_a[i + 1] != 0:
            t_a[i] = (t_a[i - 1] + t_a[i + 1]) / 2.0
        if t_a[i] == 0 and t_a[i + 1] == 0:
            ok = 'nope'
            num = 0
            while ok == 'nope':
                if i == len(t_a) - 1:
                    t_a[i] = t_a[i - 1]
                    ok = 'ok'
                if t_a[i + num] == 0 and t_a[i + num] != t_a[-1]:
                    num += 1
                if t_a[i + num] > 0:
                    t_a[i] = (t_a[i - 1] + t_a[i + num]) / 2.0
                    ok = 'ok'
                else:
                    t_a[i] = t_a[i - 1]
                    ok = 'ok'
    # print 't_a=',t_a
    return np.array(t_a)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
def what_do_i_got(time_array):
    time_diff=[]
    for i in range(len(time_array)-1):
        #print 'the time diff is: %0.2f'%(time_array[i+1]-time_array[i])
        time_diff.append(time_array[i+1]-time_array[i])
    return time_diff
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
# Determine if the values sky rocket or there's rain or virga
def it_came_from_above(bs,unit_id):
    flagged=[]
    for i in range(len(bs)):
        try:

            if unit_id=='LOGCS':
                no_flag = 38
                f = (filter(lambda x: x > -5.25, bs[i][140:]))
                if f:
                    flagged.append(f)
            else:
                no_flag = 19
                f=(filter(lambda x: x > -5.25, bs[i][70:]))
                if f:
                    flagged.append(f)
        except:
            print 'I dunno why this doesnt work'
    if len(flagged)>no_flag:
            print '^^r_flagged^^'
            return True
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
# If there is fog, tease out an array of indicies for the flagged values
def oh_that_was_Foghat(bs,unit_id):
    flagged = []
    for i in range(len(bs)):
        try:
            f = (filter(lambda x: x > -5.25, bs[i][:40]))
            if f:
                flagged.append(f)
        except:
            print 'I dunno why this doesnt work'

    if len(flagged) > 19:
        print '**f_flagged**'
        return True

def attenuation(smoothed_data,phigh,no_val):
    import numpy as np
    attenuated = []
    #looping over each time
    for i in range(len(smoothed_data)):
        # looking for saturated below hta
        at = (filter(lambda x: x > phigh, smoothed_data[i]))
        if len(at) > no_val:
            attenuated.append(i)
    return attenuated


# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
def git_api(start, end, stn, vars, ceilo=False):
    '''
    start - string format as YYYYmmddHHMM
    end - string format as YYYYmmddHHMM
    stn - can be an array or list, but needs to be a string
        ex: stn = 'naa' OR stn = 'naa,mtmet'
    vars - string of api accepted variables
        ex: vars = 'ozone_concentration' OR vars = 'ozone_concentration,air_temp'
    '''
    import json
    import urllib2
    if stn == 'LOGCS':
        URL = 'https://asn.synopticdata.com/api/v1/series?senabbr=' + stn + '&variables='+vars+'&token=e69f2cf0265d499cbae3e50dc423dd04&Begin=' + start + '&End=' + end + '&ImageDimensions=100000,5'
    else:
        URL = 'https://asn.synopticdata.com/api/v1/series?senabbr=' + stn + '&variables='+vars+'&token=e69f2cf0265d499cbae3e50dc423dd04&Begin=' + start + '&End=' + end + '&ImageDimensions=100000,5'
    print URL
    # Note the token belongs to Luke!!!!

    # Open URL and read the content
    f = urllib2.urlopen(URL,timeout=35)
    data = f.read()
    # convert that json string into some python readable format
    data = json.loads(data)
    return data, URL
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

def eradicate(array):
    import numpy as np
    temp1=[]
    for i in range(4,len(array)-1):
        if array[i]== -9999 and array[i-1]>0 and -9999 in array[i-5:i]:
           for j in range(i-5,i):
                #print j
                if array[j] < 0 and array[j+1] >= 0:
                    array[j+1]= -9999
    # np.array(temp1)
    # print 'eradiacated array:',temp1
    # for val in temp1:
    #     array[val]=-9999

    return array

def all_the_plots(mwsiteid,date_array,height,BS,threshold_array,tl,th,p_wid,dt,unit_id,lht2,gradient,multi,percent_first,percents,d_2=False,bs_smoothed=False,realtime=False):
    import numpy as np
    import time, datetime
    import matplotlib
    import matplotlib.colors as colors
    matplotlib.use('agg')
    import matplotlib.pyplot as plt
    from matplotlib import rc
    import matplotlib.dates as md
    # rc('text', usetex=True)

    print date_array[-1]
    print time.strftime('%Y-%m-%d %H%M',time.localtime(date_array[0]))
    print time.strftime('%Y-%m-%d %H%M',time.localtime(date_array[-1]))

    color=['saddlebrown','k','r','orange','green','purple','limegreen','rosybrown','g','b','orange','y','c']

    stntitlenames = {
       'MTMET': 'MTMET (1523 m ASL)',
       'USDR1': 'USDR1 (1290 m ASL)',
    }

    myd=[]
    myd_2=[]
    t_a = []
    t = []
    if not threshold_array.any:
        print "no thresholds to plot"
    else:
        #print "threshold_array",threshold_array
        for r in range(len(threshold_array)):
            v = threshold_array[r]
            #print "len(v)",v,len(v)
            t = [None if v[i] <= 60 else v[i] for i in range(len(v))]
            t_a.append(t)
    for i in range(len(date_array)):
        myd.append(md.date2num(datetime.datetime.fromtimestamp(date_array[i])))
    try:
       if d_2:
           myd_2=[]
           for i in range(len(d_2)):
                myd_2.append(md.date2num(datetime.datetime.fromtimestamp(d_2[i])))
           myd_2=np.array(myd_2)
    except:
        if any(d_2):
            myd_2=[]
            for i in range(len(d_2)):
                myd_2.append(md.date2num(datetime.datetime.fromtimestamp(d_2[i])))
            myd_2=np.array(myd_2)
    hour=md.HourLocator(interval=2)
    #fmt=md.DateFormatter('%H:%M')
    fmt=md.DateFormatter('%H')
    h=height
    bs=BS
    myd=np.array(myd)

    #print myd[-1]
    #print myd_2[-1]

    fig=plt.figure(figsize=(6.5,9.0))
    #ax1=fig.add_subplot(311)
    ax1 = plt.subplot2grid((5, 2), (0, 0), colspan=2, rowspan=2)
    #plt.pcolormesh(myd,h,bs.T,vmin=vmin,vmax=vmax,cmap=plt.cm.get_cmap('Greens_r'),norm=colors.PowerNorm(gamma=2))
    plt.pcolormesh(myd,h,bs.T,vmin=tl,vmax=th,cmap=plt.cm.get_cmap('terrain'))
#    plt.pcolormesh(myd,h,bs.T,vmin=tl,vmax=th,cmap=plt.cm.get_cmap('gist_ncar'))
    pc=plt.colorbar()
    pc.set_label(r'$log10\beta$'+r' $[m^{-1} sr^{-1}]$')
    #fig.xticks(rotation=0)
    ax1.xaxis.set_major_locator(hour)
    ax1.xaxis.set_major_formatter(fmt)
    fig.autofmt_xdate()
    for y in range(len(threshold_array)):
        if any(myd_2):
            ax1.plot(myd_2,t_a[y],color=color[y])
            #ax1.plot(myd_2,t_a[y],'o',color=color[y],markersize='1.5')
        else:
            #ax1.plot(myd,t_a[y],'o',color=color[y],markersize='1.5')
            ax1.plot(myd,t_a[y],color=color[y])

    ax1.set_ylabel('Height m (AGL)')
    ax1.set_ylim(0, np.max(h))
    try:
        plt.title(stntitlenames[str(mwsiteid)]+' - '+dt)
    except:
        plt.title(str(mwsiteid)+' - '+dt)


# plotting only half in vertical for the smoothed plot
    sbs = bs_smoothed[:,:lht2]
    ax3 = plt.subplot2grid((5, 2), (2, 0), colspan=2, rowspan=2)
    #ax3 = fig.add_subplot(312)
    plt.pcolormesh(myd_2, h[:lht2], sbs.T, vmin=tl, vmax=th, cmap=plt.cm.get_cmap('terrain'))
    #plt.pcolormesh(myd_2, h, sbs.T, vmin=vmin, vmax=vmax, cmap=plt.cm.get_cmap('gist_ncar'))
    pc = plt.colorbar()  # label=r'$\beta$'+r' $[m^{-1} sr^{-1}]$')
    pc.set_label(r'$log10\beta$' + r' $[m^{-1} sr^{-1}]$')
    fig.autofmt_xdate()
    #fig.xticks(rotation=0)

    for y in range(len(threshold_array)):
        ax3.plot(myd_2,t_a[y], color=color[y])
        #ax3.plot(myd_2,t_a[y], 'o',color=color[y],markersize='1.5')
    ax3.set_ylabel('Height m (AGL)')
    ax3.set_title('Smoothed data')
    # ax3.axis('tight')
    ax3.xaxis.set_major_locator(hour)
    ax3.xaxis.set_major_formatter(fmt)
    ax3.axis('tight')
    fig.autofmt_xdate()

    #still not generalized enough

    no_bin = np.rint((th - tl)/p_wid + 1)
    p_wid_2 = p_wid / 2.
    p_wid2 = 2. * p_wid
    p_wid3 = 3. * p_wid
    a = np.linspace(tl,th, no_bin)
    #print "percent_first",percent_first[1:-1]
    ax2 = plt.subplot2grid((5, 2), (4, 0))
    #ax2.bar(a[1:-1], percent_first[1:-1], width=p_wid, align='center')
    ax2.bar(a[1:-1], percents[1:-1], width=p_wid, align='center')
    for i in range(len(multi)):
        where_art=np.where(multi[i]== a)
        #plt.axvline(multi[i]-0.025, color=color[i], label=str(multi[i])+',%0.02f'%percents[where_art[0][0]]+r'$\%$')
        ax2.axvline(multi[i]+p_wid_2, color=color[i], label=str(multi[i]))

    #where_art=np.where(multi[0]== a)
    #print "where_art",where_art,multi[0],where_art[0][0]
    #ax2.axvline(multi[0]+p_wid_2, color=color[0], label=str(multi[0]))
    #plt.axvline(multi[0]-0.025, color=color[0], label=str(multi[0])+',%0.02f'%percents[where_art[0][0]]+r'$\%$')
    #ax2.set_ylim([0,12.5])
    ax2.set_ylim([0,20.0])
    ax2.set_ylabel(r'$\%$')
    ax2.set_xlabel(r'$log10\beta$' + r' $[m^{-1} sr^{-1}]$')
    if unit_id == 'UUCL1' or unit_id == 'UUCL0' or unit_id == 'UUCLA' or unit_id == 'UUCLB':
        alab = np.linspace(tl+p_wid2,th, 12)
        #print 'alab', alab
        ax2.set_xticks(alab)
        ax2.set_xticklabels(['-7.6','','-7.2','', '-6.8', \
        '','-6.4','','-6.0','','-5.6', ''])
    else:
        #for logan
        alab = np.linspace(-8.0,-6.2,7)
        ax2.set_xticks(alab)
        ax2.set_xticklabels(['-8.0','-7.7','-7.4','-7.1', \
         '-6.8','-6.5','-6.2']) 


    ax4 = plt.subplot2grid((5, 2), (4, 1))
    #ax4.bar(a[1:-1], percents[1:-1], width=0.05, align='center')
    ax4.bar(a[2:-1], gradient, width=p_wid, align='center')
    for i in range(len(multi)):
        where_art=np.where(multi[i]== a)
        #plt.axvline(multi[i]-0.025, color=color[i], label=str(multi[i])+',%0.02f'%percents[where_art[0][0]]+r'$\%$')
        ax4.axvline(multi[i]+p_wid_2, color=color[i], label=str(multi[i]))
    #ax4.set_ylim([0,25.0])
    ax4.set_ylim([-5.0,5.0])
    ax4.set_ylabel(r'$\%$ diff')
    ax4.set_xlabel(r'$log10\beta$' + r' $[m^{-1} sr^{-1}]$')
    if unit_id == 'UUCL1' or unit_id == 'UUCL0' or unit_id == 'UUCLA' or unit_id == 'UUCLB':
        ax4.set_xticks(alab)
        ax4.set_xticklabels(['-7.6','','-7.2','', '-6.8', \
        '','-6.4','','-6.0','','-5.6','' ])
    else:
        #for logan
        ax4.set_xticks(alab)
        alab = np.linspace(-8.0,-6.2,7)
        ax2.set_xticklabels(['-8.0','-7.7','-7.4','-7.1', \
         '-6.8','-6.5','-6.2']) 

    plt.tight_layout()

    #use last date
    timel = time.strftime('%Y%m%d',time.localtime(date_array[-1]))
    #print 'timel',timel
    pngfilename = '/uufs/chpc.utah.edu/common/home/horel-group9/uunet/ceilometer/'+str(mwsiteid)+'/'+str(mwsiteid)+'_'+timel+'_bs_plot.png'
    print "Image File Name will be "+pngfilename
    try:
        plt.savefig(pngfilename)
    except:
        print "COULD NOT SAVE FIGURE TO "+pngfilename
    if realtime:
        pngcurrentfilename = '/uufs/chpc.utah.edu/common/home/horel-group9/uunet/ceilometer/'+str(mwsiteid)+'_current.png'
        try:
            plt.savefig(pngcurrentfilename)
        except:
            print "COULD NOT SAVE FIGURE TO "+pngcurrentfilename

    # plt.show()
    plt.close()