def sanity_exampleExoplanetEU2(self):
   """
     Example of ExoplanetEU2
   """
   from PyAstronomy import pyasl
   import matplotlib.pylab as plt
   
   # Instantiate exoplanetEU2 object
   v = pyasl.ExoplanetEU2()
   
   # Show the available data
   v.showAvailableData()
   print()
   
   # Get a list of all available column names
   acs = v.getColnames()
   print("Available column names: " + ", ".join(acs))
   print()
   
   # Select data by planet name (returns a dictionary)
   print(v.selectByPlanetName("CoRoT-2 b"))
   print()
   
   # Get all data as an astropy table
   at = v.getAllDataAPT()
   
   # Export all data as a pandas DataFrame
   pd = v.getAllDataPandas()
   
   # Plot mass vs. SMA
   plt.title("Mass vs. SMA")
   plt.xlabel("[" + v.getUnitOf("mass") + "]")
   plt.ylabel("[" + v.getUnitOf("semi_major_axis") + "]")
   plt.loglog(at["mass"], at["semi_major_axis"], 'b.')

def experiment_plot( ctr, trials, success ):
	"""
	Pass in the ctr, trials and success returned
	by the `experiment` function and plot
	the Cumulative Number of Turns For Each Arm and
	the CTR's Convergence Plot side by side
	"""
	T, K = trials.shape
	n = np.arange(T) + 1
	fig = plt.figure( figsize = ( 14, 7 ) )

	plt.subplot(121)	
	for i in range(K):
		plt.loglog( n, trials[ :, i ], label = "arm {}".format(i + 1) )

	plt.legend( loc = "upper left" )
	plt.xlabel("Number of turns")
	plt.ylabel("Number of turns/arm")
	plt.title("Cumulative Number of Turns For Each Arm")

	plt.subplot(122)
	for i in range(K):
		plt.semilogx( n, np.zeros(T) + ctr[i], label = "arm {}'s CTR".format( i + 1 ) )

	plt.semilogx( n, ( success[ :, 0 ] + success[ :, 1 ] ) / n, label = "CTR at turn t" )

	plt.axis([ 0, T, 0, 1 ] )
	plt.legend( loc = "upper left" )
	plt.xlabel("Number of turns")
	plt.ylabel("CTR")
	plt.title("CTR's Convergence Plot")

	return fig

    def _test_convergence(self, fun, x, *args):

        for i in x:
            y = fun(i, *args)

            plt.figure('res')
            plt.plot(abs(i), y, '.')
            plt.loglog()

def in_degrees_dist_plot(in_degrees_dist, num_nodes):
    import matplotlib.pylab as plt
    x_axis = []
    y_axis = []
    for node, degree in in_degrees_dist.items():
        if node != 0:
            distribution = float(degree) / float(num_nodes)
            x_axis.append(node)
            y_axis.append(distribution)
    plt.loglog(x_axis, y_axis, 'ro')
    plt.xlabel('In-degrees')
    plt.ylabel('Distribution')
    plt.title('In degrees Distribution (log/log Plot)')
    plt.show()

def plotsa():
	"""
	Plots Simulated annealing example
	"""
	fig1 = pl.figure()
	data = np.loadtxt("simulatedannealing1.csv", skiprows=1, delimiter=",")
	fvals =  data[:, 0]
	nevals = range(len(fvals))
	pl.loglog(nevals, fvals, 'b-')
	pl.xlabel("function evaluations", fontsize=20)
	pl.ylabel("cost function value", fontsize=20)
	pl.ylim([50, 5000])
	ax = fig1.gca()
	ax.tick_params(axis='x', labelsize=16)
	ax.tick_params(axis='y', labelsize=16)
	pl.savefig("simulatedannealing1.png", bbox_inches='tight')

