def graphMembersOverTimeWithPlusMinusText(membertimedata, memberinfos, membership):
    ydates, yvalues = zip(*membertimedata)
    ydates = list(ydates)
    #yvalues = map(lambda x:x[0],yvalues)
    plotlabel = [u"Number of members over time","Membership Income over time"]
    plt.plot(ydates, yvalues, 'o',linewidth=2, markevery=1)
    plt.ylabel("#Members")
    plt.xlabel("Month")
    plt.grid(True)
    plt.legend(plotlabel,loc='upper left')
    plt.twiny()
    plt.ylabel("Euro")
    #plt.title(plotlabel)
    ## label with +x-y members per month
    membersinmonth = dict(membertimedata)
    #print "\n".join([ "%s:%s" % (x[0],str(x[1])) for x in extractPlusMinusFromMemberships(membership).items()])
    pm_dates, pm_fee_dates = extractPlusMinusFromMemberships(membership)
    for astrdate, tpl in pm_dates.items():
        adate = datetime.datetime.strptime(astrdate, dateformat_monthonly_).date()
        assert(adate.day==1)
        if adate in membersinmonth:
            xy = (adate, membersinmonth[adate][0])
            xytext = (xy[0], xy[1]+1)
            plt.annotate(str(tpl), xy=xy, xytext=xytext,arrowprops=dict(facecolor='gray', shrink=0.5))
    for astrdate, tpl in pm_fee_dates.items():
        adate = datetime.datetime.strptime(astrdate, dateformat_monthonly_).date()
        assert(adate.day==1)
        if adate in membersinmonth:
            xy = (adate, membersinmonth[adate][1])
            xytext = (xy[0], xy[1]+30)
            plt.annotate(str(tpl), xy=xy, xytext=xytext,arrowprops=dict(facecolor='gray', shrink=0.5))
    plt.subplots_adjust(left=0.06, bottom=0.05, right=0.99, top=0.95)

def plt_twin( axis, tick0=None, tick=None ) : 
	'''
	Add x-top or y-right axis

	axis:
		['x' | 'y']

	tick0:
		Must between [0, 1]
	'''
	if (str(axis).lower() not in ['x', 'y']) : Raise(Warning, "axis='"+str(axis)+"' not in ['x', 'y']. Do nothing !")
	axis = str(axis).lower()
	#--------------------------------------------------
	if (tick0 is not None and tick is not None) : 
		tick0 = npfmt(tick0, float)
		if (tick0.min()<0 or tick0.max()>1) : 
			Raise(Warning, 'tick0.(min,max)=(%.1f, %.1f) out of [0, 1]. Do nothing !' % (tick0.min(), tick0.max()))
		else : 
			if   (axis == 'x') : 
				plt.twiny()
				plt.xticks(tick0, tick)
			elif (axis == 'y') : 
				plt.twinx()
				plt.yticks(tick0, tick)
	#--------------------------------------------------
	elif (tick0 is None and tick is None) : 
		if   (axis == 'x') : plt.tick_params(axis='x', which='both', labeltop='on', labelbottom='on')
		elif (axis == 'y') : plt.tick_params(axis='y', which='both', labelleft='on', labelright='on')
	#--------------------------------------------------
	else : Raise(Warning, 'tick0, tick must both ==None or !=None, now one is None but the other is not. Do nothing !')

def classificate(X, bins=None, n_components=5, verbose=True,
        histogram_ylabel='X',
        histogram_xlabel='Frequency',
        gaussian_xlabel='Probability',
        scatter_xlabel='Point No.'):
    # create histogram data
    hist, bins = statistics_histogram.histogram(X, bins=bins)
    # create gaussian mixture model
    model, AIC, BIC = clustering_gaussian.fit(X, n_components)

    # plot histogram
    pl.subplot(121)
    plot_histogram.histogram(hist, bins, transposition=True, color='k', alpha=0.7)
    pl.xlabel(histogram_xlabel)
    pl.ylabel(histogram_ylabel)
    pl.grid()

    # plot gaussian
    pl.twiny()
    plot_gaussian.gaussian(X, model, transposition=True)
    pl.xlabel(gaussian_xlabel)
    pl.legend()

    # reverse y axis
    ymin, ymax = pl.ylim()
    pl.ylim(ymax, ymin)

    # plot classified scatter
    pl.subplot(122)
    plot_gaussian.gaussian_scatter(X, model)
    pl.ylim(ymax, ymin)
    pl.xlabel(scatter_xlabel)
    pl.ylabel(histogram_ylabel)
    pl.minorticks_on()
    pl.legend()
    pl.grid(which='major', alpha=0.5)
    pl.grid(which='minor', alpha=0.2)

    # print result in stdout if verbose
    if verbose:
        print histogram_xlabel
        print "N = %d, Bins = %d" % (len(X), len(bins))
        # individual fitting curve
        properties = zip(model.weights_, model.means_, model._get_covars())
        properties = sorted(properties, key=lambda x: x[1])
        for (weight, mean, covar) in properties:
            var = np.diag(covar)[0]
            # create formula of normal deviation
            formula_label = textutils.gaussian_formula(len(X), mean, var)
            formula_label = "%.1f%%: %s" % (weight * 100, formula_label)
            print formula_label.encode('utf-8')

    # show figure
    pl.show()

 def dualLIFPlot(self, y = None):
     '''
         Makes a vertical LIF photon / s rate plot using internal lif
         averages.
     '''
     if y is None:
         y = self.dye_v - self.diode_v
     l1 = plt.plot(np.array([np.mean(x) for x in self.lif_avgs1]), y,
             color='red', label='PMT1')
     plt.xticks(color='red')
     plt.xlabel('ave photons / demod sample (photon kHz)', color='red')
     plt.title('PMT1 & 2 LIF Signal', y=1.08)
     plt.ylim(np.min(y), np.max(y))
     plt.twiny()
     l2 = plt.plot(np.array([np.mean(x) for x in self.lif_avgs2]), y,
             color='red', label='PMT2')

