def pntsInterpTendon(self,nPntsFine,smoothness,kgrade=3):
     '''Generates a cubic spline (default) or a spline of grade kgrade 
     interpolation from the rough points and calculates its value and 
     the value of its derivative in nPntsFine equispaced.
     Creates the following attributes:
     
     - fineCoordMtr: matrix with coordinates of the interpolated points
     [[x1,x2, ..],[y1,y2,..],[z1,z2,..]]
     - fineDerivMtr: matrix with the vector representing the derivative
     in each interpolated  
     - tck: tuple (t,c,k) containing the vector of knots, the B-spline 
            coefficients, and the degree of the spline.
     - fineScoord: curvilinear coordinate (cummulative length of curve 
                   for in each point)
     - fineProjXYcoord: coordinate by the projection of the curve on the XY 
     plane. Matrix 1*nPntsFine whose first elements is 0 and the rest the 
     cumulative distance to the first point
     '''
     tck, u = interpolate.splprep(self.roughCoordMtr, k=kgrade,s=smoothness)
     x_knots, y_knots,z_knots = interpolate.splev(np.linspace(0, 1, nPntsFine), tck,der=0)
     self.fineCoordMtr=np.array([x_knots, y_knots,z_knots])
     x_der, y_der,z_der = interpolate.splev(np.linspace(0, 1,nPntsFine), tck,der=1)
     self.fineDerivMtr=np.array([x_der, y_der,z_der])
     self.tck=tck
     self.fineScoord=self.getCumLength()
     x0,y0=(self.fineCoordMtr[0][0],self.fineCoordMtr[1][0])
     self.fineProjXYcoord=((self.fineCoordMtr[0]-x0)**2+(self.fineCoordMtr[1]-y0)**2)**0.5
     return

def testMVCgetDerivWpt(W):
    Ndim = W.shape[0]
    Nwaypoints = W.shape[1]
    dW = np.zeros((W.shape))
    ddW = np.zeros((W.shape))

    traj, tmp = splprep(W, k=5, s=0.01)
    # L = getLengthWpt(W)
    d = 0.0
    for i in range(0, Nwaypoints - 1):
        dW[:, i] = splev(d, traj, der=1)
        ddW[:, i] = splev(d, traj, der=2)
        # dW[:,i] = dW[:,i]/np.linalg.norm(dW[:,i])

        ds = np.linalg.norm(W[:, i + 1] - W[:, i])
        dv = np.linalg.norm(dW[:, i])
        dt = ds / dv
        # ddW[:,i] = ddW[:,i]/np.linalg.norm(ddW[:,i])
        print d
        d = d + dt

    dW[:, Nwaypoints - 1] = splev(d, traj, der=1)
    ddW[:, Nwaypoints - 1] = splev(d, traj, der=2)

    return [dW, ddW]

