def plot_Bn_slice(Bn_List, q95_array, Bn_Div_Li_array, coil1_abs_array, coil1_angle_array, ROTE_value, probe):
    fig = pt.figure()
    ax = fig.add_subplot(211)
    ax2 = fig.add_subplot(212)
    color_list = ['bo-','ko-','yo-']
    color_list2 = ['bx-','kx-','yx-']
    for iii in range(0,len(Bn_List)):
        Bn_Div_Li_value = Bn_List[iii]
        q95_range = num.arange(2.5,6.6,0.1)
        Bn_Div_Li = num.ones(len(q95_range),dtype=float)*Bn_Div_Li_value
        interp_abs_line = griddata(q95_array, Bn_Div_Li_array, coil1_abs_array, q95_range,Bn_Div_Li, interp = 'linear')[0]
        interp_angle_line = griddata(q95_array, Bn_Div_Li_array, coil1_angle_array, q95_range, Bn_Div_Li, interp = 'linear')[0]

        ax.plot(q95_range, interp_abs_line,color_list[iii],label='Bn/Li='+str(Bn_Div_Li_value))
        ax2.plot(q95_range, interp_angle_line,color_list[iii],label='Bn/Li='+str(Bn_Div_Li_value))

    leg = ax.legend(loc=2, fancybox = True)
    leg.get_frame().set_alpha(0.5)
    leg = ax2.legend(loc=4, fancybox = True)
    leg.get_frame().set_alpha(0.5)
    ax.grid()
    ax2.grid()
    ax2.set_ylim([0,360])
    ax.set_xlim([2,7])
    ax2.set_xlim([2,7])
    ax2.set_xlabel('q95')
    ax.set_title(str(ROTE_value) + 'deg ' + probe[iii])
    ax.set_ylabel('abs(output)')
    ax.set_ylabel('Magnitude (G/kA)')
    ax2.set_ylabel('Phase (deg)')
    fig.canvas.draw()
    fig.show()

def plot_2d(xs,ys,zs,ax,zmax=None,logscale=False,cb=True,zmin_avg=True,shrink=0.7) :
    if zmax == None :
        zmax = max(zs)
    minx,maxx = min(xs),max(xs)
    miny,maxy = min(ys),max(ys)
    xi = np.linspace(minx,maxx,360)
    yi = np.linspace(miny,maxy,360)
    if zmax == 100 or not logscale :
        zi = griddata(xs,ys,zs,xi,yi)
    else :
        Z = []
        zmin = 1e9
        for z in zs :
            try :
                Z.append(log10(z))
                if log10(z) < zmin :
                    zmin = int(np.floor(log10(z)))
            except ValueError :
                Z.append(-20)
            except :
                print "Unknown error trying to take the log of %s"%z
                sys.exit()
        zi = griddata(xs,ys,Z,xi,yi)
    myplot = ax.pcolorfast(xi,yi,zi,cmap='RdBu')
    x0,x1 = ax.get_xlim()
    y0,y1 = ax.get_ylim()
    aspect_Ratio = (x1-x0)/(y1-y0)
    ax.set_aspect(aspect_Ratio)
    ax.set_xlim(-180,175)
    ax.set_ylim(-180,175)
    ax.set_xticks(range(-120,121,60))
    ax.set_yticks(range(-120,121,60))

    if zmax == 100 and not logscale :
        myplot.set_clim([0,100])
        if cb :
            cbar = plt.colorbar(myplot,shrink=shrink,format='%i',ticks=range(0,101,20))
            cbar.ax.set_yticklabels(['0%','20%','40%','60%','80%','>100%'])
        return ax, cbar
    elif logscale :
        # Set the minimum value to the average magnitude
        if zmin_avg :
            zmin=int(np.floor(log10(np.average(zs))))-1
        myplot.set_clim(zmin,-1)
        if cb :
            cbar = plt.colorbar(myplot,shrink=shrink,format='%i',ticks=range(zmin,0,1))
            cblabels = []
            for i in range(zmin,0,1) :
                if i == zmin :
                    cblabels.append('<1E%i'%i)
                else :
                    cblabels.append('1E%i'%i)
            cbar.ax.set_yticklabels(cblabels)
            return ax,cbar
    else :
        if cb :
            cbar = plt.colorbar(myplot,shrink=shrink,format='%1.E')
            return ax,cbar
    return ax

    def interpolateImagePoints(self, points, edges):
        t_row, t_col = np.ogrid[0:self._image.shape[0], 0:self._image.shape[1]]

        # get values for 
        xCoords = [x[0] for x in points]
        yCoords = [x[1] for x in points]

        xValues = [self.getXAt(*p) for p in points]
        yValues = [self.getYAt(*p) for p in points]
        zValues = [self.getZAt(*p) for p in points]

        # get edge points
        edgeIndices = [(points.index(e[0]),
                        points.index(e[1])) for e in edges]
        for p1, p2 in edgeIndices:
            p1v = np.array(points[p1])
            p2v = np.array(points[p2])
            mag = p2v-p1v
            mag = int(np.sqrt(mag.dot(mag)))

            interpU = np.interp(range(mag),
                                [0, mag-1],
                                [p1v[0], p2v[0]]).astype(np.int).tolist()
            interpV = np.interp(range(mag),
                                [0, mag-1],
                                [p1v[1], p2v[1]]).astype(np.int).tolist()
            interpX = np.interp(range(mag),
                                [0, mag-1],
                                [xValues[p1], xValues[p2]]).tolist()
            interpY = np.interp(range(mag),
                                [0, mag-1],
                                [yValues[p1], yValues[p2]]).tolist()
            interpZ = np.interp(range(mag),
                                [0, mag-1],
                                [zValues[p1], zValues[p2]]).tolist()
            xCoords.extend(interpU)
            yCoords.extend(interpV)
            xValues.extend(interpX)
            yValues.extend(interpY)
            zValues.extend(interpZ)

        # interpolate points
        interpX = griddata(xCoords, yCoords, xValues,
                           t_col.ravel(), t_row.ravel(),
                           "linear")
        self._image_points[t_row, t_col, 0] = interpX

        interpY = griddata(xCoords, yCoords, yValues,
                           t_col.ravel(), t_row.ravel(),
                           "linear")
        self._image_points[t_row, t_col, 1] = interpY

        interpZ = griddata(xCoords, yCoords, zValues,
                           t_col.ravel(), t_row.ravel(),
                           "linear")
        self._image_points[t_row, t_col, 2] = interpZ