def plot_age(p, form, fof_np, age_labels, symbls):
    lims = 10. ** np.linspace(9, 15, (15.25-9) / .25) #These are the limits for the mass bins
    for (i, form_i) in enumerate(form):
        (avg, med) = binning(fof_np[i], form_i, lims)
        plt.semilogx(avg[0], avg[1], marker = symbls[i], subsx = [2, 3, 4, 5, 6, 7, 8, 9], \
                     label = 'Average {0}'.format(age_labels[i]))
        plt.semilogx(med[0], med[1], linestyle = 'dashed', label = 'Median {0}'.format(age_labels[i]))
    #MS = 13.5795 - 10.3112 * np.arcsinh(0.0504329 * avg[0] ** 0.08445)
    #plt.semilogx(avg[0], MS, 'r', label = 'MS Curve')
    plt.legend()
    plt.xlabel('M [M_sun / h]')
    plt.ylabel('Formation Age Lookback time [Gyr]')
    ##Create Second x axis
    xlims = plt.xlim()
    x2 = plt.twiny()
    x2.set_xscale('log')
    x2.set_xlim((xlims[0] / .73, xlims[1] / .73))
    x2.set_xlabel('M [M_sun]')
    
    ##Create second y axis
    
    ylims = plt.ylim()
    yt = plt.yticks()[0]
    y2 = plt.twinx()
    y2.set_ylim(ylims)
    y2.set_yticks(yt)
    yl = []
    for j in yt:
        if j < 13.5795:
            yl.append('{:5.2f}'.format(1.44224957031 / mth.sinh(0.0969815 * (13.5795 - j))**(2./3.)))
        else:
            yl.append('')
    y2.set_yticklabels(yl)
    y2.set_ylabel('z + 1')

def plotBoundaries (tranNames, blocks):
    '''Plot exon boundaries as vertical lines.'''

    ticksStart  = list()     # x-axis tick positions in phony space
    ticksEnd    = list()
    labelsStart = list()     # x-axis labels are real genomic coordinates of block start
    labelsEnd   = list()

    tickSep = blocks[-1].boundary / MAX_LABELS     # separation required to avoid label overlap

    frompos = 0
    for bound in blocks:

        plt.vlines (bound.boundary, 0, len(tranNames)+1, linewidth=0.5, colors='red')    # block boundary

        if not bound.annot:            # if block does not include annotation exon(s), make it blush
            plt.axvspan(frompos, bound.boundary, facecolor='Pink', alpha=0.4)

        if len(ticksStart) == 0 or frompos - ticksStart[-1] > tickSep:     # add tick only if there's room for the label
            ticksStart.append(frompos)
            labelsStart.append(str(bound.start))
        if len(ticksEnd) == 0 or bound.boundary - ticksEnd[-1] > tickSep:
            ticksEnd.append(bound.boundary)
            labelsEnd.append(str(bound.end))

        frompos = bound.boundary           # new block start

    numberedNames = list()                 # cluster name + line number
    for ix, name in enumerate(tranNames):
        numberedNames.append('%s %3d' % (name, ix+1))

    plt.grid(axis='y')                     # turn on horizontal dotted lines
    plt.xlim (0, blocks[-1].boundary)
    plt.ylim (len(tranNames)+1, 0)
    plt.xticks (ticksStart, labelsStart, fontsize='xx-small')
    plt.yticks (xrange(1,len(tranNames)+2), numberedNames, fontsize='xx-small')

    # These lines came from an obscure posting on StackOverflow. They
    # create a second X axis at the bottom of the figure. I haven't
    # got a clue how they work, but they do.

    x2 = plt.twiny()
    x2.set_frame_on(True)
    x2.patch.set_visible(False)
    x2.xaxis.set_ticks_position('bottom')
    x2.xaxis.set_label_position('bottom')
    x2.spines['bottom'].set_position(('outward', 20))

    plt.axes(x2)

    plt.grid(axis='y')
    plt.xlim (0, blocks[-1].boundary)
    plt.ylim (len(tranNames)+1, 0)
    plt.xticks (ticksEnd, labelsEnd, fontsize='xx-small')
    plt.yticks (xrange(1,len(tranNames)+2), numberedNames, fontsize='xx-small')

    return