def interpolation_polynom(path,grad):
#    data=np.ndarray(shape=(len(path),3),dtype=float)   #create an array of float type for the input points
#    #fill the array with the Pathdata
#    a=path[0]
#    b=path[1]
#    c=path[2]
#    for i in range(len(a)):
#        data[i,0]=a[i]
#        data[i,1]=b[i]
#        data[i,2]=c[i]
#    #arrange the data to use the function
#    data = data.transpose()
    #interpolate polynom degree 1
    if grad==1:
        tck, u= interpolate.splprep(path,k=1,s=10)
        path = interpolate.splev(np.linspace(0,1,200), tck)
    #interpolate polynom degree 2
    if grad==2:
        tck, u= interpolate.splprep(path,k=2,s=10)
        path = interpolate.splev(np.linspace(0,1,200), tck)
    #interpolate polynom degree 3
    if grad==3:
        tck, u= interpolate.splprep(path, w=None, u=None, ub=None, ue=None, k=3, task=0, s=0.3, t=None, full_output=0, nest=None, per=0, quiet=1)
        path = interpolate.splev(np.linspace(0,1,200), tck)
    return path

    def romanzuniga07(wavelength, AKs, makePlot=False):
        filters = ['J', 'H', 'Ks', '[3.6]', '[4.5]', '[5.8]', '[8.0]']
        wave =      np.array([1.240, 1.664, 2.164, 3.545, 4.442, 5.675, 7.760])
        A_AKs =     np.array([2.299, 1.550, 1.000, 0.618, 0.525, 0.462, 0.455])
        A_AKs_err = np.array([0.530, 0.080, 0.000, 0.077, 0.063, 0.055, 0.059])
        
        # Interpolate over the curve
        spline_interp = interpolate.splrep(wave, A_AKs, k=3, s=0)

        A_AKs_at_wave = interpolate.splev(wavelength, spline_interp)
        A_at_wave = AKs * A_AKs_at_wave

        if makePlot:
            py.clf()
            py.errorbar(wave, A_AKs, yerr=A_AKs_err, fmt='bo', 
                        markerfacecolor='none', markeredgecolor='blue',
                        markeredgewidth=2)

            # Make an interpolated curve.
            wavePlot = np.arange(wave.min(), wave.max(), 0.1)
            extPlot = interpolate.splev(wavePlot, spline_interp)
            py.loglog(wavePlot, extPlot, 'k-')

            # Plot a marker for the computed value.
            py.plot(wavelength, A_AKs_at_wave, 'rs',
                    markerfacecolor='none', markeredgecolor='red',
                    markeredgewidth=2)
            py.xlabel('Wavelength (microns)')
            py.ylabel('Extinction (magnitudes)')
            py.title('Roman Zuniga et al. 2007')


        return A_at_wave

    def test_insert(self):
        b, b2, xx = self.b, self.b2, self.xx

        j = b.t.size // 2
        tn = 0.5*(b.t[j] + b.t[j+1])

        bn, tck_n = insert(tn, b), insert(tn, (b.t, b.c, b.k))
        assert_allclose(splev(xx, bn),
                        splev(xx, tck_n), atol=1e-15)
        assert_(isinstance(bn, BSpline))
        assert_(isinstance(tck_n, tuple))   # back-compat: tck in, tck out

        # for n-D array of coefficients, BSpline.c needs to be transposed
        # after that, the results are equivalent.
        sh = tuple(range(b2.c.ndim))
        c_ = b2.c.transpose(sh[1:] + (0,))
        tck_n2 = insert(tn, (b2.t, c_, b2.k))

        bn2 = insert(tn, b2)

        # need a transpose for comparing the results, cf test_splev
        assert_allclose(np.asarray(splev(xx, tck_n2)).transpose(2, 0, 1),
                        bn2(xx), atol=1e-15)
        assert_(isinstance(bn2, BSpline))
        assert_(isinstance(tck_n2, tuple))   # back-compat: tck in, tck out

def align_core(spks, subsmp=ALIGN_SUBSMP, maxdt=ALIGN_MAXDT,
        peakloc=ALIGN_PEAKLOC, findbwd=ALIGN_FINDBWD,
        findfwd=ALIGN_FINDFWD, cutat=ALIGN_CUTAT, outdim=ALIGN_OUTDIM,
        peakfunc=ALIGN_PEAKFUNC):
    """Alignment algorithm based on (Quiroga et al., 2004)"""

    if peakfunc == 'argmin':
        peakfunc = np.argmin
    else:
        raise ValueError('Not recognized "peakfunc"')

    R = np.empty((spks.shape[0], outdim), dtype='int16')
    n = spks.shape[1]
    x0 = np.arange(n)

    for i_spk, spk in enumerate(spks):
        tck = ipl.splrep(x0, spk, s=0)

        xn = np.arange(peakloc - findbwd, peakloc + findfwd, 1. / subsmp)
        yn = ipl.splev(xn, tck)

        dt = xn[peakfunc(yn)] - peakloc
        if np.abs(dt) > maxdt:
            dt = 0

        x = x0 + dt
        y = ipl.splev(x, tck)

        R[i_spk] = np.round(y).astype('int16')[cutat: cutat + outdim]
        #dts.append(dt)
    return R