def pcolor(
    lon,
    lat,
    lonn,
    latn,
    mag,
    nv,
    m,
    ax,
    norm,
    cmap,
    topology_type,
    fig,
    height,
    width,
    lonmin,
    latmin,
    lonmax,
    latmax,
    dataset,
    continuous,
    projection,
):
    from matplotlib.mlab import griddata

    if topology_type.lower() == "cell":
        lon, lat = m(lon, lat)
        lonn, latn = m(lonn, latn)
    else:
        lon, lat = m(lonn, latn)
        lonn, latn = lon, lat
    num = int((lonmax - lonmin) * 320)
    xi = np.arange(m.xmin, m.xmax, num)
    yi = np.arange(m.ymin, m.ymax, num)
    if topology_type.lower() == "node":
        n = np.unique(nv)
        zi = griddata(lon[n], lat[n], mag[n], xi, yi, interp="nn")
    else:
        zi = griddata(lon, lat, mag, xi, yi, interp="nn")
    fig, m, patch1 = cookie_cutter(dataset, fig, m, lonmin, latmin, lonmax, latmax, projection, continuous)

    # Should we draw anything?
    if patch1 is not None:
        m.imshow(zi, norm=norm, cmap=cmap, clip_path=patch1, interpolation="nearest")

    # from matplotlib.backends.backend_agg import FigureCanvasAgg
    # canvas = FigureCanvasAgg(fig)
    # canvas.print_png("testing_yay.png")
    return fig, m

def plot_current_field(xs,ys,vx,vy):
    ###define grid
    x_grid = np.linspace(xs.min(),xs.max(),200)
    y_grid = np.linspace(ys.min(),ys.max(),200)

    vxs = griddata(xs, ys, vx, x_grid, y_grid, interp='linear')
    vys = griddata(xs, ys, vy, x_grid, y_grid, interp='linear')


    speed = np.sqrt(vxs**2+vys**2)
    lw = 5*speed / speed.max()
    plt.figure()
    plt.streamplot(x_grid, y_grid, vxs,vys, color=lw, linewidth=2, cmap=plt.cm.autumn)
    plt.colorbar()
    plt.show()

