def function(self,E) :
        """ Calculates the number of counts in barns"""
        
        if self.delta.value != self._previous_delta:
            self._previous_delta = copy.copy(self.delta.value)
            self.integrategos(self.delta.value)
            self.calculate_knots()

        if self._previous_effective_angle != self.effective_angle.value:
            self.integrategos()
            
        factor = 4.0 * np.pi * a0 ** 2.0 * R**2 / E / self.T #to convert to m**2/bin
        Emax = self.energyaxis[-1] + self.edgeenergy + \
        self.delta.value #maximum tabulated energy
        cts = np.zeros((len(E)))
        
        if self.fs_state is True:
            if self.__knots[-1] > Emax : Emax = self.__knots[-1]
            fine_structure_indices=np.logical_and(np.greater_equal(E, 
            self.edgeenergy+self.delta.value), 
            np.less(E, self.edgeenergy + self.delta.value + self.fs_emax))
            tabulated_indices = np.logical_and(np.greater_equal(E, 
            self.edgeenergy + self.delta.value + self.fs_emax), 
            np.less(E, Emax))
            if self.fs_mode == "new_spline" :
                cts = np.where(fine_structure_indices, 
                1E-25*splev(E,(self.__knots,self.fslist.value,3),0), cts)
            elif self.fs_mode == "spline" :
                cts = np.where(fine_structure_indices, 
                cspline1d_eval(self.fslist.value, 
                E, 
                dx = self.energy_scale / self.knots_factor, 
                x0 = self.edgeenergy+self.delta.value), 
                cts)
            elif self.fs_mode == "spline_times_edge" :
                cts = np.where(fine_structure_indices, 
                factor*splev((E-self.edgeenergy-self.delta.value), 
                self.__goscoeff)*cspline1d_eval(self.fslist.value, 
                E,dx = self.energy_scale / self.knots_factor, 
                x0 = self.edgeenergy+self.delta.value), 
                cts )
        else:
            tabulated_indices = np.logical_and(np.greater_equal(E, 
            self.edgeenergy + self.delta.value), np.less(E, Emax))            
        powerlaw_indices = np.greater_equal(E,Emax)  
        cts = np.where(tabulated_indices, 
        factor * splev((E-self.edgeenergy-self.delta.value), 
        self.__goscoeff),
         cts)
        
        # Convert to barns/dispersion.
        #Note: The R factor is introduced in order to give the same value
        # as DM, although it is not in the equations.
        cts = np.where(powerlaw_indices, self.A * E**-self.r, cts) 
        return (self.__subshell_factor * self.intensity.value * self.energy_scale 
        * 1.0e28 / R) * cts       

def TransformSpectrum(args):
    flux, invvar, c0, c1, newwave = args
    smoother = 3.0

    fitc = cspline1d(flux, lamb=smoother)
    fitc_iv = cspline1d(invvar, lamb=smoother)
    newf = cspline1d_eval(fitc, newwave, dx=c1, x0=c0)
    newiv   = cspline1d_eval(fitc_iv, newwave, dx=c1, x0=c0)

    return (newf, newiv)

def splineIntLinExt(y,x,xNew,splinecoeff = 0.):
    """
    Use the scipy spline interpolation, but linearly extrapolate at the edges,
    since scipy.signal.cspline1d assumes periodic boundary conditions
    """    
    if len(x) < 4:
        return interpolateLin(y,x,xNew)

    if isinstance(xNew,float):
        wasfloat = 1
        xNew = M.array([xNew])
    else:
        wasfloat = 0

    whereSpline = N.where((xNew > x[0]) * (xNew < x[-1]))[0]
    whereLin = N.where((xNew <= x[0]) + (xNew >= x[-1]))[0]
    
    ans = xNew * 0.
    if len(whereSpline) > 0:
        if isinstance(splinecoeff,float): # not pre-calculated.
            splinecoeff = SS.cspline1d(y)
        ans[whereSpline] = SS.cspline1d_eval(splinecoeff, xNew[whereSpline], dx=x[1]-x[0], x0 = x[0])

    if len(whereLin) > 0:
        ans[whereLin] = interpolateLin(y,x,xNew[whereLin])

    if wasfloat:
        return ans[0]
    else:
        return ans