def plotga():
	"""
	Plots genetic algorithm
	"""
	fig1 = pl.figure()
	data = np.loadtxt("geneticalgorithm.csv", skiprows=1, delimiter=",")
	fvals =  data[:, 0]
	nevals = range(len(fvals))
	pl.loglog(nevals, fvals, 'b-')
	pl.xlabel("function evaluations", fontsize=20)
	pl.ylabel("cost function value", fontsize=20)
	pl.ylim([50, 3000])
	ax = fig1.gca()
	ax.tick_params(axis='x', labelsize=16)
	ax.tick_params(axis='y', labelsize=16)
	# pl.tight_layout()
	pl.savefig("geneticalgorithm.png", bbox_inches='tight')	

    def plot_convergence(self,x,y,rate=None,log=True,figname='_plot',case='',title='',tolatex=False):
        self.create_dir(self.plotdir)
        log_name  = ''
        log_label = ''
        log_title = ''
        if not tolatex:
            log_title = 'Log Log '
        if log:
            log_name  = 'loglog_'

        plt.figure()
        if rate is not None:
            p = self.p_line_range
            m = rate[0]
            c = rate[1]
            plt.hold(True)
            plt.loglog(x[p[0]:p[1]],10.0**(m*np.log10(x[p[0]:p[1]])+c-2.0),'r')
        if log:
            plt.loglog(x,y,'bo--')
        else:
            plt.plot(x,y,'o:')
        if not tolatex:
            plt.title(title+log_title+case+' convergence')
        
        plt.xlabel('$'+log_label+'mx$')
        plt.ylabel('$'+log_label+'\Delta q$')

        plt.grid(True)
        plt.draw()
        case = case.replace(" ", "_")
        fig_save_name = figname+'_'+log_name+case+'.'+self.plot_format
        figpath = os.path.join(self.plotdir,fig_save_name)
        plt.savefig(figpath,format=self.plot_format,dpi=320,bbox_inches='tight')
        plt.close()

        if tolatex:
            caption = ''
            if log:
                caption = 'Log Log '
            caption = caption + 'plot of '+case.replace("_"," ")+' convergence test'
            caption.capitalize()
            self.gen_latex_fig(figpath,caption=caption)
        return

 def sanity_example(self):
   """
     Exoplanet EU example
   """
   from PyAstronomy import pyasl
   import matplotlib.pylab as plt
   
   eu = pyasl.ExoplanetEU()
   
   # See what information is available
   cols = eu.availableColumns()
   print cols
   
   print
   # Get all data and plot planet Mass vs.
   # semi-major axis in log-log plot
   dat = eu.getAllData()
   plt.xlabel("Planet Mass [RJ]")
   plt.ylabel("Semi-major axis [AU]")
   plt.loglog(dat.plMass, dat.sma, 'b.')

def log_datafit(x, y, deg):
    z = np.polyfit(np.log10(x), np.log10(y), deg)
    p = np.poly1d(z)
    A = np.zeros(np.shape(p)[0])
    for i in range(np.shape(p)[0]):
        A[::-1][i] = p[i]

    yvals = 0.    
    for j in range(np.shape(p)[0]):
        yvals += (((np.log10(x))**j)*A[::-1][j])

    plt.ion()
    plt.loglog(x, y, 'bo', label='Data')
    plt.loglog(x, 10**(yvals), 'g--', lw=2, label='Best Fit')
    plt.legend(loc='best')
    plt.grid(which='minor')
    plt.minorticks_on()

    print "Ax+B"
    print "A = ", A[0]
    print "B =",  A[1]