def return_grid_data(q95_values, Bn_Div_Li_values,q95_array, Bn_Div_Li_array, coil1_angle_array, coil1_abs_array, xnew=None,ynew=None, interpolation='linear',deg_min=0, points = [100,100]):
    if xnew==None:
        xnew = num.linspace(q95_values[0], q95_values[1], points[0])
        ynew = num.linspace(Bn_Div_Li_values[0], Bn_Div_Li_values[1],points[1])

    for i in range(0,len(coil1_angle_array)):
        while coil1_angle_array[i]<deg_min or coil1_angle_array[i]>(deg_min+360):
            if coil1_angle_array[i]<deg_min:
                coil1_angle_array[i] += 360.
            if coil1_angle_array[i]>(deg_min+360):
                coil1_angle_array[i] -= 360.

    B1grid_data = griddata(q95_array, Bn_Div_Li_array, coil1_abs_array, xnew, ynew, interp = interpolation)
    interp_data_angle = griddata(q95_array, Bn_Div_Li_array, coil1_angle_array, xnew, ynew, interp = interpolation)
    return B1grid_data, interp_data_angle

def griddata_all(x,y,z,x1,y1,func='line_rbf'):
    '''
    把各种插值方法整合到一起
    scipy_idw
    line_rbf
    Invdisttree
    nat_grid
    '''

    xi, yi = np.meshgrid(x1, y1)


    if('nearest'==func):
        zi= griddata_nearest(x,y,z,xi,yi)

    if('griddata'==func):
        from matplotlib.mlab import griddata
        zi = griddata(x,y,z,x1,y1)

    if('kriging'==func):
        zi= griddata_kriging(x,y,z,xi,yi)

    if('scipy_idw'==func):
        zi= griddata_scipy_idw(x,y,z,xi,yi)

    if('line_rbf'==func):
        zi = griddata_linear_rbf(x,y,z,xi,yi)  #        grid3 = grid3.reshape((ny, nx))
        print(zi.shape,x.shape,y.shape,z.shape,xi.shape,yi.shape)
        #sys.exit(0)

    if('line_rbf2'==func):
        zi = griddata_linear_rbf2(x,y,z,xi,yi)  #        grid3 = grid3.reshape((ny, nx))


    if('Invdisttree'==func):
        #zi = df.griddata_Invdisttree(x,y,z,xi,yi,Nnear=15,p=3,eps=1)
        print(x.shape,y.shape,z.shape,x1.shape,y1.shape)
        #sys.exit(0)
        zi = griddata_Invdisttree(x,y,z,xi,yi,p=3)#,Nnear=10,eps=1)

    if('nat_grid'==func):
        from griddata import griddata, __version__
        zi = griddata(x,y,z,xi,yi)

    #if('test'==func):
    #    zi = griddata_scipy_spatial(x,y,z,xi,yi)

    return zi,xi,yi