def splineIntBig0LittleLog(y,x,xNew,splinecoeff = 0.):
    """
    Use the scipy spline interpolation, but linearly extrapolate at the edges,
    since scipy.signal.cspline1d assumes periodic boundary conditions
    """
    if len(x) < 4:
        return interpolateLin(y,x,xNew)

    if isinstance(xNew,float):
        wasfloat = 1
        xNew = N.array([xNew])
    else:
        wasfloat = 0

    whereSpline = N.where((xNew >= x[0]) * (xNew <= x[-1]))[0]
    whereLittle = N.where(xNew < x[0])[0]
    whereBig = N.where(xNew >= x[-1])[0]
    
    ans = xNew * 0.
    if len(whereSpline) > 0:
        if isinstance(splinecoeff,float): # not pre-calculated.
            splinecoeff = SS.cspline1d(y)
        ans[whereSpline] = SS.cspline1d_eval(splinecoeff, xNew[whereSpline], dx=x[1]-x[0], x0 = x[0])

    if len(whereLittle) > 0:
        xw = xNew[whereLittle]
        logx,logy = N.log(x[:2]),N.log(y[:2])
        ans[whereLittle] = N.exp(logy[0] + (N.log(xw)-logx[0])/(logx[1]-logx[0])*(logy[1]-logy[0]))


    if wasfloat:
        return ans[0]
    else:
        return ans

 def f( t, *args ):
     for i,arg in enumerate(args): params[ free_params[i] ] = arg
     tshift = params[-1]
     ideal = fmodel( t, *args )
     irf = cspline1d_eval( self.irf_generator, t-tshift, dx=self.irf_dt, x0=self.irf_t0 )
     convoluted = pylab.real(pylab.ifft( pylab.fft(ideal)*pylab.fft(irf) )) # very small imaginary anyway
     return convoluted

    def set_irf( self, irf=None, wraptime=None, dispersion=None ):
        """
        The detector response isn't a delta-function, meaning that what
        you measure isn't the true time-dependence of the system you are
        measuring. It's the time dependence of the system convolved with
        the response of the detector. This method sets the measured
        trace of the detector response that will be used to convolve
        the mult-exponential model before fitting (thereby taking this
        convolution into account without doing nearly-impossible
        numerical deconvolution).
        """
        if isinstance( irf, Trace ):
            self.irf = irf
        elif type(irf) == str:
            self.irf = Trace( irf )
            if wraptime is not None:
                self.irf.wrapcurves( wraptime )
            elif self.wraptime is not None:
                self.irf.wrapcurves( self.wraptime )
        
        self.irf_dispersion = dispersion
        if dispersion is not None:
            # this is meant to address chromatic dispersion within the setup
            # (e.g. optical fiber)
            # don't bother with normalization b/c it gets normalized to unit area below anyway.
            #original = self.irf.curves[0].copy()
            #self.irf.curves[0][dispersion:] += original[:-dispersion]
            #self.irf.curves[0][:-dispersion] += original[dispersion:]
            len1 = len(self.irf.curves[0])
            chain_of_three = pylab.zeros( 3*len1 ) # stack three curves end-to-end so cspline1d_eval doesn't have to extrapolate beyond data
            chain_of_three[:len1] = self.irf.curves[0][:]
            chain_of_three[len1:2*len1] = self.irf.curves[0][:]
            chain_of_three[-len1:] = self.irf.curves[0][:]
            g = cspline1d(chain_of_three)
            smoothed = pylab.zeros( len1 )
            std_dev = dispersion/1000.0
            for t0 in pylab.linspace(-2*std_dev, 2*std_dev, 50):
                weight = pylab.exp( -t0**2/2.0/std_dev**2 )
                smoothed += weight * cspline1d_eval( g, self.irf.t[0]-t0, dx=self.irf.t[0][1], x0=-self.irf.t[0][-1] )
            self.irf.curves[0] = smoothed
            
        normalized = self.irf.curves[0].astype(numpy.float)/float(sum(self.irf.curves[0])) # normalize integral to 1, just like delta-function!!!
        self.irf.curves[0] = normalized.copy()

        self.irf_generator = cspline1d(self.irf.curves[0])
        self.irf_dt = self.irf.t[0][1]-self.irf.t[0][0]
        self.irf_t0 = self.irf.t[0][0]
        
        if False:
            """not sure this matters if we do interpolation
            """
            # difference in degree of binning (e.g. 8ps vs. 4ps is bin difference of 2)
            bin_difference = pylab.np.int( self.resolution / self.irf.resolution )
            if bin_difference != 1:
                raise ValueError("Have not yet tested deconvolution with different resolution than detector trace!!!")
                d = self.irf.curves[0]
                detector_binned = pylab.zeros( len(d)/bin_difference )
                for i in range( len(detector_binned ) ):
                    detector_binned[i] = sum( d[i*bin_difference : i*bin_difference+bin_difference] )