def interpolate_1km_geolocation(lons_40km, lats_40km):
    """Interpolate AVHRR 40km navigation to 1km.

    This code was extracted from the python-geotiepoints package from the PyTroll group. To avoid adding another
    dependency to this package this simple case from the geotiepoints was copied.
    """
    cols40km = numpy.arange(24, 2048, 40)
    cols1km = numpy.arange(2048)
    lines = lons_40km.shape[0]
    # row_indices = rows40km = numpy.arange(lines)
    rows1km = numpy.arange(lines)

    lons_rad = numpy.radians(lons_40km)
    lats_rad = numpy.radians(lats_40km)
    x__ = EARTH_RADIUS * numpy.cos(lats_rad) * numpy.cos(lons_rad)
    y__ = EARTH_RADIUS * numpy.cos(lats_rad) * numpy.sin(lons_rad)
    z__ = EARTH_RADIUS * numpy.sin(lats_rad)
    along_track_order = 1
    cross_track_order = 3

    lines = len(rows1km)
    newx = numpy.empty((len(rows1km), len(cols1km)), x__.dtype)
    newy = numpy.empty((len(rows1km), len(cols1km)), y__.dtype)
    newz = numpy.empty((len(rows1km), len(cols1km)), z__.dtype)
    for cnt in range(lines):
        tck = splrep(cols40km, x__[cnt, :], k=cross_track_order, s=0)
        newx[cnt, :] = splev(cols1km, tck, der=0)
        tck = splrep(cols40km, y__[cnt, :], k=cross_track_order, s=0)
        newy[cnt, :] = splev(cols1km, tck, der=0)
        tck = splrep(cols40km, z__[cnt, :], k=cross_track_order, s=0)
        newz[cnt, :] = splev(cols1km, tck, der=0)

    lons_1km = get_lons_from_cartesian(newx, newy)
    lats_1km = get_lats_from_cartesian(newx, newy, newz)
    return lons_1km, lats_1km

def check_beta(beta):
    # now checking beta <= 1
    if max(beta)>1.:
        print('max beta!')              # TODO: remove, should be done by my_prior already
        return True
    # TODO: check smoothness of beta

    # now checking physical kappa: g(rvar, rfix, beta, dbetadr) >= 0
    if gp.usekappa == False:
        return False
    r0 = gp.xipol
    dR = r0[1:]-r0[:-1]
    r0extl = np.array([r0[0]/6., r0[0]/5., r0[0]/4., r0[0]/3., r0[0]/2., r0[0]/1.5])
    
    # extrapolation to the right (attention, could overshoot)
    dr0 = (r0[-1]-r0[-2])/8.
    r0extr = np.hstack([r0[-1]+dr0, r0[-1]+2*dr0, r0[-1]+3*dr0, r0[-1]+4*dr0])

    r0nu = np.hstack([r0extl, r0, dR/2.+r0[:-1], r0extr])
    r0nu.sort()
    tck0 = splrep(r0,beta*(r0**2+np.median(r0)**2), k=1, s=0.) # previous: k=2, s=0.1
    betanu = splev(r0nu,tck0)/(r0nu**2+np.median(r0)**2)

    drspl = splev(r0nu,tck0, der=1)
    dbetanudr = (drspl-betanu*2*r0nu)/(r0nu**2+np.median(r0)**2)
    for i in range(len(r0nu)-4):
        for j in range(i+1,len(r0nu)):
            if g(r0nu[j], r0nu[i], betanu[j], dbetanudr[j]) < 0:
                return True

    return False

def second_derivative(xdata, inds, gt=False, s=0):
    '''
    The second derivative of d^2 xdata / d inds^2

    why inds for interpolation, not log l?
    if not using something like model number instead of log l,
    the tmin will get hidden by data with t < tmin but different
    log l. This is only a problem for very low Z.
    If I find the arg min of teff to be very close to MS_BEG it
    probably means the MS_BEG is at a lower Teff than Tmin.
    '''
    tckp, _ = splprep([inds, xdata], s=s, k=3)
    arb_arr = np.arange(0, 1, 1e-2)
    xnew, ynew = splev(arb_arr, tckp)
    # second derivative, bitches.
    ddxnew, ddynew = splev(arb_arr, tckp, der=2)
    ddyddx = ddynew / ddxnew
    # not just argmin, but must be actual min...
    try:
        if gt:
            aind = [a for a in np.argsort(ddyddx) if ddyddx[a-1] < 0][0]
        else:
            aind = [a for a in np.argsort(ddyddx) if ddyddx[a-1] > 0][0]
    except IndexError:
        return -1
    tmin_ind, _ = closest_match2d(aind, inds, xdata, xnew, ynew)
    return inds[tmin_ind]