def gridxyz(xcol, ycol, zcol, xystep=None, lib='scipy', method='cubic'):
    """ grid (X, Y, Z) 1D data on a 2D regular mesh
    
    Parameters
    ----------
    xcol, ycol, zcol : 1D arrays repesenting the map (z is the intensity)
    xystep : the step size of the XY grid
    lib : library used for griddata
          [scipy]
          matplotlib
    method : interpolation method
    
    Returns
    -------
    xgrid, ygrid : 1D arrays giving abscissa and ordinate of the map
    zz : 2D array with the gridded intensity map
    
    See also
    --------
    - MultipleScanToMeshPlugin in PyMca
    """
    if xystep is None:
        xystep = 0.05
        warnings.warn("'xystep' not given: using a default value of {0}".format(xystep))
    #create the XY meshgrid and interpolate the Z on the grid
    xgrid = np.linspace(xcol.min(), xcol.max(), (xcol.max()-xcol.min())/xystep)
    ygrid = np.linspace(ycol.min(), ycol.max(), (ycol.max()-ycol.min())/xystep)
    xx, yy = np.meshgrid(xgrid, ygrid)
    if ('matplotlib' in lib.lower()):
        try:
            from matplotlib.mlab import griddata
        except ImportError:
            print("Error: cannot load griddata from Matplotlib")
            return
        if not (method == 'nn' or method == 'nearest'):
            warnings.warn("method {0} not supported by {1}".format(method, lib))
        print("Gridding data with {0}...".format(lib))
        zz = griddata(xcol, ycol, zcol, xx, yy)
        return xgrid, ygrid, zz
    elif ('scipy' in lib.lower()):
        try:
            from scipy.interpolate import griddata
        except ImportError:
            print("Error: cannot load griddata from Scipy")
            return
        print("Gridding data with {0}...".format(lib))
        zz = griddata((xcol, ycol), zcol, (xgrid[None,:], ygrid[:,None]), method=method)
        return xgrid, ygrid, zz

 def grid (self,lgrid,mgrid,freqgrid):
   from matplotlib.mlab import griddata
   l = numpy.sin(lgrid);
   m = numpy.sin(mgrid);
   gridshape = (len(l),len(m),len(self.freq));
   beamgrid = numpy.zeros(gridshape,complex);
   # figure out useful range of l/m coordinates
   lmin,lmax = l.min(),l.max();
   mmin,mmax = m.min(),m.max();
   dl = (lmax - lmin)/10;
   dm = (mmax - mmin)/10;
   lmin,lmax = lmin-dl,lmax+dl;
   mmin,mmax = mmin-dm,mmax+dm;
   dprint(3,"regridding in area from %f,%f to %f,%f deg"%(lmin/DEG,mmin/DEG,lmax/DEG,mmax/DEG));
   # loop over all per-frequency patterns
   for ifreq,(beam,beam_ampl,lcoord,mcoord,freq) in enumerate(self._beams):
     # filter out values outside the range
     select = (lcoord < lmax)&(lcoord > lmin)&(mcoord < mmax)&(mcoord > mmin); 
     dprint(3,"selection reduces beam from %d to %d points"%(len(lcoord),select.sum()));
     # regrid onto each frequency plane
     for ifreq in range(len(self.freq)):
       gbeam = beamgrid[:,:,ifreq];
       gbeam.real = griddata(lcoord[select],mcoord[select],beam.real[select],l,m);
       gbeam.imag = griddata(lcoord[select],mcoord[select],beam.imag[select],l,m);
       if beam_ampl is not None:
         ampl = griddata(lcoord[select],mcoord[select],beam_ampl[select],l,m);
         ampl1 = abs(gbeam);
         wh = ampl1 != 0;
         gbeam[wh] = (gbeam[wh]/ampl1[wh])*ampl[wh];
     # replace nans with zeroes
     beamgrid[numpy.isnan(beamgrid)] = 0;
   # interpolate onto frequency grid
   dprint(4,"interpolated data range",beamgrid.min(),beamgrid.max());
   dprint(3,"interpolating frequencies");
   l1,m1 = numpy.meshgrid(numpy.arange(len(l)),numpy.arange(len(m)));
   if len(self.freq) > 1:
     f1 = interpolate.interp1d(self.freq,range(len(self.freq)),'linear')(freqgrid);
   else:
     f1 = numpy.array(0.);
   l1,f1 = unite_shapes(l1,f1);
   m1,f1 = unite_shapes(m1,f1);
   coords = numpy.vstack((l1.ravel(),m1.ravel(),f1.ravel()));
   output_shape = (len(l),len(m),len(freqgrid));
   output = numpy.zeros(output_shape,complex);
   output.real = interpolation.map_coordinates(beamgrid.real,coords,order=self._spline_order).reshape(output_shape);
   output.imag = interpolation.map_coordinates(beamgrid.imag,coords,order=self._spline_order).reshape(output_shape);
   dprint(3,"done");
   return output;

 def polar(self, n, file='None'):
     r = np.arange(self.rmin,self.rmax,(self.rmax-self.rmin)/self.nrad)
     t = np.arange(0.,2.*np.pi,2.*np.pi/self.nsec)
     x = np.ndarray([self.nrad*self.nsec], dtype = float)
     y = np.ndarray([self.nrad*self.nsec], dtype = float)
     z = np.ndarray([self.nrad*self.nsec], dtype = float)
     k = 0
     for i in range(self.nrad):
         for j in range(self.nsec):
             x[k] = r[i]*np.cos(t[j])
             y[k] = r[i]*np.sin(t[j])
             z[k] = self.data[i,j]
             k +=1
     xx = np.arange(-self.rmax, self.rmax, (self.rmax-self.rmin)/n)
     yy = np.arange(-self.rmax, self.rmax, (self.rmax-self.rmin)/n)
     zz = griddata(x,y,z,xx,yy)
     fig_pol = plt.figure()
     ax1 = fig_pol.add_subplot(111, axisbg='k')
     ax1.set_xlabel("X")
     ax1.set_ylabel("Y")
     if(self.zmax!='None' and self.zmin!='None'):
         ax1.imshow(zz, cmap=cm.hot, origin="lower", \
                        extent=[-self.rmax,self.rmax, \
                                     -self.rmax,self.rmax])
     else:
         ax1.imshow(zz, cmap=cm.hot, origin="lower", \
                        extent=[-self.rmax,self.rmax, \
                                     -self.rmax,self.rmax])
     if(file!="None"):
         plt.savefig(file+".png",dpi=70, format="png" )
         print file+".png done"
     else:
         plt.show()

