def chickling_pd_zoom(shotno, date=time.strftime("%Y%m%d")):
	
	fname, data = file_finder(shotno,date)
	
	data_1550 = data[0]['phasediff_co2'][100:]
	plot_time = np.linspace(0,1,data_1550.size)

	
	fig, ax = plt.subplots()
	ax.plot(plot_time, data_1550)
	ax.set_ybound(max(data_1550)+0.6, min(data_1550)-0.01)
	
	plt.title("Phase Difference for shot " + str(shotno) + " Date " +  str(date))
	plt.xlabel("Time, s")
	plt.ylabel("Phase Difference, Radians") 	
	
	x_zoom_bot = int(data_1550.size*(10/100))
	x_zoom_top = int(data_1550.size*(15/100))
	

	x1, x2, y1, y2 = 0.1, 0.15, max(data_1550[x_zoom_bot:x_zoom_top])+0.01, min(data_1550[x_zoom_bot:x_zoom_top])-0.01
	
	axins = inset_axes(ax, 4.3,1, loc=9)
	axins.plot(plot_time[x_zoom_bot:x_zoom_top], data_1550[x_zoom_bot:x_zoom_top])

	axins.set_xlim(x1, x2)
	if y1 < y2:
		axins.set_ylim(y1, y2)
	else:
		axins.set_ylim(y2, y1)

	mark_inset(ax, axins, loc1=2, loc2=4, fc="none", ec="0.5",lw=2)
	
	plt.show()

def QQplot(obs, Q, pos, color=None, ax=None, axins=True):
        if not color:
            color='blue'

        if ax is None:
            ax = plt.gca()

        mean = np.mean(obs, axis=0)
        # obs = np.random.multivariate_normal(mean, S, 3000)

        md2 = np.diag(np.dot(np.dot(obs - mean, Q), (obs -mean).T))
        sorted_md2 = np.sort(md2)
        v = (np.arange(1, obs.shape[0] + 1) - 0.375) / (obs.shape[0] + 0.25)
        quantiles = scipy.stats.chi2.ppf(v, df=obs.shape[1])

        # axins = inset_axes(ax, width="60%", height=1., loc=2)
        if axins:
            axins = inset_axes(ax, width="60%", height=1., loc=2)
            axins.axis(pos)

            axins.get_xaxis().tick_bottom()
            axins.get_yaxis().tick_left()
            # Remove the tick marks; they are unnecessary with the tick lines we just plotted.
            axins.tick_params(axis="both", which="both", bottom="off", top="off", labelbottom="off", left="off", right="off", labelleft="off")
            mark_inset(ax, axins, loc1=2, loc2=4, fc="none", ec="0.5")
            axins.xaxis.set_major_locator(MaxNLocator(nbins=1, prune='lower'))
            axins.scatter(quantiles, sorted_md2, color=colorscheme[color], alpha=0.3)
            axins.plot(quantiles, quantiles, color=colorscheme['green'], lw=2.5)


        ax.scatter(quantiles, sorted_md2, color=colorscheme[color], alpha=0.3)
        ax.plot(quantiles, quantiles, color=colorscheme['green'], lw=2.5)
        _clean_axes(ax)

def plot_linear_kernel_pairs(pairs):
    xdata, ydata = zip(*pairs)
    maxval = max(max(xdata), max(ydata))

    fig, ax = plt.subplots()
    ax.plot(xdata, ydata, '.')
    ax.plot([0, maxval], [0, maxval], '--')
    plt.xlabel("Linear Ridge Mean Absolute Error (kcal/mol)")
    plt.ylabel("Kernel Ridge Mean Absolute Error (kcal/mol)")

    # 15 is the zoom, loc is nuts
    axins = zoomed_inset_axes(ax, 15, loc=5)
    axins.plot(xdata, ydata, '.')
    axins.plot([0, maxval], [0, maxval], '--')

    # sub region of the original image
    axins.set_xlim(1, 6)
    axins.set_ylim(1, 6)

    # draw a bbox of the region of the inset axes in the parent axes and
    # connecting lines between the bbox and the inset axes area
    mark_inset(ax, axins, loc1=2, loc2=4, fc="none", ec="0.5")

    plt.draw()
    plt.show()