def interp_T(T,Ts,PPs):
    xvals=zeros(1)
    xvals[0]=T
    Ps=list()
    for i in PPs:
        cj=cspline1d(i)
        Ps.append(cspline1d_eval(cj,xvals,dx=500.0,x0=Ts[0]))
    return Ps

def wiggle(x, origin=0, posFill='black', negFill=None, lineColor='black', 
        resampleRatio=10, rescale=False, ymin=0, ymax=None, ax=None):
    """Plots a "wiggle" trace
    Input:
        x: input data (1D numpy array)
        origin: (default, 0) value to fill above or below (float)
        posFill: (default, black) color to fill positive wiggles with (string 
            or None)
        negFill: (default, None) color to fill negative wiggles with (string 
            or None)
        lineColor: (default, black) color of wiggle trace (string or None)
        resampleRatio: (default, 10) factor to resample traces by before 
            plotting (1 = raw data) (float)
        rescale: (default, False) If True, rescale "x" to be between -1 and 1
        ymin: (default, 0) The minimum y to use for plotting
        ymax: (default, len(x)) The maximum y to use for plotting
        ax: (default, current axis) The matplotlib axis to plot onto
    Output:
        a matplotlib plot on the current axes
    """
    from matplotlib import pyplot as plt
    from scipy.signal import cspline1d, cspline1d_eval

    if ymax is None:
        ymax = x.size

    # Rescale so that x ranges from -1 to 1
    if rescale:
        x = x.astype(np.float)
        x -= x.min()
        x /= x.ptp()
        x *= 2
        x -= 1

    # Interpolate at resampleRatio x the previous density
    y = np.linspace(0, x.size, x.size)
    interp_y = np.linspace(0, x.size, x.size * resampleRatio)
    cj = cspline1d(x)
    interpX = cspline1d_eval(cj,interp_y) #,dx=1,x0=0
    newy = np.linspace(ymax, ymin, interp_y.size)
    if origin == None: 
        origin = interpX.mean()

    # Plot
    if ax is None:
        ax = plt.gca()
        plt.hold(True)
    if posFill is not None: 
        ax.fill_betweenx(newy, interpX, origin,
                where=interpX > origin,
                facecolor=posFill)
    if negFill is not None:
        ax.fill_betweenx(newy, interpX, origin,
                where=interpX < origin,
                facecolor=negFill)
    if lineColor is not None:
        ax.plot(interpX, newy, color=lineColor)

    def test_basic(self):
        y=array([1,2,3,4,3,2,1,2,3.0])
        x=arange(len(y))
        dx=x[1]-x[0]
        cj = signal.cspline1d(y)

        x2=arange(len(y)*10.0)/10.0
        y2=signal.cspline1d_eval(cj, x2, dx=dx,x0=x[0])

        # make sure interpolated values are on knot points
        assert_array_almost_equal(y2[::10], y, decimal=5)

