def demo():
    from pylab import hold, linspace, plot, show
    hold(True)
    #y = [9,6,1,3,8,4,2]
    #y = [9,11,13,3,-2,0,2]
    y = [9,11,2,3,8,0,2]
    #y = [9,9,1,3,8,2,2]
    xeq = linspace(0,1,len(y))
    x = xeq+0
    #x[1],x[-2] = x[0],x[-1]
    #x[1],x[-2] = x[2],x[-3]
    #x[1],x[2] = x[2],x[1]
    #x[1],x[-2] = x[2]-0.001,x[-2]+0.001
    #x[1],x[-2] = x[1]-x[1]/2,x[-1]-x[1]/2
    t = linspace(x[0],x[-1],400)
    plot(xeq,y,':oy')
    plot(t,bspline(y,t,clamp=False),'-.y') # bspline
    plot(t,bspline(y,t,clamp=True),'-y') # bspline

    xt,yt = pbs(x,y,t,clamp=False)
    plot(xt,yt,'-.b') # pbs
    xt,yt = pbs(x,y,t,clamp=True)
    plot(xt,yt,'-b') # pbs
    #xt,yt = pbs(x,y,t,clamp=True, parametric=True)
    #plot(xt,yt,'-g') # pbs
    plot(sorted(x),y,':ob')
    show()

    def zeroPaddData(self,desiredLength,paddmode='zero',where='end'):    
        #zero padds the time domain data, it is possible to padd at the beginning,
        #or at the end, and further gaussian or real zero padding is possible        
        #might not work for gaussian mode!

        desiredLength=int(desiredLength)
        #escape the function        
        if desiredLength<0:
            return 0

        #calculate the paddvectors        
        if paddmode=='gaussian':
            paddvec=py.normal(0,py.std(self.getPreceedingNoise())*0.05,desiredLength)
        else:
            paddvec=py.ones((desiredLength,self.tdData.shape[1]-1))
            paddvec*=py.mean(self.tdData[-20:,1:])
            
        timevec=self.getTimes()
        if where=='end':
            #timeaxis:
            newtimes=py.linspace(timevec[-1],timevec[-1]+desiredLength*self.dt,desiredLength)
            paddvec=py.column_stack((newtimes,paddvec))
            longvec=py.row_stack((self.tdData,paddvec))
        else:
            newtimes=py.linspace(timevec[0]-(desiredLength+1)*self.dt,timevec[0],desiredLength)
            paddvec=py.column_stack((newtimes,paddvec))
            longvec=py.row_stack((paddvec,self.tdData))
            
        self.setTDData(longvec)

def contourFromFunction(XYfunction,plotPoints=100,\
			xrange=None,yrange=None,numContours=20,alpha=1.0, contourLines=None):
	"""
	Given a 2D function, plots constant contours over the given
	range.  If the range is not given, the current plotting
	window range is used.
	"""
	
	# set up x and y ranges
	currentAxis = pylab.axis()
	if xrange is not None:
		xvalues = pylab.linspace(xrange[0],xrange[1],plotPoints)
	else:
		xvalues = pylab.linspace(currentAxis[0],currentAxis[1],plotPoints)
	if yrange is not None:
		yvalues = pylab.linspace(yrange[0],yrange[1],plotPoints)
	else:
		yvalues = pylab.linspace(currentAxis[2],currentAxis[3],plotPoints)
	
	#coordArray = _coordinateArray2D(xvalues,yvalues)
	# add extra dimension to this to make iterable?
	# bug here!  need to fix for contour plots
	z = map( lambda y: map(lambda x: XYfunction(x,y), xvalues), yvalues)
	if contourLines:
		pylab.contour(xvalues,yvalues,z,contourLines,alpha=alpha)
	else:
		pylab.contour(xvalues,yvalues,z,numContours,alpha=alpha)

def griddata( X, Y, Z, xl, yl, xr, yr, dx):
    # define grid.
    xi, yi = p.meshgrid( p.linspace(xl,xr, int((xr-xl)/dx)+1), p.linspace(yl,yr, int((yr-yl)/dx)+1))
    # grid the data.
    zi = mgriddata(X,Y,Z,xi,yi)
    New = grid( zi, xl, yl, dx)
    return New