def show_insets(image, insets_coords, size=(6, 3)):
    fig = pyplot.figure(figsize=size)

    line_colors = colors.LinearSegmentedColormap.from_list("line_colors", ["red", "blue"], len(insets_coords))

    gs = gridspec.GridSpec(len(insets_coords), 2)
    main_ax = fig.add_subplot(gs[:, 0])

    main_extent = (0, image.shape[1], image.shape[0], 0)
    main_ax_img = main_ax.imshow(image, interpolation="nearest", extent=main_extent, origin="upper")
    main_ax_img.set_cmap("binary_r")
    pyplot.setp(main_ax.get_yticklabels(), visible=False)
    pyplot.setp(main_ax.get_xticklabels(), visible=False)

    for inset_index, inset_coords in enumerate(insets_coords):
        inset_ax = fig.add_subplot(gs[inset_index, 1])
        inset_ax_img = inset_ax.imshow(image, interpolation="nearest", extent=main_extent, origin="upper")
        inset_ax_img.set_cmap("binary_r")
        inset_ax.set_ylim(inset_coords[0].stop, inset_coords[0].start)
        inset_ax.set_xlim(inset_coords[1].start, inset_coords[1].stop)
        pyplot.setp(inset_ax.get_yticklabels(), visible=False)
        pyplot.setp(inset_ax.get_xticklabels(), visible=False)
        mark_inset(main_ax, inset_ax, loc1=1, loc2=3, fc="none", ec=line_colors(inset_index))

    gs.tight_layout(fig)

    return fig

    def inset_momentum_axes(cd):

        # TODO: This plot does not refresh correctly, skip the inset
        fig = mpl.figure(cd.plotfigure.figno)
        axes = fig.add_subplot(111)

        # Plot main figure
        axes.plot(cd.x, hu_1(cd), 'b-')
        axes.plot(cd.x, hu_2(cd), 'k--')
        axes.set_xlim(xlimits)
        axes.set_ylim(ylimits_momentum)
        momentum_axes(cd)

        # Create inset plot
        from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes
        from mpl_toolkits.axes_grid1.inset_locator import mark_inset
        inset_axes = zoomed_inset_axes(axes, 0.5, loc=3)
        inset_axes.plot(cd.x, hu_1(cd), 'b-')
        inset_axes.plot(cd.x, hu_2(cd), 'k--')
        inset_axes.set_xticklabels([])
        inset_axes.set_yticklabels([])
        x_zoom = [-120e3,-30e3]
        y_zoom = [-10,10]
        inset_axes.set_xlim(x_zoom)
        inset_axes.set_ylim(y_zoom)
        mark_inset(axes, inset_axes, loc1=2, loc2=4, fc='none', ec="0.5")
        # mpl.ion()
        mpl.draw()

def plot_target_pred(train_xdata, test_xdata, train_ydata, test_ydata, loc=5):
    maxval = max(train_xdata.max(), test_xdata.max(), train_ydata.max(), test_ydata.max())
    minval = min(train_xdata.min(), test_xdata.min(), train_ydata.min(), test_ydata.min())

    fig, ax = plt.subplots()
    ax.plot(train_xdata, train_ydata, '.', label="Train")
    ax.plot(test_xdata, test_ydata, '.', label="Test")
    plt.legend(loc="best")
    ax.plot([minval, maxval], [minval, maxval], '--')
    plt.xlabel("Target Value (kcal/mol)")
    plt.ylabel("Predicted Value (kcal/mol)")

    axins = zoomed_inset_axes(ax, 30, loc=loc) # 30 is zoom, loc is .... nuts
    axins.plot(train_xdata, train_ydata, '.', label="Train")
    axins.plot(test_xdata, test_ydata, '.', label="Test")
    axins.plot([minval, maxval], [minval, maxval], '--')
    # sub region of the original image
    middle = test_xdata.mean() - 170
    x1, x2, y1, y2 = -15+middle, 15+middle, -15+middle, 15+middle
    axins.set_xlim(x1, x2)
    axins.set_ylim(y1, y2)

    mark_inset(ax, axins, loc1=2, loc2=3, fc="none", ec="0.5")
    plt.draw()
    plt.show()