def plotPatternBeyondRepeat(genomeSource1,genomeSource2,  start1,  start2, plotRange):
    f1 = open(genomeSource1,'r')
    f2 = open(genomeSource2,'r')
    
    pointerLocation1 = start1
    pointerLocation2 = start2
    
    windowSize = 10
    
    defaultsize = 24
    distanceList = []
    plotRange = plotRange/windowSize
    
    for index in range(plotRange):
        f1.seek(pointerLocation1)
        f2.seek(pointerLocation2)
    
        str1 = f1.read(windowSize)
        str2 = f2.read(windowSize)
        
        pointerLocation1 = pointerLocation1 + windowSize
        pointerLocation2 = pointerLocation2 + windowSize
        
        distance = hammingDistance(str1, str2, min(windowSize,len(str1),len(str2)))
        distanceList.append(distance)        
        
    
    #plt.subplot(111)
    plt.plot(range(0,len(distanceList)*windowSize, windowSize), distanceList)
    
    plt.tick_params(axis='both', which='major', labelsize=defaultsize)
    plt.tick_params(axis='both', which='minor', labelsize=defaultsize)
    
    plt.xlabel("Genomics location of Copy 1 ( unit : base )", fontsize=defaultsize)    
    plt.ylabel("Hamming distance within a window of length 10 ",fontsize=defaultsize)
    
    #print plt.xticks()
    locs =range(0,len(distanceList)*windowSize+1, len(distanceList)*windowSize/4)
    tick_lbls =  range(start1, start1+plotRange*windowSize+1, plotRange*windowSize/4)
    plt.xticks(locs, tick_lbls, fontsize=defaultsize)
    
    ax2 = plt.twiny()
    plt.tick_params(axis='both', which='major', labelsize=defaultsize)
    plt.tick_params(axis='both', which='minor', labelsize=defaultsize)
    
    plt.xlabel("Genomics location of Copy 2 ( unit : base )", fontsize=defaultsize)    
    plt.ylabel("Hamming distance within a window of length 10 ",fontsize=defaultsize)
    
    tick_locs = range(0, plotRange*windowSize+1, plotRange*windowSize/4)
    tick_lbls = range(start2, start2+plotRange*windowSize+1, plotRange*windowSize/4)
    plt.xticks(tick_locs, tick_lbls,fontsize=defaultsize)
  
    plt.show()

    f1.close()
    f2.close()

def make_twiny(offset):
    """
    """
    
    ax3 = plt.twiny()
    ax3.xaxis.set_ticks_position('bottom')
    ax3.spines['bottom'].set_position(('axes',-offset))
    ax3.set_frame_on(True)
    ax3.patch.set_visible(False)
    ax3.spines["bottom"].set_visible(True)
    
    return ax3

def case_three(filename):
    with open(filename) as f:
        data = map(lambda line: map(int, line.split(',')), f.readlines())

    plots = {}
    for d in data:
        lock = d[3]
        if lock in plots:
            plots[lock][0].append(d[2])  # c
            plots[lock][1].append(d[4])  # time
        else:
            plots[lock] = [[d[2]], [d[4]]]

    for (lock, (x, y)) in plots.items():
        if lock == 1:
            lock_name = 'Pthread Mutex Lock'
        elif lock == 2:
            lock_name = 'Pthread Spin Lock'
        elif lock == 3:
            lock_name = 'TAS Spin Lock'
        elif lock == 4:
            lock_name = 'TTAS Spin Lock'
        elif lock == 5:
            lock_name = 'Exponential Backoff Lock'
        else:
            lock_name = 'Queue Lock'
        pyplot.plot(range(len(x)), y, label='%s' % lock_name)

    pyplot.legend(loc='best')
    pyplot.xticks(range(11), (1000, 900, 800, 700, 600, 500, 400, 300, 200, 100, 0))
    pyplot.xlabel('Operations Inside Critical Section')
    pyplot.twiny()
    pyplot.xticks(range(11), (0, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000))
    pyplot.xlabel('Operations Outside Critical Section')
    pyplot.ylabel('Running Time (ms)')
    pyplot.savefig('case_three.png')
    pyplot.clf()