def calculate_speed_up():
    N = np.array([80,160,320,640,1280])
    h = 2./N
    cpu_time = np.array([8.97, 33.28 , 141.02 , 554.96, 2164.34])
    gpu_time = np.array([1.4877, 3.1019, 10.46, 33.3699, 125.660])
    speedup = np.array([6.25, 11.43, 13.48, 16.63, 17.22])

    plt.figure()
    plt.loglog(N, cpu_time,'ro-', label='Xeon E5-2650 2.60GHz')
    plt.loglog(N, gpu_time,'b+-', label='Tesla K40')
    plt.loglog(N, (cpu_time[0]*1.2)*(N/N[0])**2, 'k-', label='2rd order')
    #legend(handles=[numerical_solution, order5])
    grid(True)
    ylabel('Runtime for 100 time step (second)')
    xlabel('Domain size')
    title("CPU and GPU Performance")
    xlim([80, 2560])
    xticks([20,40,80,160,320,640,1280, 2560], [20,40,80,160,320,640,1280, 2560])
    legend(bbox_to_anchor=(0., 1, 0.5, .0))

    plt.figure()
    plt.plot(N, speedup,'b+-', label='Tesla K40')
    #loglog(N, (speedup[0]*1.2)*(N/N[0])**2, 'k-', label='2rd order')
    #legend(handles=[numerical_solution, order5])
    grid(True)
    ylabel('Speedup ($t_{cpu}/t_{gpu}$)')
    xlabel('domain size')
    title("GPU Speedup")
    xticks([80,160,320,640,1280], [80,160,320,640,1280])
    legend(bbox_to_anchor=(0., 1, 0.3, .0))
    show()

    def plot_regularization_curve(self, name_add=''):
        # Matplotlib is loaded selectively as it is requires
        # libraries that are often not installed on clusters
        import matplotlib.pylab as plt

        an_name = self.settings.get('an_name', 'xx')
        filename = an_name + '_regu_curve' + name_add + '.png'

        sort_order = np.argsort(self.omega2_list)
        omega2_arr = np.array(self.omega2_list).take(sort_order)
        epe_arr = np.array(self.epe_list).take(sort_order)
        err_arr = np.array(self.err_list).take(sort_order)
        ERR_arr = np.array(self.ERR_list).take(sort_order)

        plt.title('omega2: %.2e Neff: %.1f' % (self.opt_omega2, self.n_eff))

        plt.loglog(omega2_arr, epe_arr, '-', label='EPE')
        plt.loglog(omega2_arr, np.sqrt(err_arr), '--', label='err')
        plt.loglog(omega2_arr, np.sqrt(ERR_arr), '-.', label='ERR')

        plt.legend(loc=2)
        plt.xlabel('Omega2')
        plt.ylabel('eV')
        # plt.show()
        plt.savefig(filename)
        plt.clf()

def plot(out_filepath, graph):

    fig = pylab.figure(1)

    pylab.subplot(211)
    networkx.draw_spring(graph, node_size = 8, with_labels = False)

    deg = {}
    for d in [ len(graph.edges(n)) for n in graph.nodes() ]:
        try:
            deg[d] += 1
        except KeyError:
            deg[d] = 1

    print(deg)
    
    pylab.subplot(212)
    plot = deg.items()

    pylab.loglog([ x[0] for x in plot ], [ x[1] for x in plot ], '.')

    pylab.savefig(out_filepath + '.png')

def many_spectra():
    import healpy as hp
    chunks = ['f','g','h']
    thresh = {'f':1.0, 'g':1.0, 'h':1.8}
    fwhm_deg = {'f':7.0, 'g':7.0, 'h':5.0}
    final_fwhm_deg = {'f':7.0, 'g':7.0, 'h':7.0}
    cl={}
    pl.figure(1)
    pl.clf()
    xlim = [10,700]
    ylim = [1e-7, 3e-4]
    for k in chunks:
        print k
        nmap = get_hpix(quick=True,name=k)
        mask = mask_from_map(nmap, fwhm_deg=fwhm_deg[k], final_fwhm_deg=final_fwhm_deg[k], thresh=thresh[k])