def plot_us(lats, lons, save_name=None):
    fig = plt.figure(figsize=(10, 10))
    ax = fig.add_subplot(111)
    big_map = Basemap(resolution='h',
                      lat_0=36, lon_0=-107.5,
                      llcrnrlat=32, llcrnrlon=-125,
                      urcrnrlat=43, urcrnrlon=-110)
    big_map.drawcoastlines()
    big_map.drawstates()
    big_map.drawcountries()
    big_map.drawmapboundary(fill_color='#7777ff')
    big_map.fillcontinents(color='#ddaa66', lake_color='#7777ff', zorder=0)
    x, y = big_map(lons, lats)
    big_map.plot(x[0], y[0], 'ro', markersize=2)

    axins = zoomed_inset_axes(ax, 20, loc=1)
    ll_lat, ll_lon = 37.8, -122.78
    ur_lat, ur_lon = 38.08, -122.43

    axins.set_xlim(ll_lon, ur_lon)
    axins.set_ylim(ur_lon, ur_lat)

    small_map = Basemap(resolution='h',
                        llcrnrlat=ll_lat, llcrnrlon=ll_lon,
                        urcrnrlat=ur_lat, urcrnrlon=ur_lon,
                        ax=axins)
    small_map.drawcoastlines()
    small_map.drawmapboundary(fill_color='#7777ff')
    small_map.fillcontinents(color='#ddaa66', lake_color='#7777ff', zorder=0)
    x, y = small_map(lons, lats)
    small_map.plot(x, y, 'ro', markersize=3)

    mark_inset(ax, axins, loc1=2, loc2=4, fc="none", ec="0.5")
    if save_name: 
        fig.savefig(save_name)

def render_input_spectrum():

    folder = "/cosma/home/durham/rhgk18/ls_structure/pks"
    bao    = folder+"/wig.txt"
    nbao   = folder+"/nowig.txt"

    pkbao  = np.loadtxt(bao)
    pknbao = np.loadtxt(nbao)

    fig, ax = plt.subplots()

    ax.loglog(pkbao[:,0]  ,pkbao[:,1]  , color='r', label='With BAO')
    ax.loglog(pknbao[:,0] ,pknbao[:,1] , color='b', label='Without BAO')
    ax.set_xlabel("Wavenumber, $k$ [$h/Mpc$]", size=28)
    ax.set_ylabel("Power, $P(k)$ [$(Mpc/h)^3$]", size=28)
    ax.legend(prop={'size':28})
    
    axins = zoomed_inset_axes(ax, 5, loc=3)
    axins.loglog(pkbao[:,0]  ,pkbao[:,1]  , color='r', label='iWith BAO')
    axins.loglog(pknbao[:,0] ,pknbao[:,1] , color='b', label='iWithout BAO')

    x1, x2, y1, y2 = 0.02, 0.1, 5000, 30000
    axins.set_xlim(x1, x2)
    axins.set_ylim(y1, y2)

    plt.xticks(visible=False)
    plt.yticks(visible=False)

    mark_inset(ax, axins, loc1=2, loc2=1, fc="none", ec="0.5")

    
    plt.show()

def plotBinned(cls_in,dcls_in,l_out,bins,output_prefix,title=None,theory=None,dtheory=None,delta=None,cosmic=None):
	fig=plt.figure(1)
	plt.clf()
        ax=fig.add_subplot(111)
        good_l=np.logical_and( l_out > 25 , l_out <= 250)

	if not (theory is None) :
		ax.plot(l_out,theory,'r-')
	if not (cosmic is None) :
		ax.fill_between(l_out,(theory-cosmic),(theory+cosmic),alpha=.5,facecolor='red')
		#plt.fill_between(l_out,(theory-dtheory),(theory+dtheory),alpha=.5,facecolor='red')
	if not (dtheory is None) :
		ax.errorbar(l_out,theory,yerr=dtheory,color='red')
	ax.errorbar(l_out,cls_in,yerr=dcls_in,color='black',fmt='k.',linestyle='None')
	if not (delta is None) :
		ax.fill_between(l_out,cls_in-delta,cls_in+delta,color='gray',alpha=0.5)
	#plt.xlim([0,np.max(l_out+bins)])
	#plt.ylim([np.min(b_cl['llcl']-b_dcl['llcl']),np.max(b_cl['llcl']+b_dcl['llcl'])])
	ax.set_xlabel('$\ell$')
	ax.set_ylabel('$\\frac{\ell(\ell+1)}{2\pi}C_{\ell}\ \\frac{\mu K^{2}}{m^{4}}$')
        #ax.set_xlim([0,np.max(l_out+bins)])