def labelnai(color='gray', zorder=1, twins=True, **kwargs):
    importmpl()
    import gbm
    for (i, (detp, dett)) in enumerate(gbm.naipt.T):
        plt.text(detp, dett, gbm.nlist[i], horizontalalignment='center', verticalalignment='center', color=color, zorder=zorder, *kwargs)
    plt.xlabel(r'spacecraft $\phi$ [rad]')
    plt.ylabel(r'spacecraft $\theta$ [rad]')
    # plt.subplots_adjust(top=0.875)
    plt.setp(plt.gca(), xlim=[0, 2*np.pi], ylim=[np.pi, 0])
    plt.gca().xaxis.tick_bottom()
    plt.gca().yaxis.tick_left()
    if twins:
        ax2 = plt.twinx()
        plt.setp(ax2, ylim=[180, 0])
        ay2 = plt.twiny()
        plt.setp(ay2, xlim=[0, 360])

    def setup_axes(self):
        self.ax = plt.subplot(111)

        # Second Y axis on the right for the ZHR 
        self.ax_zhr = plt.twinx(ax=self.ax)
        self.ax_zhr.set_ylabel("ZHR (r={0:.1f}, $\gamma$={1:.2f})".format(
                                                self.profile.popindex,
                                                self.profile.gamma))
        self.ax_zhr.yaxis.set_major_formatter(plt.FuncFormatter(self.zhr_formatter))
        
        # Second X axis as top for the solar longitude
        self.ax2 = plt.twiny(ax=self.ax)
        self.ax2.set_xlabel("Solar longitude (J2000.0)")
        self.ax2.xaxis.set_major_formatter(plt.FuncFormatter(self.sollon_formatter))
          
        self.ax.grid(which="both")

 def _create_figure():
     # create the figure with four subplots with different size
     # - 1st is for the predominant melody and performed notes
     # - 2nd is the pitch distribution and note models, it shares the y
     # axis with the 1st
     # - 3rd is the melodic progression, it shares the x axis with the 1st
     # - 4th is for the sections, it is on top the 3rd
     fig = plt.figure()
     gs = gridspec.GridSpec(2, 2, width_ratios=[6, 1], height_ratios=[4, 1])
     ax1 = fig.add_subplot(gs[0])  # pitch and notes
     ax2 = fig.add_subplot(gs[1], sharey=ax1)  # pitch dist. and note models
     ax4 = fig.add_subplot(gs[2])  # sections
     ax5 = fig.add_subplot(gs[3])  # makam, tempo, tonic, ahenk annotations
     ax3 = plt.twiny(ax4)  # melodic progression
     ax1.get_shared_x_axes().join(ax1, ax3)
     fig.subplots_adjust(hspace=0, wspace=0)
     return fig, ax1, ax2, ax3, ax4, ax5

def main():


	plt.figure(figsize=(6, 5), dpi=80)
	
	wyk1 = plt.subplot(121)
	
	rysowanie('rsel.csv', 'b')
	rysowanie('cel-rs.csv', 'g')
	rysowanie('2cel-rs.csv', 'r')
	rysowanie('cel.csv', 'k')
	rysowanie('2cel.csv', 'm')  

	plt.grid(linewidth=0.5, color='grey')
	plt.axis((0,500,60,100), size="small")
	wyk1.tick_params(labelsize="small")

	plt.xlabel('Rozegranych gier (x1000)', size="small")
	plt.ylabel('Odsetek wygranych gier [%]', size="small")	

	legenda = plt.legend(['1-Evol-RS','1-Coev-RS','2-Coev-RS','1-Coev','2-Coev'], loc="lower right", prop={'size':'x-small'}, fancybox=True)

	legenda.get_frame().set_alpha(0.75)
	legenda.get_frame().set_linewidth(0.5)

	wyk2 = plt.twiny()
	plt.axis([0,200,60,100], size="small")
	wyk2.set_xticks([0,40,80,120,160,200])
	wyk2.tick_params(labelsize="small")
	plt.xlabel('Pokolenie', size="small")


	names_list=['1-Evol-RS','1-Coev-RS','2-Coev-RS','1-Coev','2-Coev']
	
	plt.subplot(122)
	#plt.boxplot(box_list, notch = True, labels = names_list, showmeans=True, meanprops=dict(marker='o'))
	frame = plt.gca()
	plt.grid()
	frame.axes.get_yaxis().set_visible(False)
	plt.twinx()
	plt.grid()
	plt.axis([0,6,60,100])
	
	
	plt.savefig('myplot.jpg')
	plt.close()