def data2fig(data, X, options, legend_title, xlabel, ylabel=r'Reachability~$\reachability$'):
    if options['grayscale']:
        colors = options['graycm'](pylab.linspace(0, 1, len(data.keys())))
    else:
        colors = options['color'](pylab.linspace(0, 1, len(data.keys())))
    fig = MyFig(options, figsize=(10, 8), xlabel=r'Sources~$\sources$', ylabel=ylabel, grid=False, aspect='auto', legend=True)
    for j, nhdp_ht in enumerate(sorted(data.keys())):
        d = data[nhdp_ht]
        try:
            mean_y = [scipy.mean(d[n]) for n in X]
        except KeyError:
            logging.warning('key \"%s\" not found, continuing...', nhdp_ht)
            continue
        confs_y = [confidence(d[n])[2] for n in X]
        poly = [conf2poly(X, list(numpy.array(mean_y)+numpy.array(confs_y)), list(numpy.array(mean_y)-numpy.array(confs_y)), color=colors[j])]
        patch_collection = PatchCollection(poly, match_original=True)
        patch_collection.set_alpha(0.3)
        patch_collection.set_linestyle('dashed')
        fig.ax.add_collection(patch_collection)
        fig.ax.plot(X, mean_y, label='$%d$' % nhdp_ht, color=colors[j])
    fig.ax.set_xticks(X)
    fig.ax.set_xticklabels(['$%s$' % i for i in X])
    fig.ax.set_ylim(0,1)
    fig.legend_title = legend_title
    return fig

 def plot_interfaces(current_data):
     from pylab import linspace, plot
     xl = linspace(xp1,xp2,100)
     yl = linspace(yp1,yp2,100)
     plot(xl,yl,'g')
     xl = linspace(xlimits[0],xlimits[1],100)
     plot(xl,0.0*xl,'b')

def test_operation_approx():
    def flux_qubit_potential(phi_m, phi_p):
        return 2 + alpha - 2 * pl.cos(phi_p)*pl.cos(phi_m) - alpha * pl.cos(phi_ext - 2*phi_p)
    alpha = 0.7
    phi_ext = 2 * np.pi * 0.5
    phi_m = pl.linspace(0, 2*np.pi, 100)
    phi_p = pl.linspace(0, 2*np.pi, 100)
    X,Y = pl.meshgrid(phi_p, phi_m)
    Z = flux_qubit_potential(X, Y).T

    # the diagram creatinos
    from diagram.operations.computations import multiply
    from diagram.ternary import AEV3DD

    aevdd = AEV3DD()
    diagram3 = aevdd.create(Z, 0, True)
    diagram4 = aevdd.create(Z, 0, True)

    aevdd_mat = multiply(diagram3, diagram4, 9).to_matrix(77, True)
    aevdd_mat_approx = multiply(diagram3, diagram4, 9, approximation_precision=1, in_place='1').to_matrix(27, True)

    pl.plt.figure()
    fig, ax = pl.plt.subplots()
    p = ax.pcolor(X/(2*pl.pi), Y/(2*pl.pi), Z, cmap=pl.cm.RdBu, vmin=abs(Z).min(), vmax=abs(Z).max())
    cb = fig.colorbar(p, ax=ax)
    p = ax.pcolor(X/(2*pl.pi), Y/(2*pl.pi), Z, cmap=pl.cm.RdBu, vmin=abs(Z).min(), vmax=abs(Z).max())
    cb = fig.colorbar(p, ax=ax)
    # cnt = ax.contour(Z, cmap=pl.cm.RdBu, vmin=abs(Z).min(), vmax=abs(Z).max(), extent=[0, 1, 0, 1])
    pl.show()

 def plotDistribution(self):  
     # plot frequency count for the entire vocabulary
     threshold = 1000
     size = len(self.listOfDict[0])
     x = linspace(1, size, size)
     y = sorted(self.listOfDict[0].values(), reverse=True)
     fig = plt.figure()
     axes = fig.add_axes([0.1, 0.1, 0.8, 0.8])
     axes.plot(x, y, 'r')
     axes.set_xlabel('x')
     axes.set_ylabel('y')
     axes.set_title('title');
     plt.show()
     #plot frequency count for all words with count over a given threshold (e.g. number of files read)
     threshold = self.numFilesRead
     size = len(self.listOfDict[0])
     size_relevant = sum(1 for i in self.listOfDict[0].values() if i>threshold)
     x = linspace(1, size_relevant, size_relevant)
     y = sorted([i for i in self.listOfDict[0].values() if i>threshold], reverse=True)
     fig = plt.figure()
     axes = fig.add_axes([0.1, 0.1, 0.8, 0.8])
     axes.plot(x, y, 'r')
     axes.set_xlabel('x')
     axes.set_ylabel('y')
     axes.set_title('title');
     plt.show()