#        ax.autoscale(axis='y',tight=True)

	if title:
		ax.set_title(title)
	else:
		ax.set_title('Binned Cls {:02d}'.format(bins))

        axins = inset_axes(ax,width="50%",height="30%",loc=9)

	if not (theory is None) :
		axins.plot(l_out[good_l],theory[good_l],'r-')
	if not (cosmic is None) :
		axins.fill_between(l_out[good_l],(theory-cosmic)[good_l],(theory+cosmic)[good_l],alpha=.5,facecolor='red')
		#plt.fill_between(l_out,(theory-dtheory),(theory+dtheory),alpha=.5,facecolor='red')
	if not (dtheory is None) :
		axins.errorbar(l_out[good_l],theory[good_l],yerr=dtheory[good_l],color='red')
	axins.errorbar(l_out[good_l],cls_in[good_l],yerr=dcls_in[good_l],color='black',fmt='k.',linestyle='None')
	if not (delta is None) :
		axins.fill_between(l_out[good_l],(cls_in-delta)[good_l],(cls_in+delta)[good_l],color='gray',alpha=0.5)
        axins.set_xlim(25,255)
        axins.set_yscale('log')
        #axins.set_ylim(ymax=1e5)
        #axins.set_ylim(-1e5,1e5)

        mark_inset(ax,axins,loc1=2,loc2=4,fc='none',ec='0.1')
        #axins.set_xlim([25,250])
        axins.autoscale('y',tight=True)


        plt.draw()
        #plt.show(block=True)
	fig.savefig(output_prefix+'_linear_{:02d}.eps'.format(bins),format='eps')
	fig.savefig(output_prefix+'_linear_{:02d}.png'.format(bins),format='png')
        ax.set_yscale('log')
	fig.savefig(output_prefix+'_log_{:02d}.eps'.format(bins),format='eps')
	fig.savefig(output_prefix+'_log_{:02d}.png'.format(bins),format='png')

def plot_default_topo(bmap,x, y, data):
    assert isinstance(bmap, Basemap)
    fig = plt.figure(figsize=(10, 6))
    gs = GridSpec(1, 2, width_ratios=[4, 1], height_ratios=[3, 1])
    ax = fig.add_subplot(gs[0, 1])

    bmap.etopo()
    img = bmap.contourf(x, y, data, norm=LogNorm(), alpha=0.8, ax=ax)
    cb = bmap.colorbar(img)
    cb.ax.set_title(r"${\rm km^2}$")
    bmap.drawcoastlines(linewidth=0.3)

    bmap.drawmeridians(np.arange(-180, 180, 20))
    bmap.drawparallels(np.arange(-90, 110, 20))
    bmap.readshapefile("data/shp/wri_basins2/wribasin", "basin", color="k", linewidth=2, ax=ax)



    lower_left_lat_lon = (-50, 45)
    axins = fig.add_subplot(gs[0, 0])  # plt.axes([-0.15, 0.1, 0.6, 0.8])  # zoomed_inset_axes(ax, 2, loc=1)  # zoom = 6


    #bmap.etopo(ax=axins)

    bm_zoom = plot_changes_in_seasonal_means.get_arctic_basemap_nps(round=False, resolution="l")
    lower_left_xy = bmap(*lower_left_lat_lon)
    upper_right_xy = bmap(-160, 55)

    bm_zoom.etopo(ax=axins)
    #axins.set_xlim(lower_left_xy[0], lower_left_xy[0] + 400000)
    #axins.set_ylim(lower_left_xy[1], lower_left_xy[1] + 5000000)

    print(lower_left_xy)
    print(upper_right_xy)


    #bm_zoom.etopo(ax=axins)
    assert isinstance(axins, Axes)

    img = bm_zoom.contourf(x, y, data, norm=LogNorm(), alpha=0.8, ax=axins)
    bm_zoom.drawcoastlines(linewidth=0.3)
    bm_zoom.readshapefile("data/shp/wri_basins2/wribasin", "basin", color="k", linewidth=2, ax=axins)


    axins.set_xlim(lower_left_xy[0], upper_right_xy[0])
    axins.set_ylim(lower_left_xy[1], upper_right_xy[1])




    # draw a bbox of the region of the inset axes in the parent axes and
    # connecting lines between the bbox and the inset axes area
    mark_inset(ax, axins, loc1=1, loc2=4, fc="none", ec="0.5")
    fig.tight_layout()
    fig.savefig("simple2.png", bbox_inches="tight")

    plt.show()

def test_gettightbbox():
    fig, ax = plt.subplots(figsize=(8, 6))

    l, = ax.plot([1, 2, 3], [0, 1, 0])

    ax_zoom = zoomed_inset_axes(ax, 4)
    ax_zoom.plot([1, 2, 3], [0, 1, 0])

    mark_inset(ax, ax_zoom, loc1=1, loc2=3, fc="none", ec='0.3')

    remove_ticks_and_titles(fig)
    bbox = fig.get_tightbbox(fig.canvas.get_renderer())
    np.testing.assert_array_almost_equal(bbox.extents,
                                         [-17.7, -13.9, 7.2, 5.4])