def doresample(orig_x, orig_y, new_x, method='cubic', padlen=0, antialias=False):
    """
    Resample data from one spacing to another.  By default, does not apply any antialiasing filter.

    Parameters
    ----------
    orig_x
    orig_y
    new_x
    method
    padlen

    Returns
    -------

    """
    pad_y = tide_filt.padvec(orig_y, padlen=padlen)
    tstep = orig_x[1] - orig_x[0]
    if padlen > 0:
        pad_x = np.concatenate((np.arange(orig_x[0] - padlen * tstep, orig_x[0], tstep),
                                orig_x,
                                np.arange(orig_x[-1] + tstep, orig_x[-1] + tstep * (padlen + 1), tstep)))
    else:
        pad_x = orig_x
    if padlen > 0:
        print('padlen=', padlen)
        print('tstep=', tstep)
        print(pad_x)

    # antialias and ringstop filter
    init_freq = len(pad_x) / (pad_x[-1] - pad_x[0])
    final_freq = len(new_x) / (new_x[-1] - new_x[0])
    if antialias and (init_freq > final_freq):
        aafilterfreq = final_freq / 2.0
        aafilter = tide_filt.noncausalfilter(filtertype='arb', usebutterworth=False)
        aafilter.setarb(0.0, 0.0, 0.95 * aafilterfreq, aafilterfreq)
        pad_y = aafilter.apply(init_freq, pad_y)

    if method == 'cubic':
        cj = signal.cspline1d(pad_y)
        return tide_filt.unpadvec(
            np.float64(signal.cspline1d_eval(cj, new_x, dx=(orig_x[1] - orig_x[0]), x0=orig_x[0])), padlen=padlen)
    elif method == 'quadratic':
        qj = signal.qspline1d(pad_y)
        return tide_filt.unpadvec(
            np.float64(signal.qspline1d_eval(qj, new_x, dx=(orig_x[1] - orig_x[0]), x0=orig_x[0])), padlen=padlen)
    elif method == 'univariate':
        interpolator = sp.interpolate.UnivariateSpline(pad_x, pad_y, k=3, s=0)  # s=0 interpolates
        return tide_filt.unpadvec(np.float64(interpolator(new_x)), padlen=padlen)
    else:
        print('invalid interpolation method')
        return None

 def run(self):
     _, fileExtension = os.path.splitext(self.fileName)
     if fileExtension == '.gmv':
         print('Geomagnetic variation')
         with open(self.fileName, 'rt') as csvdata:
             date = []
             value = []
             for row in csv.reader(csvdata):
                 if ('#' in row[0]):
                     self.header.append(row)
                 else:
                     date.append(row[0])
                     value.append(row[1])
         self.notifyProgress.emit(20)
     elif fileExtension == '.ske':
         print('Kp estimation')
         with open(self.fileName, 'rt') as csvdata:
             date = []
             value = []
             for row in csv.reader(csvdata, delimiter=' '):
                 if ('#' in row[0]):
                     self.header.append(row)
                 else:
                     print(row)
                     if int(row[7]) < 2:
                         date.append(
                             dt.datetime.strptime(
                                 ''.join((row[0], row[1], row[2],
                                         row[4])),
                                 '%Y%m%d%H%M')),
                         value.append(float(row[-1])-float(row[-14]))  # 4h
                         # value.append(float(row[-1])-float(row[19]))  # 1h
         self.notifyProgress.emit(20)
     signal_src = np.array((date, value), dtype=np.dtype('a25'))
     signal = signal_src[:, np.logical_not(
         np.isnan(signal_src[1, :].astype(np.float)))]
     self.notifyProgress.emit(60)
     if self.interpolate:
         self.time = signal_src[0, :].astype(np.datetime64).astype(
             dt.datetime)
         dx = dates.date2num(self.time[1])-dates.date2num(self.time[0])
         cj = cspline1d(signal[1, :].astype(float))
         self.value = cspline1d_eval(cj, dates.date2num(self.time),
                                     dx=dx,
                                     x0=dates.date2num(self.time[0]))
     else:
         self.time = dates.signal[0, :].astype(np.datetime64).astype(
             dt.datetime)
         self.value = signal[1, :].astype(np.float)
     self.notifyProgress.emit(80)
     self.loaded.emit()

def interp_around(X_sc,s_fracpeak,s_before,s_after,kind='cubic'):
    n_c = X_sc.shape[1]
    n_s = s_before+s_after
    Out_sc = np.empty((n_s,n_c),dtype=np.float32)
    for i_c in xrange(n_c):
        if kind == 'cubic':
            coeffs = cspline1d(X_sc[:,i_c])
            Out_sc[:,i_c] = cspline1d_eval(coeffs,
                                       newx=np.arange(s_fracpeak - s_before,s_fracpeak+s_after,dtype=np.float32))
        elif kind == "linear":
            Out_sc[:,i_c] = interp1d(np.arange(X_sc.shape[0]),X_sc[:,i_c],
                                     bounds_error=True,kind=kind)(np.arange(s_fracpeak - s_before,s_fracpeak+s_after,dtype=np.float32))
        else: raise Exception("kind must be 'linear' or 'cubic'")
    return Out_sc

def interpol(x1,y1,x_out,plot=False):
	from scipy.signal import cspline1d, cspline1d_eval
	#assumes that data points are evenly spaced

	dx=x1[1]-x1[0]
	cj=cspline1d(y1)
	y_out=cspline1d_eval(cj,x_out,dx=dx,x0=x1[0])
	if plot:
		from pylab import plot, show,legend
		plot(x_out,y_out,'ob',x1,y1,'xg',ms=2)
		legend(['interpolated','original'])
		show()

	return y_out