#        mean_nmap = np.mean(nmap[np.where(mask!=0)[0]])
        mean_nmap = np.mean(nmap[np.where(mask>0.5)[0]])
        delta = (nmap-mean_nmap)/mean_nmap
        hp.mollview(hp.smoothing(delta*mask,fwhm=1.*np.pi/180.),title=k,min=-0.3,max=0.3)
        continue
        this_cl = hp.anafast(delta*mask)/np.mean(mask**2.)
        l=np.arange(len(this_cl))
        # smooth in l*cl
        lcl = this_cl*l
        sm_lcl = np.zeros_like(lcl)
        rbox = 15
        for i in l:
            imin = np.max([0,i-rbox])
            imax = np.min([np.max(l), i+rbox])
            sm_lcl[i]=np.mean(lcl[imin:imax])
#        cl[k:this_cl]
#        pl.loglog(l,this_cl,linewidth=2)
        pl.loglog(l,sm_lcl,linewidth=2)
        pdb.set_trace()
    pl.xlim(xlim)
#    pl.ylim(ylim)
    pl.legend(chunks)
    pl.xlabel('L')
    pl.ylabel('L*CL')

def expt1():
	"""
	Experiment 1: Chooses the result files and generates figures
	"""
	# filename = sys.argv[1]
	result_file = "./expt1.txt"
	input_threads, input_sizes, throughputs, resp_times \
		= parse_output(result_file) 

	throughputs_MiB = [tp/2**20 for tp in throughputs]

	fig1 = pl.figure()
	fig1.set_tight_layout(True)
	fig1.add_subplot(221)
	
	pl.semilogx(input_sizes, throughputs_MiB, 
				'bo-', ms=MARKER_SIZE, mew=0, mec='b')
	pl.xlabel("fixed file size (Bytes)")
	pl.ylabel("throughput (MiB/sec)")
	pl.text(2E3, 27, "(A)")

	fig1.add_subplot(222)
	pl.loglog(input_sizes, resp_times, 
			  'mo-', ms=MARKER_SIZE, mew=0, mec='m')
	pl.xlabel("fixed file size (Bytes)")
	pl.ylabel("response time (sec)")
	pl.text(2E3, 500, "(B)")

	fig1.add_subplot(223)
	pl.semilogx(resp_times, throughputs_MiB, 
				'go-', ms=MARKER_SIZE, mew=0, mec='g')
	pl.xlabel("response time(sec)")
	pl.ylabel("throughput (MiB/sec)")
	pl.text(0.2, 27, "(C)")

	pl.tight_layout()
	pl.savefig("./figures/%s" % result_file.replace(".txt", ".pdf"))

def plot_reusage(db, keynames, save_path, attr_name = [ 'linkWeightDistr_x', 'linkWeightDistr_y' ]):
	plt.clf()
	plt.figure(figsize = (8, 5))
	plt.loglog(db[keynames['mog']][attr_name[0]], db[keynames['mog']][attr_name[1]], 'b-', lw = 5, label = 'fariyland')
	plt.loglog(db[keynames['mblg']][attr_name[0]], db[keynames['mblg']][attr_name[1]], 'r:', lw = 5, label = 'twitter')
	plt.loglog(db[keynames['im']][attr_name[0]], db[keynames['im']][attr_name[1]], 'k--', lw = 5, label = 'yahoo')
	plt.xlabel('Usage (days)')
	plt.ylabel('CCDF')
	plt.title('Usage of Links')
	plt.grid(True)
	plt.legend(('fairyland', 'twitter', 'yahoo'), loc = 'best')

	plt.savefig(os.path.join(save_dir, save_path))

def plot_ngals_vs_match_radius():
    f = open(datapath+'num_sources_v_match_radius.txt','r')
    f.next()
    rad = []
    n = []
    for line in f:
        tmp = line.split()
        rad.append(float(tmp[0]))
        n.append(int(tmp[1]))
    pl.clf()
    pl.loglog(rad,n,'bo')
    pl.loglog(rad,n,'b')
    pl.loglog(rad, 550.*(np.array(rad)/15.)**2. + 990.,'r')
    pl.ylim(800,np.max(n))
    f.close()