def makeplot(refl, winds, w, stride, map, gs, title, file_name, box=None):
    pylab.figure()
    axmain = pylab.axes((0, 0.025, 1, 0.9))

    gs_x, gs_y = gs
    nx, ny = refl.shape
    xs, ys = np.meshgrid(gs_x * np.arange(nx), gs_y * np.arange(ny))

    pylab.contourf(xs, ys, refl, levels=np.arange(10, 80, 10))
    pylab.colorbar()
    pylab.contour(xs, ys, w, levels=np.arange(-10, 0, 2), colors='#666666', style='--')
    pylab.contour(xs, ys, w, levels=np.arange(2, 12, 2), colors='#666666', style='-')
    
    u, v = winds
    wind_slice = tuple([ slice(None, None, stride) ] * 2)
    pylab.quiver(xs[wind_slice], ys[wind_slice], u[wind_slice], v[wind_slice])

    if box:
        lb_y, lb_x = [ b.start for b in box ]
        ub_y, ub_x = [ b.stop for b in box ]

        box_xs = gs_x * np.array([ lb_x, lb_x, ub_x, ub_x, lb_x])
        box_ys = gs_y * np.array([ lb_y, ub_y, ub_y, lb_y, lb_y])

        map.plot(box_xs, box_ys, '#660099')

        axins = zoomed_inset_axes(pylab.gca(), 4, loc=4)
        pylab.sca(axins)

        pylab.contourf(xs[box], ys[box], refl[box], levels=np.arange(10, 80, 10))
        pylab.contour(xs[box], ys[box], w[box], levels=np.arange(-10, 0, 2), colors='#666666', style='--')
        pylab.contour(xs[box], ys[box], w[box], levels=np.arange(2, 12, 2), colors='#666666', style='-')
        pylab.quiver(xs[box], ys[box], u[box], v[box])

        drawPolitical(map)

        pylab.xlim([lb_x * gs_x, ub_x * gs_x - 1])
        pylab.ylim([lb_y * gs_y, ub_y * gs_y - 1])

        mark_inset(axmain, axins, loc1=1, loc2=3, fc='none', ec='k')

    pylab.sca(axmain)
    drawPolitical(map)

    pylab.suptitle(title)
    pylab.savefig(file_name)
    pylab.close()
    return

def test_inset_axes():
    def get_demo_image():
        from matplotlib.cbook import get_sample_data
        import numpy as np
        f = get_sample_data("axes_grid/bivariate_normal.npy", asfileobj=False)
        z = np.load(f)
        # z is a numpy array of 15x15
        return z, (-3, 4, -4, 3)

    fig, ax = plt.subplots(figsize=[5, 4])

    # prepare the demo image
    Z, extent = get_demo_image()
    Z2 = np.zeros([150, 150], dtype="d")
    ny, nx = Z.shape
    Z2[30:30 + ny, 30:30 + nx] = Z

    # extent = [-3, 4, -4, 3]
    ax.imshow(Z2, extent=extent, interpolation="nearest",
              origin="lower")

    # creating our inset axes with a bbox_transform parameter
    axins = inset_axes(ax, width=1., height=1., bbox_to_anchor=(1, 1),
                       bbox_transform=ax.transAxes)

    axins.imshow(Z2, extent=extent, interpolation="nearest",
                 origin="lower")
    axins.yaxis.get_major_locator().set_params(nbins=7)
    axins.xaxis.get_major_locator().set_params(nbins=7)
    # sub region of the original image
    x1, x2, y1, y2 = -1.5, -0.9, -2.5, -1.9
    axins.set_xlim(x1, x2)
    axins.set_ylim(y1, y2)

    plt.xticks(visible=False)
    plt.yticks(visible=False)

    # draw a bbox of the region of the inset axes in the parent axes and
    # connecting lines between the bbox and the inset axes area
    mark_inset(ax, axins, loc1=2, loc2=4, fc="none", ec="0.5")

    asb = AnchoredSizeBar(ax.transData,
                          0.5,
                          '0.5',
                          loc='lower center',
                          pad=0.1, borderpad=0.5, sep=5,
                          frameon=False)
    ax.add_artist(asb)