def gfe4():
  x2=plt.linspace(1e-20,.13,90000)
  xmin2=((4*np.pi*(SW.R)**3)/3)*1e-20
  xmax2=((4*np.pi*(SW.R)**3)/3)*.13
  xff2 = plt.linspace(xmin2,xmax2,90000)
  thigh=100
  plt.figure()
  plt.title('Grand free energy per volume vs ff @ T=%0.4f'%Tlist[thigh])
  plt.ylabel('Grand free energy per volume')
  plt.xlabel('filling fraction')  
  plt.plot(xff2,SW.phi(Tlist[thigh],x2,nR[thigh]),color='#f36118',linewidth=3)
  #plt.axvline(nL[thigh])
  #plt.axvline(nR[thigh])
  #plt.axhline(SW.phi(Tlist[thigh],nR[thigh]))
  #plt.plot(x2,x2-x2,'c')
  plt.plot(nL[thigh]*((4*np.pi*(SW.R)**3)/3),SW.phi(Tlist[thigh],nL[thigh],nR[thigh]),'ko')
  plt.plot(nR[thigh]*((4*np.pi*(SW.R)**3)/3),SW.phi(Tlist[thigh],nR[thigh],nR[thigh]),'ko')
  plt.axhline(SW.phi(Tlist[thigh],nR[thigh],nR[thigh]),color='c',linewidth=2)
  print(Tlist[100])
  print(nL[100],nR[100])
  plt.savefig('figs/gfe_cotangent.pdf')

  plt.figure()
  plt.plot(xff2,SW.phi(Tlist[thigh],x2,nR[thigh]),color='#f36118',linewidth=3)
  plt.plot(nL[thigh]*((4*np.pi*(SW.R)**3)/3),SW.phi(Tlist[thigh],nL[thigh],nR[thigh]),'ko')
  plt.plot(nR[thigh]*((4*np.pi*(SW.R)**3)/3),SW.phi(Tlist[thigh],nR[thigh],nR[thigh]),'ko')
  plt.axhline(SW.phi(Tlist[thigh],nR[thigh],nR[thigh]),color='c',linewidth=2)
  plt.xlim(0,0.0003)
  plt.ylim(-.000014,0.000006)
  print(Tlist[100])
  print(nL[100],nR[100])
  plt.savefig('figs/gfe_insert_cotangent.pdf')

def draw_bandstructure(
    jobname, kspace, band, ext=".csv", format="pdf", filled=True, levels=15, lines=False, labeled=False, legend=False
):
    # clf()
    fig = figure(figsize=fig_size)
    ax = fig.add_subplot(111, aspect="equal")
    x, y, z = loadtxt(jobname + ext, delimiter=", ", skiprows=1, usecols=(1, 2, 4 + band), unpack=True)
    if kspace.dimensions == 1:
        pylab.plot(x, y, z)
    elif kspace.dimensions == 2:
        xi = linspace(-0.5, 0.5, kspace.x_res)
        yi = linspace(-0.5, 0.5, kspace.y_res)
        zi = griddata(x, y, z, xi, yi)
        if filled:
            cs = ax.contourf(xi, yi, zi, levels, **contour_filled)
            legend and colorbar(cs, **colorbar_style)
            cs = lines and ax.contour(xi, yi, zi, levels, **contour_lines)
            labeled and lines and clabel(cs, fontsize=8, inline=1)
        else:
            cs = ax.contour(xi, yi, zi, levels, **contour_plain)
            legend and colorbar(cs, **colorbar_style)
            labeled and clabel(cs, fontsize=8, inline=1)
        ax.set_xlim(-0.5, 0.5)
        ax.set_ylim(-0.5, 0.5)
    savefig(jobname + format, format=format, transparent=True)