def draw_left(m_list, ox_list, series):
    s1 = plt.subplot(1, 2, 1)
    [i.set_linewidth(0.5) for i in s1.spines.itervalues()]
    plt.xlim(xmax=500)
    plt.axis(siez='small')
    plt.grid(True, alpha=0.5, linewidth=0.3)
    plt.tick_params(labelsize='x-small')
    i = 0
    for c, l, mark in series:
        plt.plot(ox_list[i], m_list[i], color=c, label=l, marker=mark, markersize=5, markevery=25, alpha=0.8)
        i += 1
    legend = plt.legend(loc = 'lower right', framealpha=0.5, fontsize='small', fancybox=True, prop={'size': 'x-small'})
    legend.get_frame().set_linewidth(0.5)
    plt.xlabel('Rozegranych gier ($\\times$ 1000)', size='x-small')
    plt.ylabel('Odsetek wygranych gier [$\%$]', size='x-small')
    s2 = plt.twiny()
    plt.xlabel('Pokolenie', size='x-small')
    s2.set_xticks([0, 40, 80, 120, 160, 200])
    s2.tick_params(labelsize='x-small')

 def figure1(self, count=50):
   tha = self.theta_a
   thc = self.theta_c
   plt.xlabel(r'$t_c \mathrm{[nm]}$', size=15)
   plt.ylabel(r'$I_a/I_c$', size=15)
   plt.ylim(0, 60)
   plt.text(0.05, 50,
            '$T_{\mathrm{Fix}}=$'+str(self.T)+' $\mathrm{nm}$',
            horizontalalignment='left',
            verticalalignment='bottom',
            fontsize=20)
   ta = np.linspace(0, self.T, count)
   x = ta
   y = []
   for t1 in ta:
     val = self._get_I(t1, self.T-t1, tha, thc)
     y.append(val)
   plt.plot(x, y)
   ax = plt.twiny()
   ax.set_xlabel(r'$t_a \mathrm{[nm]}$', size=15)
   ax.set_xlim(self.T, 0)
   plt.show()

 def makePlot(self):
     minx = np.min(self.powers)
     maxx = np.max(self.powers)
     xvals = np.linspace(minx, maxx, self.powers.size*100)
     yvals = self.fitfunc(xvals)
     pyplot.figure()
     pyplot.plot(self.powers, self.dataFluor, 'o', markersize=4)
     pyplot.plot(xvals,yvals)
     pyplot.ylim(ymin = 0)
     s = "Pulsed Axial Power Scan {}\n".format(self.name)
     s = s + "Detuning {0:.2f} Gamma\n".format(self.detuning)
     maxsat = maxx / self.xscalefit
     s = s +  "Max Saturation {0:.0f}\n".format(maxsat)
     s = s + "Power for s = 1 : {0:.2f} microwatt\n".format(1000 * self.xscalefit)
     Isat = 4.33*10**-7 #mW / micron^2, #see e.g aarhus phd peter staanum
     waist = np.sqrt(self.xscalefit/ Isat) / 2.0
     s = s + "Waist of ~{0:.0f} micron".format(waist)
     pyplot.text(0.5, 10000.0 , s )#,  verticalalignment='top')
     pyplot.xlabel('397 Power mW')
     pyplot.ylabel('Fluorescence Counts/Sec')
     ax = pyplot.twiny()
     pyplot.xlabel("Saturation =  Omega / Gamma")
     pyplot.xlim(xmin = minx / self.xscalefit, xmax = maxx/ self.xscalefit)
     pyplot.show()

def sfh(sim, filename=None, massform=True, clear=False, legend=False,
        subplot=False, trange=False, bins=100, **kwargs):
    '''
    star formation history

    **Optional keyword arguments:**

       *trange*: list, array, or tuple
         size(t_range) must be 2. Specifies the time range.

       *bins*: int
         number of bins to use for the SFH

       *massform*: bool
         decides whether to use original star mass (massform) or final star mass

       *subplot*: subplot object
         where to plot SFH

       *legend*: boolean
         whether to draw a legend or not

       *clear*: boolean
         if False (default), plot on the current axes. Otherwise, clear the figure first.

    By default, sfh will use the formation mass of the star.  In tipsy, this will be
    taken from the starlog file.  Set massform=False if you want the final (observed)
    star formation history

    **Usage:**

    >>> import pynbody.plot as pp
    >>> pp.sfh(s,linestyle='dashed',color='k')


    '''

    if subplot:
        plt = subplot
    else:
        import matplotlib.pyplot as plt

    if "nbins" in kwargs:
        bins = kwargs['nbins']

    if 'nbins' in kwargs:
        bins = kwargs['nbins']
        del kwargs['nbins']

    if ((len(sim.g)>0) | (len(sim.d)>0)): simstars = sim.star
    else: simstars = sim

    if trange:
        assert len(trange) == 2
    else:
        trange = [simstars['tform'].in_units(
            "Gyr").min(), simstars['tform'].in_units("Gyr").max()]
    binnorm = 1e-9 * bins / (trange[1] - trange[0])

    trangefilt = filt.And(filt.HighPass('tform', str(trange[0]) + ' Gyr'),
                          filt.LowPass('tform', str(trange[1]) + ' Gyr'))
    tforms = simstars[trangefilt]['tform'].in_units('Gyr')

    if massform:
        try:
            weight = simstars[trangefilt][
                'massform'].in_units('Msol') * binnorm
        except (KeyError, units.UnitsException):
            warnings.warn(
                "Could not load massform array -- falling back to current stellar masses", RuntimeWarning)
            weight = simstars[trangefilt]['mass'].in_units('Msol') * binnorm
    else:
        weight = simstars[trangefilt]['mass'].in_units('Msol') * binnorm

    if clear:
        plt.clf()
    sfhist, thebins, patches = plt.hist(tforms, weights=weight, bins=bins,
                                        histtype='step', **kwargs)
    if not subplot:
        # don't set the limits
        #plt.ylim(0.0, 1.2 * np.max(sfhist))
        plt.xlabel('Time [Gyr]', fontsize='large')
        plt.ylabel('SFR [M$_\odot$ yr$^{-1}$]', fontsize='large')
    else:
        plt.set_ylim(0.0, 1.2 * np.max(sfhist))

    # Make both axes have the same start and end point.
    if subplot:
        x0, x1 = plt.get_xlim()
    else:
        x0, x1 = plt.gca().get_xlim()


    # add a z axis on top if it has not been already done by an earlier plot:
    from pynbody.analysis import pkdgrav_cosmo as cosmo
    c = cosmo.Cosmology(sim=sim)
    old_axis = plt.gca()

    if len(plt.gcf().axes)<2:
        pz = plt.twiny()