def plotConfusion(confusion, grid, title, file_name, inset=None, fudge=16):
    pylab.figure()
    axmain = pylab.axes()

    tick_labels = [ "Missing", "Correct\nNegative", "False\nAlarm", "Miss", "Hit" ]
    min_label = -1

    xs, ys = grid.getXY()

    pylab.pcolormesh(xs, ys, confusion, cmap=confusion_cmap, vmin=min_label, vmax=(min_label + len(tick_labels) - 1))

    tick_locs = np.linspace(-1, min_label + len(tick_labels) - 2, len(tick_labels))
    tick_locs += (tick_locs[1] - tick_locs[0]) / 2
    bar = pylab.colorbar()
    bar.locator = FixedLocator(tick_locs)
    bar.formatter = FixedFormatter(tick_labels)
    pylab.setp(pylab.getp(bar.ax, 'ymajorticklabels'), fontsize='large')
    bar.update_ticks()

    grid.drawPolitical()

    if inset:
        lb_y, lb_x = [ b.start for b in inset ]
        ub_y, ub_x = [ b.stop + fudge for b in inset ]

        inset_exp = (slice(lb_y, ub_y), slice(lb_x, ub_x))

        axins = zoomed_inset_axes(pylab.gca(), 2, loc=4)
        pylab.sca(axins)

        pylab.pcolormesh(xs[inset_exp], ys[inset_exp], confusion[inset_exp], cmap=confusion_cmap, vmin=min_label, vmax=(min_label + len(tick_labels) - 1))
        grid.drawPolitical()

        gs_x, gs_y = grid.getGridSpacing()

        pylab.xlim([lb_x * gs_x, (ub_x - 1) * gs_x])
        pylab.ylim([lb_y * gs_y, (ub_y - 1) * gs_y])

        mark_inset(axmain, axins, loc1=1, loc2=3, fc='none', ec='k')

    pylab.sca(axmain)
    pylab.suptitle(title)
    pylab.savefig(file_name)
    pylab.close()
    return

def draw_inset(plt, m, pos, lat, lon):
    ax = plt.subplot(111)
    zoom = 50
    axins = zoomed_inset_axes(ax, zoom, loc=1)
    m.plot(lon, lat, '.b--', zorder=10, latlon=True)
    m.scatter(lon, lat,  # longitude first!
              latlon=True,  # lat and long in degrees
              zorder=11)  # on top of all
    x1, y1 = m(lon[1] - 0.005, lat[0] - 0.0025)
    x2, y2 = m(lon[1] + 0.005, lat[0] + 0.0025)
    axins.set_xlim(x1, x2)
    axins.set_ylim(y1, y2)

    plt.xticks(visible=False)
    plt.yticks(visible=False)
    # draw a bbox of the region of the inset axes in the parent axes and
    # connecting lines between the bbox and the inset axes area
    mark_inset(ax, axins, loc1=2, loc2=4, fc="none", ec="0.5")

def plot_site(options):
    options['prefix'] = 'site'
    fig = MyFig(options, figsize=(10, 8), legend=True, grid=False, xlabel=r'Probability $p_s$', ylabel=r'Reachability~$\reachability$', aspect='auto')
    fig_vs = MyFig(options, figsize=(10, 8), legend=True, grid=False, xlabel=r'Fraction of Forwarded Packets~$\forwarded$', ylabel=r'Reachability~$\reachability$', aspect='auto')

    if options['grayscale']:
        colors = options['graycm'](pylab.linspace(0, 1.0, 3))
    else:
        colors = fu_colormap()(pylab.linspace(0, 1.0, 3))
    nodes = 105

    axins = None
    if options['inset_loc'] >= 0:
        axins = zoomed_inset_axes(fig_vs.ax, options['inset_zoom'], loc=options['inset_loc'])
        axins.set_xlim(options['inset_xlim'])
        axins.set_ylim(options['inset_ylim'])
        axins.set_xticklabels([])
        axins.set_yticklabels([])
        mark_inset(fig_vs.ax, axins, loc1=2, loc2=3, fc="none", ec="0.5")

    for i, (name, label) in enumerate(names):
        data = parse_site_only(options['datapath'][0], name)
        rs = numpy.array([r for p, r, std, conf, n, fw, fw_std, fw_conf in data if r <= options['limit']])
        fws = numpy.array([fw for p, r, std, conf, n, fw, fw_std, fw_conf in data if r <= options['limit']])/(nodes-1)
        ps = numpy.array([p for p, r, std, conf, n, fw, fw_std, fw_conf in data]) #[0:len(rs)])
        yerr = numpy.array([conf for p,r,std,conf, n, fw, fw_std, fw_conf in data]) #[0:len(rs)])

        patch_collection = PatchCollection([conf2poly(ps, list(rs+yerr), list(rs-yerr), color=colors[i])], match_original=True)
        patch_collection.set_alpha(0.3)
        patch_collection.set_linestyle('dashed')
        fig.ax.add_collection(patch_collection)

        fig.ax.plot(ps, rs, label=label, color=colors[i])
        fig_vs.ax.plot(fws, rs, label=label, color=colors[i])
        if axins:
            axins.plot(fws, rs, color=colors[i])

    fig.ax.set_ylim(0, options['limit'])
    fig.legend_title = 'Graph'
    fig.save('graphs')
    fig_vs.legend_title = 'Graph'
    fig_vs.ax.set_xlim(0,1)
    fig_vs.ax.set_ylim(0,1)
    fig_vs.save('vs_pb')