def nishiyama09(wavelength, AKs, makePlot=False):
    # Data pulled from Nishiyama et al. 2009, Table 1

    filters = ['V', 'J', 'H', 'Ks', '[3.6]', '[4.5]', '[5.8]', '[8.0]']
    wave =      np.array([0.551, 1.25, 1.63, 2.14, 3.545, 4.442, 5.675, 7.760])
    A_AKs =     np.array([16.13, 3.02, 1.73, 1.00, 0.500, 0.390, 0.360, 0.430])
    A_AKs_err = np.array([0.04,  0.04, 0.03, 0.00, 0.010, 0.010, 0.010, 0.010])

    # Interpolate over the curve
    spline_interp = interpolate.splrep(wave, A_AKs, k=3, s=0)

    A_AKs_at_wave = interpolate.splev(wavelength, spline_interp)
    A_at_wave = AKs * A_AKs_at_wave

    if makePlot:
        py.clf()
        py.errorbar(wave, A_AKs, yerr=A_AKs_err, fmt='bo', 
                    markerfacecolor='none', markeredgecolor='blue',
                    markeredgewidth=2)
        
        # Make an interpolated curve.
        wavePlot = np.arange(wave.min(), wave.max(), 0.1)
        extPlot = interpolate.splev(wavePlot, spline_interp)
        py.loglog(wavePlot, extPlot, 'k-')

        # Plot a marker for the computed value.
        py.plot(wavelength, A_AKs_at_wave, 'rs',
                markerfacecolor='none', markeredgecolor='red',
                markeredgewidth=2)
        py.xlabel('Wavelength (microns)')
        py.ylabel('Extinction (magnitudes)')
        py.title('Nishiyama et al. 2009')

    
    return A_at_wave

def expand_traj_dim_with_derivative(data, dt=0.01):
    augmented_trajs = []
    for traj in data:
        time_len = len(traj)
        t = np.linspace(0, time_len * dt, time_len)
        if time_len > 3:
            if len(traj.shape) == 1:
                """
                mono-dimension trajectory, row as the entire trajectory...
                """
                spl = interpolate.splrep(t, traj)
                traj_der = interpolate.splev(t, spl, der=1)
                tmp_augmented_traj = np.array([traj, traj_der]).T
            else:
                """
                multi-dimensional trajectory, row as the state variable...
                """
                tmp_traj_der = []
                for traj_dof in traj.T:
                    spl_dof = interpolate.splrep(t, traj_dof)
                    traj_dof_der = interpolate.splev(t, spl_dof, der=1)
                    tmp_traj_der.append(traj_dof_der)
                tmp_augmented_traj = np.vstack([traj.T, np.array(tmp_traj_der)]).T

            augmented_trajs.append(tmp_augmented_traj)

    return augmented_trajs

def interp_data_fixed_num(data_set, num=100):
    """
    interpolate data with fixed number of points
    """
    interp_data = dict()
    for key in data_set:
        interp_data[key] = []
        for l in data_set[key]:
            interp_letter = []
            for s in l:
                time_len = s.shape[0]
                if time_len > 3:
                    #interpolate each dim, cubic
                    t = np.linspace(0, 1, time_len)
                    spl_x = interpolate.splrep(t, s[:, 0])
                    spl_y = interpolate.splrep(t, s[:, 1])
                    #resample, 4 times more, vel is also scaled...
                    t_spl = np.linspace(0, 1, num)
                    x_interp = interpolate.splev(t_spl, spl_x, der=0)
                    y_interp = interpolate.splev(t_spl, spl_y, der=0)
                    # #construct new stroke
                    data = np.concatenate([x_interp, y_interp])
                    dt = float(time_len)/num
                    interp_letter.append(np.concatenate([data, [dt]]))
                else:
                    #direct copy if no sufficient number of points
                    interp_letter.append(s)
            interp_data[key].append(interp_letter)
    return interp_data

def interp_data(data_set):
    """
    interpolate data
    """
    interp_data = dict()
    for key in data_set:
        interp_data[key] = []
        for l in data_set[key]:
            interp_letter = []
            for s in l:
                time_len = s.shape[0]
                if time_len > 3:
                    #interpolate each dim, cubic
                    t = np.linspace(0, 1, time_len)
                    spl_x = interpolate.splrep(t, s[:, 0])
                    spl_y = interpolate.splrep(t, s[:, 1])
                    #resample, 4 times more, vel is also scaled...
                    t_spl = np.linspace(0, 1, 4 * len(t))
                    x_interp = interpolate.splev(t_spl, spl_x, der=0)
                    y_interp = interpolate.splev(t_spl, spl_y, der=0)
                    # #construct new stroke
                    interp_letter.append(np.concatenate([[x_interp], [y_interp]], axis=0).transpose())
                else:
                    #direct copy if no sufficient number of points
                    interp_letter.append(s)
            interp_data[key].append(interp_letter)
    return interp_data