def SplineResample(y, n_x, x_range = None, smoother = 3):
    ''' Resample a signal y using spline.

    y: the signal must be sampled on uniform x
    n_x: number of points for the new signal
    x_range: tuple (start x, end x) of the signal
    smoother: spline smoothing strength
    '''

    if x_range is None: x_range = (0., 1.)
    spcoeffs = cspline1d(y, lamb = smoother)
    return cspline1d_eval(
        spcoeffs, linspace(x_range[0], x_range[1], n_x), 
        dx = (x_range[1]-x_range[0])/(len(y)-1.0), x0 = x_range[0])

def evalBeams(arrBeams,guessShearLoad,spc,debugFlag = False):
    
    lastBeam = arrBeams[-1]
    arrStrainEnergy = np.zeros(np.size(arrBeams))
    arrSurfEnergy = np.zeros(np.size(arrBeams))
    lastIndex = np.size(arrBeams)
    maxSurfEnergy = -gamma * lastBeam.w*(lastBeam.Lt - lastBeam.L)
    arrResults = np.array([])
    
    for j,beam in enumerate(arrBeams):
        if debugFlag:
            print('Beam Length: %f of %f'%(beam.L*scale,beam.Lt*scale))
        
        dispSearch(beam=beam,initLoad = guessShearLoad,goal=spc/2/scale,tol=1e-7,right=0,left=0)
        guessShearLoad = beam.shearLoad
        
        if debugFlag:
            print('Solved Beam -- TipDisp: %s Goal: %s Force: %s' % (beam.yTipDisplacement()*scale,spc/2,beam.shearLoad))
        
        arrStrainEnergy[j] = beam.calculateStrainEnergy()
        arrSurfEnergy[j] = -gamma * beam.w *(beam.Lt - beam.L)
        if arrStrainEnergy[j] >= np.abs(maxSurfEnergy): 
            # since there is more bending energy than surface energy stop computing 
            print('Super stiff beam')
            lastIndex = j
            break
    
    if lastIndex > 0:   # This ensures that we have more than one data point before trying to interpolate
        interpLens = np.linspace(arrBeamLens[0],arrBeamLens[lastIndex-1],num=100,endpoint=True) # Generate x values for which to interpolate
        csFit = cspline1d((arrStrainEnergy[0:lastIndex]+arrSurfEnergy[0:lastIndex]))    # Generate cubic spline fit to the sub dataset
        interpTotalEnergy = cspline1d_eval(csFit,interpLens,dx=(arrBeamLens[1]-arrBeamLens[0]), x0 = arrBeamLens[0])    # Generate the interpolated values from the fit and x points
        finalLen = interpLens[interpTotalEnergy.argmin()]   # find the minimum of the energy balance and grab index to choose the appropriate length
        
        if debugFlag:
            print('beamLens shape: %s arrStrain: %s'%(arrBeamLens[0:lastIndex].shape,arrStrainEnergy[0:lastIndex].shape))
            mpl.figure()
            mpl.hold(True)
            mpl.plot(arrBeamLens[0:lastIndex]*scale,arrStrainEnergy[0:lastIndex]*scale)
            mpl.plot(arrBeamLens[0:lastIndex]*scale,arrSurfEnergy[0:lastIndex]*scale)
            mpl.plot(interpLens*scale,interpTotalEnergy*scale,arrBeamLens[0:lastIndex]*scale,(arrStrainEnergy+arrSurfEnergy)[0:lastIndex]*scale,'o')
        arrResults = np.array([arrBeamLens[0:lastIndex],arrStrainEnergy[0:lastIndex]])
    else:   # since there is only one datapoint then use that as the value
        finalLen = arrBeamLens[lastIndex]
        arrResults = np.array([arrBeamLens[lastIndex],arrStrainEnergy[lastIndex]])
    
    
    
    return (finalLen,arrResults)

def sr_interpol2(x,y,ytarget,doplot=0,factor=10):
    dx = x[1]-x[0]
    newx = linspace(min(x),max(x),factor*len(x))

    cj = cspline1d(y)
    newy = cspline1d_eval(cj, newx, dx=dx,x0=x[0])

    ysq = (ytarget-newy)**2
    index = where(ysq == min(ysq))

    if doplot:
        clf()
        plot(x,y,'o')
        plot(newx,newy)
        plot(newx[index],newy[index],'o')
        show()

    return newx[index[0][0]]