def main():
    """
    This shows the use of SynChan with Izhikevich neuron. This can be
    used for creating a network of Izhikevich neurons.
    """
    
    simtime = 200.0
    stepsize = 10.0
    model_dict = make_model()
    vm, inject, gk, spike = setup_data_recording(model_dict['neuron'],
                                          model_dict['pulse'],
                                          model_dict['synapse'],
                                          model_dict['spike_in'])
    mutils.setDefaultDt(elecdt=0.01, plotdt2=0.25)
    mutils.assignDefaultTicks(solver='ee')
    moose.reinit()
    mutils.stepRun(simtime, stepsize)
    pylab.subplot(411)
    pylab.plot(pylab.linspace(0, simtime, len(vm.vector)), vm.vector, label='Vm (mV)')
    pylab.legend()
    pylab.subplot(412)
    pylab.plot(pylab.linspace(0, simtime, len(inject.vector)), inject.vector, label='Inject (uA)')
    pylab.legend()
    pylab.subplot(413)
    pylab.plot(spike.vector, pylab.ones(len(spike.vector)), '|', label='input spike times')
    pylab.legend()
    pylab.subplot(414)
    pylab.plot(pylab.linspace(0, simtime, len(gk.vector)), gk.vector, label='Gk (mS)')
    pylab.legend()
    pylab.show()

def Plot_field_gp():
    r = pl.linspace(-1.*mm,1.*mm,50)
    z = pl.linspace(0,g,50)
    
    X, Y = np.meshgrid(z, r)
    for z,r in zip(np.ravel(X), np.ravel(Y)):
        print z*1000,r*1000,SpaceChargeField(r,0,z,0,0,0.5*mm)*0.001

def mk_grid(llx, ulx, nx, lly, uly, ny):
    # Get the Galaxy info
    #galaxies = mk_galaxy_struc()
    galaxies = pickle.load(open('galaxies.pickle','rb'))
    galaxies = filter(lambda galaxy: galaxy.ston_I > 30., galaxies)
    galaxies = pyl.asarray(filter(lambda galaxy: galaxy.ICD_IH < 0.5, galaxies))

    # Make the low mass grid first
    x = [galaxy.Mass for galaxy in galaxies]
    y = [galaxy.ICD_IH *100 for galaxy in galaxies]
    bins_x =pyl.linspace(llx, ulx, nx)
    bins_y = pyl.linspace(uly, lly, ny)

    grid = []

    for i in range(bins_x.size-1):
        xmin = bins_x[i]
        xmax = bins_x[i+1]
        for j in  range(bins_y.size-1):
            ymax = bins_y[j]
            ymin = bins_y[j+1]

            cond=[cond1 and cond2 and cond3 and cond4 for cond1, cond2, cond3,
                cond4 in zip(x>=xmin, x<xmax, y>=ymin, y<ymax)]

            grid.append(galaxies.compress(cond))

    return grid

def demo():
    from pylab import hold, linspace, subplot, plot, legend, show
    hold(True)
    #y = [9,6,1,3,8,4,2]
    #y = [9,11,13,3,-2,0,2]
    y = [9, 11, 2, 3, 8, 0]
    #y = [9,9,1,3,8,2,2]
    x = linspace(0, 1, len(y))
    t = linspace(x[0], x[-1], 400)
    subplot(211)
    plot(t, bspline(y, t, clamp=False), '-.y',
         label="unclamped bspline")  # bspline
    # bspline
    plot(t, bspline(y, t, clamp=True), '-y', label="clamped bspline")
    plot(sorted(x), y, ':oy', label="control points")
    legend()
    #left, right = _derivs(t, bspline(y, t, clamp=False))
    #print(left, (y[1] - y[0]) / (x[1] - x[0]))

    subplot(212)
    xt, yt = pbs(x, y, t, clamp=False)
    plot(xt, yt, '-.b', label="unclamped pbs")  # pbs
    xt, yt = pbs(x, y, t, clamp=True)
    plot(xt, yt, '-b', label="clamped pbs")  # pbs
    #xt,yt = pbs(x,y,t,clamp=True, parametric=True)
    # plot(xt,yt,'-g') # pbs
    plot(sorted(x), y, ':ob', label="control points")
    legend()
    show()