def Interpo(spectra) :
    wave_min = 1000
    wave_max = 20000
    pix = 2
    #wavelength = np.linspace(wave_min,wave_max,(wave_max-wave_min)/pix+1)  #creates N equally spaced wavelength values
    wavelength = np.arange(ceil(wave_min), floor(wave_max), dtype=int, step=pix)
    fitted_flux = []
    fitted_error = []
    new = []
    #new = Table()
    #new['col0'] = Column(wavelength,name = 'wavelength')
    new_spectrum=spectra	#declares new spectrum from list
    new_wave=new_spectrum[:,0]	#wavelengths
    new_flux=new_spectrum[:,1]	#fluxes
    new_error=new_spectrum[:,2]   #errors
    lower = new_wave[0] # Find the area where interpolation is valid
    upper = new_wave[len(new_wave)-1]
    lines = np.where((new_wave>lower) & (new_wave<upper))	#creates an array of wavelength values between minimum and maximum wavelengths from new spectrum
    indata=inter.splrep(new_wave[lines],new_flux[lines])	#creates b-spline from new spectrum
    inerror=inter.splrep(new_wave[lines],new_error[lines]) # doing the same with the errors
    fitted_flux=inter.splev(wavelength,indata)	#fits b-spline over wavelength range
    fitted_error=inter.splev(wavelength,inerror)   # doing the same with errors
    badlines = np.where((wavelength<lower) | (wavelength>upper))
    fitted_flux[badlines] = 0  # set the bad values to ZERO !!! 
    new = Table([wavelength,fitted_flux],names=('col1','col2')) # put the interpolated data into the new table    
    #newcol = Column(fitted_flux,name = 'Flux')  
    #new.add_column(newcol,index = None)
    return new

    def dX_dv(self, y, v, other_variables, classical_source, solver):

        alpha = self.alpha
        r = y[0]
        sigma = y[1]
        f = y[2]
        g = y[3]
        h = y[4]
        I = y[5]
        J = y[6]

        u = other_variables[0]
        k = other_variables[1]
        F = splev(v, self.F_tck)
        d2A_du2 = splev(u, self.A_tck, 2)
        dA_du = splev(u, self.A_tck, 1)

        df_du = 2*f*I - alpha/r*(J - I**2. + d2A_du2 + dA_du**2.)

        dg_du = -(r/(r**2.- alpha))*(f*g + exp(2*sigma)/4.)

        dr_dv = g

        dsigma_dv = h

        df_dv = dg_du

        dg_dv = 2*g*h - F/r - alpha/r*(k - h**2.)

        dh_dv = k

        dI_dv = (1./(r**2.-alpha))*(f*g + np.exp(2*sigma)/4.) # = dH_du

        dJ_dv =   (-2.*f/(r**2.-alpha)**2.)*(f*g + np.exp(2*sigma)/4.) 
                + (1./(r**2.-alpha))*(df_du*g + f*dg_du + np.exp(2*sigma)*I/2.)