def peakdetect_spline(y_axis, x_axis, pad_len=20):
    """
    Performs a b-spline interpolation on the data to increase resolution and
    send the data to the 'peakdetect_zero_crossing' function for peak 
    detection.
    
    Omitting the x_axis is forbidden as it would make the resulting x_axis
    value silly if it was returned as the index 50.234 or similar.
    
    will find the same amount of peaks as the 'peakdetect_zero_crossing'
    function, but might result in a more precise value of the peak.
    
    keyword arguments:
    y_axis -- A list containing the signal over which to find peaks
    
    x_axis -- A x-axis whose values correspond to the y_axis list and is used
        in the return to specify the position of the peaks. 
        x-axis must be equally spaced.
    
    pad_len -- By how many times the time resolution should be increased by,
        e.g. 1 doubles the resolution.
        (default: 20)
    
    
    return: two lists [max_peaks, min_peaks] containing the positive and
        negative peaks respectively. Each cell of the lists contains a tuple
        of: (position, peak_value) 
        to get the average peak value do: np.mean(max_peaks, 0)[1] on the
        results to unpack one of the lists into x, y coordinates do: 
        x, y = zip(*max_peaks)
    """
    # check input data
    x_axis, y_axis = _datacheck_peakdetect(x_axis, y_axis)
    # could perform a check if x_axis is equally spaced
    #if np.std(np.diff(x_axis)) > 1e-15: raise ValueError
    # perform spline interpolations
    dx = x_axis[1] - x_axis[0]
    x_interpolated = np.linspace(x_axis.min(), x_axis.max(), len(x_axis) * (pad_len + 1))
    cj = cspline1d(y_axis)
    y_interpolated = cspline1d_eval(cj, x_interpolated, dx=dx,x0=x_axis[0])
    # get peaks
    max_peaks, min_peaks = peakdetect_zero_crossing(y_interpolated, x_interpolated)
    
    return [max_peaks, min_peaks]

    def fitSpline(self,degree=2):
        """
        **SUMMARY**

        A function to generate a spline curve fitting over the points in LineScan with
        order of precision given by the parameter degree

        **PARAMETERS**

        * *degree* - the precision of the generated spline 

        **RETURNS**

        The spline as a LineScan fitting over the initial values of LineScan

        **EXAMPLE**

        >>> import matplotlib.pyplot as plt
        >>> img = Image("lenna")
        >>> ls = img.getLineScan(pt1=(10,10)),pt2=(20,20)).normalize()
        >>> spline = ls.fitSpline()
        >>> plt.plot(ls)
        >>> plt.show()
        >>> plt.plot(spline)
        >>> plt.show()
        
        **NOTES**

        Implementation taken from http://www.scipy.org/Cookbook/Interpolation  

        """
        if degree > 4:
            degree = 4  # No significant improvement with respect to time usage
        if degree < 1:
            warnings.warn('LineScan.fitSpline - degree needs to be >= 1')
            return None
        retVal = None
        y = np.array(self)
        x = np.arange(0,len(y),1)
        dx = 1
        newx = np.arange(0,len(y)-1,pow(0.1,degree))
        cj = sps.cspline1d(y)
        retVal = sps.cspline1d_eval(cj,newx,dx=dx,x0=x[0])
        return retVal

def _spline_interpolate(oldx, oldy, newx, smoothing=0.001,fast=True, **kw):
    """
    cubic splines for axis alignment using
    scipy.signal and/or scipy.interpolate

    newy = _spline_interpolate(oldx, oldy, newx, fast=True)
    if fast = True
       1-dimensional cubic spline for cases where
       oldx and newx are on a uniform grid.
    else
       handles multi-dimensional data, non-uniform x-grids, but is
       much slower for 1d cubic splines
    """
    from scipy.interpolate import splrep, splev
    from scipy.signal import cspline1d, cspline1d_eval
    if fast:
        return cspline1d_eval(cspline1d(oldy), newx, dx=oldx[1]-oldx[0],x0=oldx[0])
    else:
        rep = splrep(oldx,oldy,s=smoothing,full_output=False,**kw)
        return splev(newx, rep)

def sr_interpol3(x,y,ytarget,doplot=0,factor=10):
    s = numpy.sign(numpy.diff(y)[0])
    if s==1:
        y[numpy.argmax(y)+1:] = 2*abs(max(y))
    else:
        y[numpy.argmin(y)+1:] = -2*abs(max(y))
    dx = x[1]-x[0]
    newx = linspace(min(x),max(x),factor*len(x))

    cj = cspline1d(y)
    newy = cspline1d_eval(cj, newx, dx=dx,x0=x[0])

    ysq = (ytarget-newy)**2
    index = where(ysq == min(ysq))

    if doplot:
        clf()
        plot(x,y,'o')
        plot(newx,newy)
        plot(newx[index],newy[index],'o')
        show()

    return newx[index[0][0]]