#.........这里部分代码省略.........

#.........这里部分代码省略.........
                      frameon=False)
    #Add the lopsided and ellipticity constraints from Rix/Zaritsky
    if _ADDRIX:
        pyplot.errorbar([1./16.],[5.6],
                        yerr=[5.6/2.],
                        marker='d',color='0.6',
                        mew=1.5,mfc='none',mec='0.6')
        pyplot.errorbar([1./8./numpy.sqrt(2.)],[6.4],
                        yerr=numpy.reshape(numpy.array([6.4/0.045*0.02,6.4/0.045*0.03]),(2,1)),
                        marker='d',color='0.6',mec='0.6',
                        mew=1.5,mfc='none')
    if _ADDGCS:
        ks_gcs, psd_gcs, e_psd_gcs, gcsp= plot_psd_gcs()
    if _ADDRAVE:
        ks_rave, psd_rave, e_psd_rave, ravep= plot_psd_rave()
    if _ADDRED:
        plot_psd_red()
    l2= pyplot.legend((apop[0],ravep[0],gcsp[0]),
                      (r'$\mathrm{APOGEE}$',
                       r'$\mathrm{RAVE}$',
                       r'$\mathrm{GCS}$'),
                      loc='upper right',bbox_to_anchor=(.95,.750),
                      numpoints=1,
                      prop={'size':14},
                      frameon=False)
    pyplot.gca().add_artist(l1)
    pyplot.gca().add_artist(l2)
    #Plot an estimate of the noise, based on looking at the bands
    nks= numpy.linspace(2.,120.,2)
    if not 'mba23' in socket.gethostname():
        pyplot.fill_between(nks,[2.,2.],hatch='/',color='k',facecolor=(0,0,0,0),
                            lw=0.)
    #                        edgecolor=(0,0,0,0))
    pyplot.plot(nks,[2.,2.],color='k',lw=1.5)
    nks= numpy.linspace(0.17,2.,2)
    def linsp(x):
        return .8/(numpy.log(0.17)-numpy.log(2.))*(numpy.log(x)-numpy.log(0.17))+2.8
    if not 'mba23' in socket.gethostname():
        pyplot.fill_between(nks,linsp(nks),hatch='/',color='k',facecolor=(0,0,0,0),
                            lw=0.)
    #                        edgecolor=(0,0,0,0))   
    #Plot lines
    pyplot.plot([0.17,0.17],[0.,2.8],color='k',lw=1.5)
    pyplot.plot(nks,linsp(nks)+0.01,color='k',lw=1.5)
    bovy_plot.bovy_text(0.19,.5,r'$95\,\%\,\mathrm{noise\ range}$',
                        bbox=dict(facecolor='w',edgecolor='w'),fontsize=14.)
    if _INTERP:
        interpks= list(ks[:-5])
        interppsd= list((scale*numpy.sqrt(psd1d[1][1:-3]-_SUBTRACTERRORS*numpy.median(noisepsd,axis=0)))[:-5])
        interppsd_w= (scale*0.5*psd1d[2][1:-3]/numpy.sqrt(psd1d[1][1:-3]))[:-5]
        interppsd_w[:8]*= 0.025 #fiddling to get a decent fit
        interppsd_w[3:5]*= 0.001
        interppsd_w= list(interppsd_w)
        interpks.append(0.025)
        interppsd.append(10.**-5.)
        interppsd_w.append(0.00001)
        interpks.append(110.)
        interppsd.append(1.)
        interppsd_w.append(0.00001)
        if _ADDGCS:
            interpks.extend(ks_gcs)
            interppsd.extend(psd_gcs)
            interppsd_w.extend(e_psd_gcs)
        if _ADDRAVE:
            interpks.extend(ks_rave[5:])
            interppsd.extend(psd_rave[5:])
            interppsd_w.extend(e_psd_rave[5:])
        interpks= numpy.array(interpks)
        sortindx= numpy.argsort(interpks)
        interpks= interpks[sortindx]
        interppsd= numpy.array(interppsd)[sortindx]
        interppsd_w= interppsd/numpy.array(interppsd_w)[sortindx]
        interpindx= True-numpy.isnan(interppsd)
    #interpspec= interpolate.InterpolatedUnivariateSpline(interpks[interpindx],
        interpspec= interpolate.UnivariateSpline(numpy.log(interpks[interpindx]),
                                                 numpy.log(interppsd[interpindx]/3.),
                                                 w=interppsd_w,
                                                 k=3,s=len(interppsd_w)*0.8)
        pks= numpy.linspace(interpks[0],interpks[-1],201)
        #bovy_plot.bovy_plot(pks,
        #                    3.*numpy.exp(interpspec(numpy.log(pks))),