def plot_horizon(xyz_data, xyz_titles, az_lims, el_lims):
    """Calculate and plot horizon.

    Parameters
    ----------
    xyz_data : tuple of (azimuth, elevation, power_db) data lists
        Data to plot.
    xyz_titles : tuple of titles for (azimuth, elevation and power_db) axes
        Titles for axes.
    az_lims : tuple of (min azimuth, max azimuth)
        Azimuth limits for the plot.
    el_lims : tuple of (min elevation, max elevation)
        Elevation limits for the plot.
    pow_lims : tuple of (min power, max power)
        Power limits for the plot.
    """
    azimuth, elevation, power_db = xyz_data
    az_title, el_title, pow_title = xyz_titles
    az_min, az_max = az_lims
    el_min, el_max = el_lims
    #pow_min, pow_max = pow_lims
    az_pos = np.linspace(az_min, az_max, (az_max - az_min) / 0.1)
    el_pos = np.linspace(el_min, el_max, (el_max - el_min) / 0.1)
    power_db_pos = mlab.griddata(azimuth[:,0], elevation[:,0], power_db[:,0], az_pos, el_pos)
    plt.imshow(power_db_pos, aspect='auto', origin='lower')
    #cs = plt.contour(az_pos, el_pos, power_db_pos, pow_levels)
    #plt.contourf(az_pos,el_pos, power_db_pos, pow_levels, antialiased=True)
    plt.colorbar()

    if az_title:
        plt.xlabel(az_title)
    if el_title:
        plt.ylabel(el_title)
    if pow_title:
        plt.title(pow_title)

def ribbon3(spectra, alpha):
    matplotlib.rcParams.update({'font.size':10})
    fig=figure()
    ax=fig.gca(projection='3d')
    for i in range(0,spectra.shape[1]):
        y=spectra[:,i]
        x=sorted(range(1,len(y)+1)*2)
        a=[i,i+1]*len(y)
        b=list(itertools.chain(*zip(y,y)))
        xi=np.linspace(min(x),max(x))
        yi=np.linspace(min(a),max(a))
        X,Y=np.meshgrid(xi/(len(x)*0.5),yi)
        Z=griddata(x,a,b,xi,yi)
        ax.plot_surface(X,Y,Z, rstride=50, cstride=1, cmap='Spectral')
        ax.set_zlim3d(np.min(Z),np.max(Z))

    ax.grid(False)
    ax.w_xaxis.pane.set_visible(False)
    ax.w_yaxis.pane.set_visible(False)
    ax.w_zaxis.pane.set_color('gainsboro')
    ax.set_title('Pyramid Gradient MFS')
    ax.set_xlim3d(0,1)
    ax.set_xticks(alpha)
    ax.set_xticklabels(alpha)
    ax.set_xlabel(r'$\alpha$')
    #ax.set_yticks([0.5,1.5,2.5,3.5,4.5])
    #ax.set_yticklabels(['1','2','3','4','5'])
    ax.set_ylabel('Resolution')
    ax.set_zlim3d(0,3)
    ax.set_zlabel(r'$f(\alpha)$')
    show()