def createAndSaveLogLogPlot(xData, yData, figFileRoot, xLabel="", yLabel="",
                            fileExt='.png', xMin=-1, xMax=-1, yMin=-1,
                            yMax=-1, plotType='bo', axisFontSize=20,
                            tickFontSize=16, svgFlag=0):

    figFileName = figFileRoot + fileExt
    xData = convert_list_to_array(xData)
    yData = convert_list_to_array(yData)
    tempPlot = plt.loglog(xData, yData, plotType, hold="False")
    plt.xlabel(xLabel, fontsize=axisFontSize)
    plt.ylabel(yLabel, fontsize=axisFontSize)
    ax = plt.gca()
    for tick in ax.xaxis.get_major_ticks():
        tick.label1.set_fontsize(tickFontSize)
    for tick in ax.yaxis.get_major_ticks():
        tick.label1.set_fontsize(tickFontSize)
    if xMin == -1:
        xMin = min(xData.tolist())
    if xMax == -1:
        xMax = max(xData.tolist())
    if yMin == -1:
        yMin = min(yData.tolist())
        if isnan(yMin):
            yMin = 0
    if yMax == -1:
        yMax = 0
        yDataList = yData.tolist()
        for item in yDataList:
            if not isnan(item):
                if item > yMax:
                    yMax = item
    plt.xlim(xMin, xMax)
    plt.ylim(yMin, yMax)
    plt.savefig(figFileName, dpi=150)
    if svgFlag == 1:
        figFileName = figFileRoot + '.svg'
        plt.savefig(figFileName, dpi=150)
    plt.clf()

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.mlab as mlab
import os
import matplotlib.pylab as plt
from math import*
os.chdir('/Users/ronaldholt/Desktop/ORNL/SANS_NCBD_Reduced')
for filename in os.listdir("."):
    if filename.endswith(".txt"):
    	x,y=np.loadtxt(filename, skiprows=1, usecols=(0,1), delimiter=",", unpack=True)
    	x=x[27:91]
    	z=y[27:91]
    	plt.loglog(x,z,label=filename)
plt.legend(loc='lower left')
plt.title("I(Q) vs. Q")
plt.xlabel("Q")
plt.xlim([.1,.3])
plt.ylabel("I(Q)")
plt.show()