def plot_elecs_and_neurons(neuron_dict, ext_sim_dict, neural_sim_dict):
    pl.close('all')
    fig_all = pl.figure(figsize=[15,15])
    ax_all = fig_all.add_axes([0.1, 0.1, 0.8, 0.8], frameon=False)
    for elec in xrange(len(ext_sim_dict['elec_z'])):
        ax_all.plot(ext_sim_dict['elec_z'][elec], ext_sim_dict['elec_y'][elec], color='b',\
                marker='$E%i$'%elec, markersize=20 )    
    legends = []
    for i, neur in enumerate(neuron_dict):
        folder = os.path.join(neural_sim_dict['output_folder'], neuron_dict[neur]['name'])
        coor = np.load(os.path.join(folder,'coor.npy'))
        x,y,z = coor
        n_compartments = len(x)
        fig = pl.figure(figsize=[10, 10])
        ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], frameon=False)
        # Plot the electrodes
        for elec in xrange(len(ext_sim_dict['elec_z'])):
            ax.plot(ext_sim_dict['elec_z'][elec], ext_sim_dict['elec_y'][elec], color='b',\
                   marker='$%i$'%elec, markersize=20 )
        # Plot the neuron
        xmid, ymid, zmid = np.load(folder + '/coor.npy')
        xstart, ystart,zstart = np.load(folder + '/coor_start.npy')
        xend, yend, zend = np.load(folder + '/coor_end.npy')
        diam = np.load(folder + '/diam.npy')
        length = np.load(folder + '/length.npy')
        n_compartments = len(diam)
        for comp in xrange(n_compartments):
            if comp == 0:
                xcoords = pl.array([xmid[comp]])
                ycoords = pl.array([ymid[comp]])
                zcoords = pl.array([zmid[comp]])
                diams = pl.array([diam[comp]])    
            else:
                if zmid[comp] < 0.400 and zmid[comp] > -.400:  
                    xcoords = pl.r_[xcoords, pl.linspace(xstart[comp],
                                                         xend[comp], length[comp]*3*1000)]   
                    ycoords = pl.r_[ycoords, pl.linspace(ystart[comp],
                                                         yend[comp], length[comp]*3*1000)]   
                    zcoords = pl.r_[zcoords, pl.linspace(zstart[comp],
                                                         zend[comp], length[comp]*3*1000)]   
                    diams = pl.r_[diams, pl.linspace(diam[comp], diam[comp],
                                                length[comp]*3*1000)]
        argsort = pl.argsort(-xcoords)
        ax.scatter(zcoords[argsort], ycoords[argsort], s=20*(diams[argsort]*1000)**2,
                       c=xcoords[argsort], edgecolors='none', cmap='gray')
        ax_all.plot(zmid[0], ymid[0], marker='$%i$'%i, markersize=20, label='%i: %s' %(i, neur))
        #legends.append('%i: %s' %(i, neur))
        ax.axis(ext_sim_dict['plot_range'])
        ax.axis('equal')
        ax.axis(ext_sim_dict['plot_range'])
        ax.set_xlabel('z [mm]')
        ax.set_ylabel('y [mm]')
        fig.savefig(os.path.join(neural_sim_dict['output_folder'],\
                  'neuron_figs', '%s.png' % neur))
    ax_all.axis('equal')
    ax.axis(ext_sim_dict['plot_range'])
    ax_all.set_xlabel('z [mm]')
    ax_all.set_ylabel('y [mm]')
    ax_all.legend()
    fig_all.savefig(os.path.join(neural_sim_dict['output_folder'], 'fig.png'))

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