def regularlyGrid(xarr, yarr, zarr, xstart=None, xfinal=None, xstep=None, ystart=None, yfinal=None, ystep=None):
    "Returns the regularly grided xi, yi, and zi arrays from the initial data."
    # if xstart,xfinal,xstep,ystart,yfinal,ystep are NOT given, they are derived from the data
    if xstart == None:
        xstart = xarr.min() 
    if xfinal == None:
        xfinal = xarr.max()
    if xstep == None:
        xstep = 1.0 * (xfinal - xstart) / len(xarr)
    if ystart == None:
        ystart = yarr.min()
    if yfinal == None:
        yfinal = yarr.max()
    if ystep == None:
        ystep = 1.0 * (yfinal - ystart) / len(yarr)

    xi = N.arange(xstart, xfinal, xstep)
    yi = N.arange(ystart, yfinal, ystep)
    '''
	Programming note:
		linspace does (start, final, how many steps)
		arange does (start, final, step)
	'''

    # grid the data.
    print len(xarr), len(yarr), len(zarr), len(xi), len(yi)

    xarr,yarr,zarr=remove_duplicates(xarr,yarr,zarr)
    zi = griddata(xarr, yarr, zarr, xi, yi)

    print "done gridding"
    return xi, yi, zi

def window3(x,y,z):
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    x = x
    y = y
    z = z
    
    
    xi = np.linspace(1, 19)# axis values go here
    yi = np.linspace(1, 19)
    
    N, Y = np.meshgrid(xi, yi)
    Z = griddata(x, y, z, xi, yi)
    
    surf = ax.plot_surface(N, Y, Z, rstride=3, cstride=3, cmap=cm.jet,
                           linewidth=0, antialiased=True)
    
    ax.set_zlim3d(np.min(Z), np.max(Z))
    fig.colorbar(surf)
    def getMarker(i):
    # Use modulus in order not to have the index exceeding the lenght of the list (markers)
        return "$"+'\gamma'+"$"
    plt.xlabel('N', fontsize=18)
    plt.ylabel('Y', fontsize=18)
    plt.title("Intersection Error figure. \n View from top for heatmap")
    
    plt.show()

  def compareVars(self,histories):
    #find x,y in histories
    nameX=histories['vars'][0]
    nameY=histories['vars'][1]
    pathX=histories['varPaths'][0]
    pathY=histories['varPaths'][1]
    xs=np.zeros(histories['nRun'][-1]+1)
    ys=np.zeros(histories['nRun'][-1]+1)
    for e,edict in enumerate(histories['varVals']):
      xs[e]=edict[pathX]
      ys[e]=edict[pathY]
    zs=np.array(histories['soln'])

    xi=np.linspace(np.min(xs),np.max(xs),100)
    yi=np.linspace(np.min(ys),np.max(ys),100)

    plt.figure()
    zi=mlab.griddata(xs,ys,zs,xi,yi)
    plt.contour(xi,yi,zi,15,linewidth=0.5,colors='k')
    plt.pcolormesh(xi,yi,zi,cmap=plt.get_cmap('rainbow'))

    plt.colorbar()
    plt.scatter(xs,ys,marker='o',c='b',s=1,zorder=10)
    plt.xlim(np.min(xs),np.max(xs))
    plt.ylim(np.min(ys),np.max(ys))

    plt.xlabel(nameX)
    plt.ylabel(nameY)
    plt.title('Total Runs: '+str(len(xs)))