#.........这里部分代码省略.........
    r = np.power(np.multiply(3./4*N/(np.pi*rho),np.array([random.uniform(0,1) for _ in range(N)])),1./3)   #create a power law distribution of neuron radii
    r.sort()
    if efield:
        rijk = [[random.uniform(0,1)-0.5 for _ in range(N)],[random.uniform(0,1)-0.5 for _ in range(N)],[random.uniform(0,1)-0.5 for _ in range(N)]] #create vector direction of field 
        #if plotAll:
        #    vi = pylab.plot(rijk[0])
        #    vj = pylab.plot(rijk[1])
        #    vk = pylab.plot(rijk[2])
        #    pylab.show()
    R3 = 0.96e3
    C3 = 2.22e-6
    C2 = 9.38e-9
    C3 = 1.56e-6
    C2 = 9.38e-9
    R4 = 100.e6
    R2N = np.multiply(1./(4*np.pi*epsilon),r)
    R1 = 2100.;
    t_impulse = np.array([dt*n for n in range(100)])

    log.info('initialization complete')

    Vt = pylab.zeros(len(times))
    Vi = Vt
    Vj = Vt
    Vk = Vt

    # start simulation
    #-------------------------------------------------------------------------------#
    for neuron in range(N):
        R2 = R2N[neuron]
        ppwave = pylab.zeros(len(times))
        if BGsim:
            absoluteTimes = np.random.exponential(1./(maxrate*STNdata[0]),1)
        else:
            if len(distributionParameter)  ==  1:
                absoluteTimes = np.random.exponential(1./(distributionParameter[0]),1)
            else:
                absoluteTimes = [random.weibullvariate(distributionParameter[0],distributionParameter[1])]
        while absoluteTimes[-1] < times[-1]-currentLength*dt:
            wave_start = int(absoluteTimes[-1]/dt)
            wave_end = wave_start+currentLength
            if wave_end > len(times):
                break
            ppwave[wave_start:wave_end] = np.add(ppwave[wave_start:wave_end],It)
            if BGsim:
                isi = np.random.exponential(1./(maxrate*STNdata[int(absoluteTimes[-1]/BGdt)]),1)
            else:
                if len(distributionParameter)  ==  1:
                    isi = np.random.exponential(1./(distributionParameter[0]),1)
                else:
                    isi = random.weibullvariate(distributionParameter[0],distributionParameter[1])
            absoluteTimes = np.append(absoluteTimes,[absoluteTimes[-1]+isi])
        # calculate neuron contribution
        #------------------------------------------------------------------------------#
        extracellular_impulse_response = np.multiply(np.multiply(np.exp(np.multiply(t_impulse,-20*17*((C2*R1*R2 + C2*R1*R3 + C2*R1*R4 - C3*R1*R3 + C3*R2*R3 + C3*R3*R4))/(2*C2*C3*R1*R3*(R2 + R4)))),(np.add(np.cosh(np.multiply(t_impulse,(C2**2*R1**2*R2**2 + 2*C2**2*R1**2*R2*R3 + 2*C2**2*R1**2*R2*R4 + C2**2*R1**2*R3**2 + 2*C2**2*R1**2*R3*R4 + C2**2*R1**2*R4**2 + 2*C2*C3*R1**2*R2*R3 - 2*C2*C3*R1**2*R3**2 + 2*C2*C3*R1**2*R3*R4 - 2*C2*C3*R1*R2**2*R3 - 2*C2*C3*R1*R2*R3**2 - 4*C2*C3*R1*R2*R3*R4 - 2*C2*C3*R1*R3**2*R4 - 2*C2*C3*R1*R3*R4**2 + C3**2*R1**2*R3**2 - 2*C3**2*R1*R2*R3**2 - 2*C3**2*R1*R3**2*R4 + C3**2*R2**2*R3**2 + 2*C3**2*R2*R3**2*R4 + C3**2*R3**2*R4**2)**(1/2)/(2*C2*C3*R1*R3*(R2 + R4)))),np.divide(np.sinh(np.multiply(t_impulse,(C2**2*R1**2*R2**2 + 2*C2**2*R1**2*R2*R3 + 2*C2**2*R1**2*R2*R4 + C2**2*R1**2*R3**2 + 2*C2**2*R1**2*R3*R4 + C2**2*R1**2*R4**2 + 2*C2*C3*R1**2*R2*R3 - 2*C2*C3*R1**2*R3**2 + 2*C2*C3*R1**2*R3*R4 - 2*C2*C3*R1*R2**2*R3 - 2*C2*C3*R1*R2*R3**2 - 4*C2*C3*R1*R2*R3*R4 - 2*C2*C3*R1*R3**2*R4 - 2*C2*C3*R1*R3*R4**2 + C3**2*R1**2*R3**2 - 2*C3**2*R1*R2*R3**2 - 2*C3**2*R1*R3**2*R4 + C3**2*R2**2*R3**2 + 2*C3**2*R2*R3**2*R4 + C3**2*R3**2*R4**2)**(1/2)/(2*C2*C3*R1*R3*(R2 + R4))))*(C2*R1*R2 - C2*R1*R3 + C2*R1*R4 + C3*R1*R3 - C3*R2*R3 - C3*R3*R4),(C2**2*R1**2*R2**2 + 2*C2**2*R1**2*R2*R3 + 2*C2**2*R1**2*R2*R4 + C2**2*R1**2*R3**2 + 2*C2**2*R1**2*R3*R4 + C2**2*R1**2*R4**2 + 2*C2*C3*R1**2*R2*R3 - 2*C2*C3*R1**2*R3**2 + 2*C2*C3*R1**2*R3*R4 - 2*C2*C3*R1*R2**2*R3 - 2*C2*C3*R1*R2*R3**2 - 4*C2*C3*R1*R2*R3*R4 - 2*C2*C3*R1*R3**2*R4 - 2*C2*C3*R1*R3*R4**2 + C3**2*R1**2*R3**2 - 2*C3**2*R1*R2*R3**2 - 2*C3**2*R1*R3**2*R4 + C3**2*R2**2*R3**2 + 2*C3**2*R2*R3**2*R4 + C3**2*R3**2*R4**2)**(1/2))))),-R4/(C2*(R2 + R4)));
        electrode_ppwave = np.convolve(ppwave,extracellular_impulse_response,'same');
        if efield:  #add fields
            amp = 1/np.sqrt((np.square(rijk[0][neuron])+np.square(rijk[1][neuron])+np.square(rijk[2][neuron])))
            rijk[0][neuron] = rijk[0][neuron]*amp
            rijk[1][neuron] = rijk[1][neuron]*amp
            rijk[2][neuron] = rijk[2][neuron]*amp
            Vi = np.add(Vi,np.multiply(electrode_ppwave,rijk[0][neuron]))
            Vj = np.add(Vj,np.multiply(electrode_ppwave,rijk[1][neuron]))
            Vk = np.add(Vk,np.multiply(electrode_ppwave,rijk[2][neuron]))
        else:       #add scalar
            Vt = np.add(Vt,electrode_ppwave)
        if np.mod(neuron,1000) == 999:
            log.info(str(neuron+1)+" neurons calculated")
    #------------------------------------------------------------------------------#        
    # end simulation
    
    log.info('neuron contribution to MER complete')
    
    #remove bias
    if efield:
        Vt = np.sqrt(np.add(np.square(Vi),np.square(Vj),np.square(Vk)))   
    Vt = np.subtract(Vt,np.mean(Vt))

    #apply hardware filters
    flow = 5500*2.
    fhigh = 500.
    b,a = signal.butter(18,flow*dt,'low')
    Vt = signal.lfilter(b, a, Vt)
    b,a = signal.butter(1,fhigh*dt,'high')
    Vt = signal.lfilter(b, a, Vt)

    #produce plots
    if plotAll:
        volts = pylab.plot(times,Vt)
        if BGsim:
            stnrate = pylab.plot(Ratetime,np.multiply(STNdata,200))
        pylab.show()
        nfft=2**int(math.log(len(Vt),2))+1
        sr = 1/dt
        Pxi,freqs=pylab.psd(x=Vt,Fs=sr,NFFT=nfft/10,window=pylab.window_none, noverlap=100)
        pylab.show()
        return freqs, Pxi
        psd = pylab.loglog(freqs, Pxi)
        pylab.show()
    return Vt, times

                                           dx=dx, detrend=True, tap=tap)

        try:

          SS_obs = SS_obs + numpy.interp(f0, ffo, PSD_obs)
          SS_model = SS_model + numpy.interp(f0, ffm, PSD_model)
        except:
          SS_obs = numpy.interp(f0, ffo, PSD_obs)
          SS_model = numpy.interp(f0,ffm, PSD_model)
        nr += 1

SS_obs /= nr
SS_model /= nr

ff = f0
dff = ff[1]-ff[0]

plt.close()
plt.figure()
plt.loglog(ff, SS_model, color='red',lw=2, label='ssh_model')
plt.loglog(ff, SS_obs, color='k', label='ssh_obs')
plt.grid()
#plt.axis([5e-3,0.25,1e-4,1e3])
plt.xlabel(u'cy/km')
plt.ylabel(u'm²/(cy/km)')
plt.legend()
plt.title('SSH spectra for SWOT-like data and model data interpolated on the swath')
plt.savefig('{}_ssh_spectra.png'.format(p.config))
plt.show()