def plot_calibration(c, insert = True):
    fig = plt.figure()
    ax = plt.gca()
    for baseline, correlation in c.time_domain_correlations_values.items():
        correlation_max_val = correlation[np.argmax(correlation)]
        correlation_max_time = c.time_domain_correlations_times[baseline][np.argmax(correlation)]
        lines = ax.plot(
            c.time_domain_correlations_times[baseline] * 1e9,
            correlation / correlation_max_val,
            label = baseline)
        #ax.plot([correlation_max_time, correlation_max_time], [0, correlation_max_val], color = lines[0].get_color())
    ax.set_ylim(top=1.2)
    ax.xaxis.set_ticks(np.arange(
        -200,
        200,
        2))
    if insert == True:
        #axins = zoomed_inset_axes(ax, 5, loc=1)
        axins = zoomed_inset_axes(ax, 9, loc=1)
        for baseline, correlation in c.time_domain_correlations_values.items():
            correlation_max_val = correlation[np.argmax(correlation)]
            correlation_max_time = c.time_domain_correlations_times[baseline][np.argmax(correlation)]
            lines = axins.plot(
                c.time_domain_correlations_times[baseline] * 1e9,
                correlation / correlation_max_val,
                label = baseline,
                linewidth=2)
            #axins.plot([correlation_max_time, correlation_max_time], [0, correlation_max_val], color = lines[0].get_color())
        #axins.set_xlim(-0.4, 2.9)
        axins.set_xlim(-0.4, 0.4)
        #axins.set_ylim(0.90, 1.04)
        axins.set_ylim(0.96, 1.03)
        #axins.xaxis.set_ticks(np.arange(-0.4, 2.9, 0.4))
        axins.xaxis.set_ticks(np.arange(-0.4, 0.4, 0.2))
        mark_inset(ax, axins, loc1=2, loc2=3, fc='none', ec='0.5')

    plt.xticks(visible=True)
    plt.yticks(visible=False)
    ax.set_title("Time domain cross correlations with broad band noise\n arriving through full RF chain AFTER calibration")
    ax.set_xlabel("Time delay (ns)")
    ax.set_ylabel("Cross correlation value (normalised)")
    ax.legend(loc=2)
    #ax.legend()
    plt.show()

 def plot_figure(self, ax, xaxis, yaxis, title='', legend='', 
         drawline=True, legend_outside=False, marker=None, 
         linestyle=None, xlabel='', ylabel='', xlog=False, ylog=False, 
         zoom=False, zoom_ax=None, zoom_xaxis=[], zoom_yaxis=[], 
         legend_pos='best', xlabel_format=0, xlimit=0, ylimit=0, legend_markscale=1.0):
     if drawline:
         ax.plot(xaxis, yaxis, marker=marker if marker else '+', ls=linestyle if linestyle else '-', label=legend)
     else: #draw dots 
         ax.plot(xaxis, yaxis, marker=marker if marker else 'o', markerfacecolor='r', ms=4, ls='None', label=legend)
         #ax.vlines(xaxis, [0], yaxis)
     if xlog:
         ax.set_xscale('log')
     if ylog:
         ax.set_yscale('log')
     if xlimit > 0:
         ax.set_xlim(0, ax.get_xlim()[1] if ax.get_xlim()[1]<xlimit else xlimit)
     if ylimit > 0:
         ax.set_ylim(0, ax.get_ylim()[1] if ax.get_ylim()[1]<ylimit else ylimit)
     ax.set_title(title)
     ax.legend(loc=legend_pos, markerscale=legend_markscale)
     ax.set_xlabel(xlabel)
     ax.set_ylabel(ylabel)
     ax.ticklabel_format(axis='y', style='sci', scilimits=(0,0))
     # zoom
     if zoom:
         new_zoom_ax = False
         if zoom_ax is None:
             new_zoom_ax = True  
             zoom_ax = inset_axes(ax,
                    width="70%",  # width = 50% of parent_bbox
                    height="70%",  # height : 1 inch
                    loc=7) # center right
         if drawline:
             zoom_ax.plot(zoom_xaxis, zoom_yaxis, marker=marker if marker else '+', ls=linestyle if linestyle else '-')
         else: #draw dots 
             zoom_ax.plot(zoom_xaxis, zoom_yaxis, marker=marker if marker else 'o', markerfacecolor='r', ms=4, ls='None')
             #zoom_ax.vlines(zoom_xaxis, [0], zoom_yaxis)
         if new_zoom_ax:
             zoom_ax.ticklabel_format(axis='y', style='sci', scilimits=(0,0))
             mark_inset(ax, zoom_ax, loc1=2, loc2=4, fc="none", ec="0.5")
     return ax, zoom_ax