def printcontourn (name,outdir,surf,num):
    '''
    FUNCTION FOR PRINTING A CONTOURN OUT OF A FILE CONTAINING 
    X, Y AND POTENTIAL COLUMNS
    '''
    x = np.loadtxt(outdir+name+"_"+surf+".txt", usecols=[0], skiprows=1)
    y = np.loadtxt(outdir+name+"_"+surf+".txt", usecols=[1], skiprows=1)
    z = np.loadtxt(outdir+name+"_"+surf+".txt", usecols=[2], skiprows=1)

    dimx=int(x[:1])
    dimy=int(y[:1])
    uniques=len(np.unique(x))-1
    xi = np.linspace(dimx,-dimx,uniques)
    yi = np.linspace(dimy,-dimy,uniques*2)
    zi = griddata(x,y,z,xi,yi,interp='nn')

    #limits=[]
    #limits.append(int(min(z)))
    #limits.append(int(max(z)))
    plt.subplot(num)
    #normal=matcol.Normalize(vmin=min(limits),vmax=max(limits))
    normal=matcol.Normalize(vmin=-80,vmax=80)
    cmap=plt.cm.RdBu
    CS1 = plt.contour(xi,yi,zi,20,linewidths=0.5,colors='k')
    CS2 = plt.contourf(xi,yi,zi,20,cmap=cmap,norm=normal)
    cb=plt.colorbar(drawedges='TRUE') 
    #plt.scatter(x,y,marker='o',c='b',s=0.1,zorder=10)
    plt.xlim(dimx,-dimx)
    plt.ylim(dimy,-dimy)
    plt.title(surf, fontsize=11, fontweight='bold')

    def computeGrid(self):
        x, y, z = self._data

        self.gridSpacing = math.sqrt((max(x) - min(x)) *
                                     (max(y) - min(y)) / len(x)) / 2.0
        gridSpacing = self.gridSpacing

        if gridSpacing <= 0.0:
            raise ContourError("Grid spacing must be greater than 0")

        # make grid
        x0 = math.floor(min(x) / gridSpacing) * gridSpacing
        nx = int(math.floor((max(x) - x0) / gridSpacing)) + 1
        gx = np.linspace(x0, x0 + gridSpacing * nx, nx)

        y0 = math.floor(min(y) / gridSpacing) * gridSpacing
        ny = int(math.floor((max(y) - y0) / gridSpacing)) + 1
        gy = np.linspace(y0, y0 + gridSpacing * ny, ny)

        try:
            # interpolate values on grid
            gz = griddata(x, y, z, gx, gy)
        except:
            raise ContourError("Unable to generate a grid for this data set")

        self._gridData = (gx, gy, gz)

def grid(x, y, z, resX=100, resY=100):
  "Convert 3 column data to matplotlib grid"
  xi = np.linspace(min(x), max(x), resX)
  yi = np.linspace(min(y), max(y), resY)
  Z = griddata(x, y, z, xi, yi)
  X, Y = np.meshgrid(xi, yi)
  return X, Y, Z

def interp_2d_xy(x,y,data,xi,yi):
	try:
		di = griddata(x.reshape(x.size),y.reshape(y.size),data.reshape(data.size),xi,yi)
	except TypeError:
		di = np.zeros(xi.size)
		if x.ndim ==2 and y.ndim ==2:
			x_vec = x[0,:]
			y_vec = y[:,0]
		elif x.ndim == 3 and y.ndim == 3:
			x_vec = x[0,0,:]
			y_vec = y[0,:,0]
		else:
			x_vec = x
			y_vec = y
		xl = np.nonzero(x_vec <= xi)[0][-1]
		xh = np.nonzero(xi <= x_vec)[0][0]
		yl = np.nonzero(y_vec <= yi)[0][-1]
		yh = np.nonzero(yi <= y_vec)[0][0]
				
		if not x_vec[xl] == x_vec[xh]:
			xd = (xi-x_vec[xh])/(x_vec[xl]-x_vec[xh])
		else:
			xd = 1.
		if not y_vec[yl] == y_vec[yh]:
			yd = (yi-y_vec[yh])/(y_vec[yl]-y_vec[yh])
		else:
			yd = 1.

		w0 = data[yl,xl]*(1-yd) + data[yh,xl]*yd
		w1 = data[yl,xh]*(1-yd) + data[yh,xh]*yd
	
		di = w0*(1-xd) + w1*xd
	return di

def pcolorRandom():
    "Makes a pcolormesh plot of randomly generated data pts."
    # make up some randomly distributed data
    npts = 100
    x = uniform(-3, 3, npts)
    y = uniform(-3, 3, npts)
    z = x * N.exp(-x ** 2 - y ** 2)

    # define grid.
    xi = N.arange(-3.1, 3.1, 0.05)
    yi = N.arange(-3.1, 3.1, 0.05)

    # grid the data.
    zi = griddata(x, y, z, xi, yi)

    # contour the gridded data, plotting dots at the randomly spaced data points.
    plt.pcolormesh(xi, yi, zi)	
    #CS = plt.contour(xi,yi,zi,15,linewidths=0.5,colors='k')
    #CS = plt.contourf(xi,yi,zi,15,cmap=plt.cm.jet)
    plt.colorbar() # draw colorbar

    # plot data points.
    plt.scatter(x, y, marker='o', c='b', s=5)
    plt.xlim(-3, 3)
    plt.ylim(-3, 3)
    plt.title('griddata test (%d points)' % npts)
    plt.show()