#                    'k-',overplot=True)
    def my_formatter(x, pos):
        return r'$%g$' % x
    def my_formatter2(x, pos):
        return r'$%g$' % (1./x)
    major_formatter = FuncFormatter(my_formatter)
    major_formatter2 = FuncFormatter(my_formatter2)
    ax= pyplot.gca()
    ax.xaxis.set_major_formatter(major_formatter)
    if not _PROPOSAL:
        ax2= pyplot.twiny()
        xmin, xmax= ax.xaxis.get_view_interval()
        ax2.set_xscale('log')
        ax2.xaxis.set_view_interval(1./xmin,1./xmax,ignore=True)
        ax2.set_xlabel('$\mathrm{Approximate\ scale}\,(\mathrm{kpc})$',
                       fontsize=12.,ha='center',x=0.5)
        ax2.xaxis.set_major_formatter(major_formatter)
    bovy_plot.bovy_end_print(plotfilename,dpi=300)
    return None

#.........这里部分代码省略.........
            xlabel='Temperature (K)'
            xvalues=self.Temperature
            desc_str=' Variable Density'

        # put all actual plotting here
        plt.ion()
        #  maxAll is an array
        ymax = np.max(intensity[:, topLines[0]]/maxAll)
        ymin = ymax
        plt.figure()
        ax = plt.subplot(111)
        nxvalues=len(xvalues)
        # reversing is necessary - otherwise, get a ymin=ymax and a matplotlib error
        for iline in range(top-1, -1, -1):
            tline=topLines[iline]
            plt.loglog(xvalues,intensity[:, tline]/maxAll)
            if np.min(intensity[:, tline]/maxAll) < ymin:
                ymin = np.min(intensity[:, tline]/maxAll)
            if np.max(intensity[:, tline]/maxAll) > ymax:
                ymax = np.max(intensity[:, tline]/maxAll)
            skip=2
            start=divmod(iline,nxvalues)[1]
            for ixvalue in range(start,nxvalues,nxvalues//skip):
                if ionNum == 1:
                    text = '%10.4f'%(wvl[tline])
                else:
                    text = '%s %10.4f'%(ionS[tline], wvl[tline])
                plt.text(xvalues[ixvalue], intensity[ixvalue, tline]/maxAll[ixvalue], text)
        plt.xlim(xvalues.min(),xvalues.max())
        plt.xlabel(xlabel,fontsize=fontsize)
        plt.ylabel(ylabel,fontsize=fontsize)
        if ndens == ntemp and ntemp > 1:
            plt.text(0.07, 0.5,title, horizontalalignment='left', verticalalignment='center', fontsize=fontsize,  transform = ax.transAxes)
            ax2 = plt.twiny()
            xlabelDen=r'Electron Density (cm$^{-3}$)'
            plt.xlabel(xlabelDen, fontsize=fontsize)
            plt.loglog(eDensity,intensity[:, topLines[top-1]]/maxAll, visible=False)
            ax2.xaxis.tick_top()
            plt.ylim(ymin/1.2, 1.2*ymax)
        else:
            plt.ylim(ymin/1.2, 1.2*ymax)
            plt.title(title+desc_str,fontsize=fontsize)
        plt.draw()
        #  need time to let matplotlib finish plotting
        time.sleep(0.5)

        # get line selection
        selectTags = []
        for itop in topLines:
            if ionNum == 1:
                selectTags.append(str(wvl[itop]))
            else:
                selectTags.append(ionS[itop]+ ' '+ str(wvl[itop]))
        numden = chGui.gui.choice2Dialog(selectTags)

        # num_idx and den_idx are tuples
        num_idx=numden.numIndex
        if len(num_idx) == 0:
            print(' no numerator lines were selected')
            return
        den_idx=numden.denIndex
        if len(den_idx) == 0:
            print(' no denominator lines were selected')
            return
        numIntens=np.zeros(len(xvalues),'Float64')
        for aline in num_idx:

 def _createPlot(self):
     bins = self._fluxdata.getBins()
     
     self._fig = plt.figure(figsize=(11,6), dpi=80) # 11*80 = 880 pixels wide !
     ax = plt.subplot(111)
     self._fig.subplots_adjust(0.1,0.17,0.92,0.87)
                    
     ax_zhr = plt.twinx(ax=ax)
     ax_zhr.set_ylabel("ZHR (r=%.1f, $\gamma$=%.2f)" % (self._fluxdata._popindex, self._fluxdata._gamma), fontsize=16)
     ax_zhr.yaxis.set_major_formatter(plt.FuncFormatter(self.zhr_formatter))
     
     ax2 = plt.twiny(ax=ax)
     ax2.set_xlabel("Solar longitude (J2000.0)", fontsize=16)
     ax2.xaxis.set_major_formatter(plt.FuncFormatter(self.sollon_formatter))
       
     ax.grid(which="both")
     
     if len(bins) > 0 and len(bins['time']) > 0:
         ax.errorbar(bins['time'], bins['flux'], yerr=bins['e_flux'], fmt="s", ms=4, lw=1.0, c='red' )    #fmt="+", ms=8    
     
     ax.set_xlim([self._begin, self._end])
     ax2.set_xlim([self._begin, self._end])
     ax_zhr.set_xlim([self._begin, self._end])
     
     # Determine the limit of the Y axis
     if self._ymax:
         my_ymax = self._ymax
     elif len(bins) == 0 or len(bins['time']) == 0:
     	my_ymax = 100
     else:
         my_ymax = 1.1*max(bins['flux']+bins['e_flux'])
     
     if len(bins) > 0:
         ax.set_ylim([0, my_ymax])
         ax2.set_ylim([0, my_ymax])
         ax_zhr.set_ylim([0, my_ymax])        
     
     
     if self._timespan > 90*24*3600: # 90 days
         """ More than 5 days: only show dates """
         majorLocator = mpl.dates.AutoDateLocator(maxticks=10)
         sollonLocator = majorLocator
         majorFormatter = mpl.dates.DateFormatter('%d %b')
         xlabel = "Date (UT, %s)" % self._begin.year   
                 
     elif self._timespan > 5*24*3600: # 5 days
         """ More than 5 days: only show dates """
         majorLocator = mpl.dates.AutoDateLocator(maxticks=10)
         sollonLocator = majorLocator
         majorFormatter = mpl.dates.DateFormatter('%d %b')
         xlabel = "Date (UT, %s)" % self._begin.year   
     
     elif self._timespan > 1*24*3600:
         """ Between 1 and 5 days: show hours """
         majorLocator = mpl.dates.HourLocator(byhour=[0])
         majorFormatter = mpl.dates.DateFormatter('%d %b')
         
         if self._timespan > 3*24*3600:
             t = 12
         elif self._timespan > 1.5*24*3600:
             t = 6
         else:
             t = 3
         byhour = np.arange(t, 24, t)
         
         minorLocator = mpl.dates.HourLocator(byhour=byhour)
         ax.xaxis.set_minor_locator(minorLocator)
         sollonLocator = mpl.dates.HourLocator(byhour=np.append(0, byhour))
         fmt2 = mpl.dates.DateFormatter('%H:%M')        
         ax.xaxis.set_minor_formatter(plt.FuncFormatter(fmt2))
         xlabel = "Date (UT, %s)" % self._begin.year   
         
     else:
         if self._timespan > 18*3600:
             majorLocator = mpl.dates.HourLocator( np.arange(0, 24, 3) ) 
             minorLocator = mpl.dates.HourLocator( np.arange(0, 24, 1) ) 
         elif self._timespan > 12*3600:
             majorLocator = mpl.dates.HourLocator( np.arange(0, 24, 2) ) 
             minorLocator = mpl.dates.HourLocator( np.arange(1, 24, 2) ) 
         elif self._timespan > 6*3600:
             majorLocator = mpl.dates.HourLocator( np.arange(0, 24, 1) ) 
             minorLocator = mpl.dates.MinuteLocator(  np.arange(0,60,30)  )
         elif self._timespan > 3*3600:
             majorLocator = mpl.dates.MinuteLocator(  np.arange(0,60,30)  )
             minorLocator = mpl.dates.MinuteLocator(  np.arange(0,60,15)  )
         elif self._timespan > 1*3600:
             majorLocator = mpl.dates.MinuteLocator( np.arange(0,60,15) ) 
             minorLocator = mpl.dates.MinuteLocator(  np.arange(0,60,5) )
         else:
             majorLocator = mpl.dates.MinuteLocator( np.arange(0,60,10) ) 
             minorLocator = mpl.dates.MinuteLocator(  np.arange(0,60,2) )
             
         
         majorFormatter = mpl.dates.DateFormatter('%H:%M')  
         ax.xaxis.set_minor_locator(minorLocator)
         sollonLocator = majorLocator
         xlabel = "Time (UT, %s)" % self._begin.strftime('%d %b %Y')
     
     ax.xaxis.set_major_formatter(plt.FuncFormatter(majorFormatter))    
     ax2.xaxis.set_major_locator(sollonLocator)
#.........这里部分代码省略.........