def graph(data_, j_size, title=None):
    fig, ax = plt.subplots()
    data = pd.read_csv(data_, index_col = 0)
    jitter_ = jitter(-j_size, j_size, len(list(data.iterrows())[0][1]))
    max_1 = 0
    max_2 = 0
    for i, d in enumerate(data.iterrows()):
        if max(list(d[1])) > max_1 and i > 2:
            max_1 = max(list(d[1]))
        elif max(list(d[1])) > max_2 and i < 3:
            max_2 = max(list(d[1]))
        plt.scatter(x = np.array([i]*len(jitter_)+jitter_), y = list(d[1]))
    plt.xlim(-0.1, 5.1) # apply the x-limits
    c1 = 1
    while max_1/(c1*10) > 1:
        c1 *= 10
    c2 = 1
    while max_2/(c2*10) > 1:
        c2 *= 10
    cb = max(c1,c2)
    cs = min(c1,c2)
    max1 = max(max_1, max_2)
    max2 = min(max_1, max_2)
    plt.ylim(-cb, (max1/cb+2)*cb)
    from mpl_toolkits.axes_grid1.inset_locator import inset_axes
    if c2 > c1:
        axins = inset_axes(ax, 4,3, loc=5)
    else:
        axins = inset_axes(ax, 4,3, loc=6)        
        axins.yaxis.tick_right()
    for i, d in enumerate(data.iterrows()):
        plt.scatter(x = np.array([i]*len(jitter_)+jitter_), y = list(d[1]))
    if c2 > c1:
        axins.set_xlim(2.9, 5.1) # apply the x-limits
    else:
        axins.set_xlim(-0.1, 2.1)
    plt.ylim(-cs, (max2/cs+2)*cs)
    from mpl_toolkits.axes_grid1.inset_locator import mark_inset
    mark_inset(ax, axins, loc1=2, loc2=4, fc="none", ec="0.5")
    plt.show()        

def nice_plot_e_of_T():
    Ts1 = np.arange(273.15-10, 273.15-5, 1.)
    Ts2 = np.arange(273.15-3, 273.15-1, .1)
    Ts3 = np.arange(273.15-1, 273.15+1, .0001)
    Ts4 = np.arange(273.15+1, 273.15+5, .1)
    Ts5 = np.arange(273.15+5, 273.15+50, 1.)

    Ts = np.concatenate((Ts1, Ts2, Ts3, Ts4, Ts5))

    # water content at T=273.65, p=85000.
    wc = water_content(273.65, 85000., saturations_A)
    ie = []

    for T in Ts:
        # for each T, determine the p required to maintain same mass
        (sl,si,sg),p = saturation_with_wc(T,wc, saturations_A)
        ie.append(internal_energy_per_mol(T,p, saturations_A))

    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.plot(Ts - 273.15,np.array(ie), 'b')
    ax.set_xlim(-10,20)
    ax.set_ylim(-8000,2000)
    plt.xlabel("temperature [C]")
    plt.ylabel("energy density (water) [J/mol]")

    axins = inset_axes(ax, 3, 3, loc=4)
    axins.plot(Ts - 273.15,np.array(ie), 'b')
    axins.set_xlim(-.01,.01)
    axins.set_ylim(-7000,500)
    axins.set_xticks([-0.01, 0, 0.01])
    axins.tick_params(labeltop=1, labelbottom=0)
    axins.set_yticks([-6000,-4000,-2000,0])

    mark_inset(ax, axins, loc1=1, loc2=3, fc="none", ec="0.5")
    plt.show()