def bistability_analysis():
    f2_range = linspace(0, 0.4, 41)
    t = linspace(0, 50000, 1000)
    ion()
    ss_aBax_vals_up = []
    ss_aBax_vals_down = []

    for f2 in f2_range:
        model.parameters['Bid_0'].value = f2 * 1e-1
        bax_total = 2e-1

        # Do "up" portion of hysteresis plot
        model.parameters['aBax_0'].value = 0
        model.parameters['cBax_0'].value = bax_total
        x = odesolve(model, t)
        figure('up')
        plot(t, x['aBax_']/bax_total)
        ss_aBax_vals_up.append(x['aBax_'][-1]/bax_total)

        # Do "down" portion of hysteresis plot
        model.parameters['aBax_0'].value = bax_total
        model.parameters['cBax_0'].value = 0
        x = odesolve(model, t)
        figure('down')
        plot(t, x['aBax_']/bax_total)
        ss_aBax_vals_down.append(x['aBax_'][-1]/bax_total)

    figure()
    plot(f2_range, ss_aBax_vals_up, 'r')
    plot(f2_range, ss_aBax_vals_down, 'g')

    def _plot_bar3d(self):
        logging.debug('')
        if self.dimension > 0:
            return
        array = self._data_to_array()

        # width, depth, and height of the bars (array == height values == dz)
        dx = list(numpy.array([1.0/len(array[0]) for x in range(0, array.shape[1])]))
        dy = list(numpy.array([1.0/len(array) for x in range(0, array.shape[0])]))
        dx *= len(array)
        dy *= len(array[0])
        dz = array.flatten()+0.00000001 # dirty hack to cirumvent ValueError

        # x,y,z position of each bar
        x = pylab.linspace(0.0, 1.0, len(array[0]), endpoint=False)
        y = pylab.linspace(0.0, 1.0, len(array), endpoint=False)
        xpos, ypos = pylab.meshgrid(x, y)
        xpos = xpos.flatten()
        ypos = ypos.flatten()
        zpos = numpy.zeros(array.shape).flatten()

        fig = MyFig(self.options, xlabel='Probability p', ylabel='Fraction of Nodes', zlabel='Fraction of Executions', ThreeD=True)
        fig.ax.set_zlim3d(0.0, 1.01)
        fig.ax.set_xlim3d(0.0, 1.01)
        fig.ax.set_ylim3d(0.0, 1.01)
        fig.ax.set_autoscale_on(False)

        assert(len(dx) == len(dy) == len(array.flatten()) == len(xpos) == len(ypos) == len(zpos))
        fig.ax.bar3d(xpos, ypos, zpos, dx, dy, dz, color=['#CCBBDD'])
        try:
            self.figures['bar3d'] = fig.save('bar3d' + str(self.data_filter))
        except ValueError, err:
            logging.warning('%s', err)

    def plot_wire_surface_pcolor(self):
        """
        Plot the fraction of executions as a function of the fraction of nodes for
        each source. Three plots are created: wireframe, surface, and pseudocolor.

        """
        logging.debug('')
        if self.dimension > 0:
            return
        array = self._data_to_array()

        x = pylab.linspace(0.0, 1.0, len(array[0])+1)
        y = pylab.linspace(0.0, 1.0, len(array)+1)
        X, Y = pylab.meshgrid(x, y)

        #fig_wire = MyFig(self.options, xlabel='Probability p', ylabel='Fraction of Nodes', ThreeD=True)
        #fig_surf = MyFig(self.options, xlabel='Probability p', ylabel='Fraction of Nodes', ThreeD=True)
        fig_pcol = MyFig(self.options, xlabel='Probability p', ylabel='Fraction of Nodes')

        #fig_wire.ax.plot_wireframe(X, Y, array)
        #fig_surf.ax.plot_surface(X, Y, array, rstride=1, cstride=1, linewidth=1, antialiased=True)
        pcolor = fig_pcol.ax.pcolor(X, Y, array, cmap=cm.jet, vmin=0.0, vmax=1.0)
        cbar = fig_pcol.fig.colorbar(pcolor, shrink=0.8, aspect=10)
        cbar.ax.set_yticklabels(pylab.linspace(0.0, 1.0, 11), fontsize=0.8*self.options['fontsize'])

        #for ax in [fig_wire.ax, fig_surf.ax]:
            #ax.set_zlim3d(0.0, 1.01)
            #ax.set_xlim3d(0.0, 1.01)
            #ax.set_ylim3d(0.0, 1.01)

        #self.figures['wireframe'] = fig_wire.save('wireframe_' + str(self.data_filter))
        #self.figures['surface'] = fig_surf.save('surface_' + str(self.data_filter))
        self.figures['pcolor'] = fig_pcol.save('pcolor_' + str(self.data_filter))

def showSF( N, label = '' ):
    fig = P.figure()
    ax1 = fig.add_subplot(1, 2, 1, projection='3d')
    ax2 = fig.add_subplot(1, 2, 2)
    
    Nx = 21
    Ny = 21
    nLevels = 12

    tix = P.linspace( 0.0, 1.0, Nx )
    tiy = P.linspace( 0.0, 1.0, Ny )
    (X,Y) = P.meshgrid( tix, tiy )
    z = g.RVector( len( X.flat ) )
    
    for i, x in enumerate( X.flat ):
        p = g.RVector3( X.flat[ i ], Y.flat[ i ] )
        z[ i ] = N(p)
    
    Z = P.ma.masked_where( z == -99., z )
    Z = Z.reshape( Ny, Nx )
    
    ax2.contourf( X, Y, Z )
    ax2.set_aspect( 'equal' )
    surf = ax1.plot_surface( X, Y, Z, rstride = 1, cstride = 1, cmap=P.cm.jet,linewidth=0 )
    
    ax2.set_title( label + N.__str__() )
    fig.colorbar( surf )