def ladybug_interp_data( data ):
    '''Interpolates between each 'new' GPS coordinate. Repeated coordinates are
    assumed to be incorrect/inaccurate. This method replaces all repetitions
    with interpolated coordinates in place. This method uses cubic spline
    interpolation.
    '''
    for i in reversed(xrange(1,len(data['lon']))):
        if data['lon'][i] == data['lon'][i-1]:
            data['lon'][i] = 1000.0
            data['valid'][i] = False
    select = where(data['lon'] < 999)

    # SPLINE VERSION
    data['alt'] = interpolate.splev(data['seqid'],
                                    interpolate.splrep(data['seqid'][select],
                                                       data['alt'][select],
                                                       s=0, k=2  ),
                                    der=0)
    data['lon'] = interpolate.splev(data['seqid'],
                                    interpolate.splrep(data['seqid'][select],
                                                       data['lon'][select],
                                                       s=0, k=2  ),
                                    der=0)
    data['lat'] = interpolate.splev(data['seqid'],
                                    interpolate.splrep(data['seqid'][select],
                                                       data['lat'][select],
                                                       s=0, k=2  ),
                                    der=0)
    return data

    def dX_du(self, y, u, other_variables, solver):

        alpha = self.alpha
        r = y[0]
        sigma = y[1]
        f = y[2]
        g = y[3]
        h = y[4]
        I = y[5]
        k = y[6]

        v = other_variables[0]
        J = other_variables[1]
        F = splev(v, self.F_tck)
        d2A_du2 = splev(u, self.A_tck, 2)
        dA_du = splev(u, self.A_tck, 1)

        dr_du = f

        dsigma_du = I

        df_du = 2*f*I - alpha/r*(J - I**2. + d2A_du2 + dA_du**2.)

        dg_du = -(r/(r**2.- alpha))*(f*g + exp(2*sigma)/4.)

        dh_du = (1./(r**2.-alpha))*(f*g + np.exp(2*sigma)/4.) 

        dI_du = J

        df_dv = dg_du

        dg_dv = 2*g*h - F/r - alpha/r*(k - h**2.)

        dk_du =   (-2.*g/(r**2.-alpha)**2.)*(f*g + np.exp(2*sigma)/4.) 
                + (1./(r**2.-alpha))*(df_dv*g + f*dg_dv + np.exp(2*sigma)*h/2.)

    def _interp1d(self):
        """Interpolate in one dimension."""
        lines = len(self.hrow_indices)

        self.newx = np.empty((len(self.hrow_indices),
                              len(self.hcol_indices)),
                             self.x__.dtype)

        self.newy = np.empty((len(self.hrow_indices),
                              len(self.hcol_indices)),
                             self.y__.dtype)

        self.newz = np.empty((len(self.hrow_indices),
                              len(self.hcol_indices)),
                             self.z__.dtype)

        for cnt in range(lines):
            tck = splrep(self.col_indices, self.x__[cnt, :], k=self.ky_, s=0)
            self.newx[cnt, :] = splev(self.hcol_indices, tck, der=0)

            tck = splrep(self.col_indices, self.y__[cnt, :], k=self.ky_, s=0)
            self.newy[cnt, :] = splev(self.hcol_indices, tck, der=0)

            tck = splrep(self.col_indices, self.z__[cnt, :], k=self.ky_, s=0)
            self.newz[cnt, :] = splev(self.hcol_indices, tck, der=0)

   def __call__(self, x):
      '''Interpolate at point [x].  Returns a 3-tuple: (y, mask) where [y]
      is the interpolated point, and [mask] is a boolean array with the same
      shape as [x] and is True where interpolated and False where extrapolated'''
      if not self.setup:
         self._setup()

      if len(num.shape(x)) < 1:
         scalar = True
      else:
         scalar = False

      x = num.atleast_1d(x)
      if self.realization:
         evm = num.atleast_1d(splev(x, self.realization))
         mask = num.greater_equal(x, self.realization[0][0])*\
                num.less_equal(x,self.realization[0][-1])
      else:
         evm = num.atleast_1d(splev(x, self.tck))
         mask = num.greater_equal(x, self.tck[0][0])*num.less_equal(x,self.tck[0][-1])

      if scalar:
         return evm[0],mask[0]
      else:
         return evm,mask

 def __init__(self, freq,
              nripples=2,
              **keywords):
     nripples_max = 2
     if nripples not in range(nripples_max + 1):
         raise ValueError(
             'Input nripples is not a non-negative integer less than {}'.
             format(nripples_max + 1))
     self.nripples = nripples
     with open(PATH + 'sb_peak_plus_two_ripples_150HGz.pkl', 'r') as f:
         fl = load(f)
     fl /= fl.max()
     if freq == 150e9:
         fl_ = fl
     else:
         corr1 = [  1.65327594e-02, -2.24216210e-04, 9.70939946e-07, -1.40191824e-09]
         corr2 = [ -3.80559542e-01, 4.76370274e-03, -1.84237511e-05, 2.37962542e-08]
         def f(x, p):
             return p[0] + p[1] * x + p[2] * x**2 + p[3] * x**3
         ell = np.arange(len(fl)) + 1
         spl = splrep(ell * freq / 150e9, fl)
         if freq > 150e9:
             fl_ = splev(ell, spl) * (1 + f(freq / 1e9, corr2) + ell * f(freq / 1e9, corr1))
         else:
             fl_ = np.zeros(len(ell))
             fl_ = splev(ell[ell < ell.max() * freq / 150e9], spl) * \
                    (1 + f(freq / 1e9, corr2) + ell[ell < ell.max() * freq / 150e9] * f(freq / 1e9, corr1))
     self.fl = np.sqrt(fl_)
     self.fl[np.isnan(self.fl)] = 0.