Python scipy.interpolate 模块，splrep() 实例源码

def compute_beta_splines(TT, minT, splinesoptions={}):
    """
    This function computes the volumetric thermal expansion as a numerical
    derivative of the volume as a function of temperature V(T) given in the
    input array *minT*. This array can obtained
    from the free energy minimization which should be done before.
    This function uses splines for smoothing the derivative.
    """

    betaT = np.zeros(len(TT))

    x = np.array(TT)
    y0 = np.array(minT)

    if (splinesoptions == {}):
        tck0 = interpolate.splrep(x, y0)
    else:
        tck0 = interpolate.splrep(x, y0, k=splinesoptions['k0'], s=splinesoptions['s0'])

    ynew0 = interpolate.splev(x, tck0, der=0)
    betaT = interpolate.splev(x, tck0, der=1)

    return betaT/minT

def resample_photosphere(opacities, photosphere, opacity_index):
    """ Resample photospheric quantities onto a new opacity scale. """

    if opacity_index is None:
        return photosphere

    resampled_photosphere = np.zeros(photosphere.shape)
    n_quantities = photosphere.shape[1]
    for i in range(n_quantities):
        if i == opacity_index: continue
        # Create spline function.
        tk = \
            interpolate.splrep(photosphere[:, opacity_index], photosphere[:, i])

        # Evaluate photospheric quantities at the new opacities
        resampled_photosphere[:, i] = interpolate.splev(opacities.flatten(), tk)

    # Update photosphere with new opacities
    resampled_photosphere[:, opacity_index] = opacities
    return resampled_photosphere

def resample_1d(data, npoints, smooth = 0, periodic = False, derivative = 0):
  """Resample 1d data using n equidistant points

  Arguments:
    data (array): data points
    npoints (int): number of points in equidistant resampling
    smooth (number): smoothness factor
    periodic (bool): if True assumes the curve is a closed curve

  Returns:
    (array): resampled data points
  """

  x = np.linspace(0, 1, data.shape[0]);
  dinterp = splrep(x, data, s = smooth, per = periodic);
  if npoints is all:
    npoints = data.shape[0];
  x2 = np.linspace(0, 1, npoints);
  return splev(x2, dinterp, der = derivative)

def phi(self, points = None):
    """Returns a Spline representing the tangent angle along a 2d curve

    Arguments:
      points (int, array or None): sample points used to determine the phi

    Returns:
      Spline: spline of phi
    """
    if self.ndim != 2:
      raise RuntimeError('phi angle can only be computed for 2d curves');

    points = self.get_points(points, error = 'cannot determine sample points needed for the calculation of phi');    

    #get the tangents  and tangent angles
    tgs = splev(points, self.tck(), der = 1);
    tgs = np.vstack(tgs).T;
    phi = np.arctan2(tgs[:,1], tgs[:,0]);
    phi = np.mod(phi + np.pi, 2 * np.pi) - np.pi;
    #return Spline(phi, points = self.points, knots = self.knots, degree = self.degree - 1);
    tck = splrep(points, phi, s = 0.0, k = self.degree + 1);
    return Spline(tck = tck, points = points);

def set_parameter_from_values(self, values, points = None):
    """Change the parameter to represent the values and resmaple on current sample points if necessary

    Arguments: 
      values (array): values
      points (array or None): sample points for the values
    """

    if points is None:
      self.set_values(values);
    else:
      if points is self._points or (self._points is not None and self._points.shape[0] == points.shape[0] and np.allclose(points, self._points)):
        self.set_values(values);
      else:
        tck = splrep(points, values, t = self._knots, task = -1, k = self.degree);
        self._parameter = tck[1][:self._nparameter];
        # values changed
        self._values = None;

def set_parameter_from_values(self, values, points = None):
    """Change the parameter to represent the values and resmaple on current sample points if necessary

    Arguments: 
      values (array): values
      points (array or None): sample points for the values
    """

    if points is None:
      self.set_values(values);
    else:
      if points is self._points or (self._points is not None and self._points.shape[0] == points.shape[0] and np.allclose(points, self._points)):
        self.set_values(values);
      else:
        tck = splrep(points, values, t = self._knots, task = -1, k = self.degree);
        self._parameter = tck[1][:self._nparameter];
        # values changed
        self._values = None;

def optc(self):
        """
        Optimize the coeffs, don't touch the knots
        This is the fast guy, one reason to use splines :-)
        Returns the chi2 in case you want it (including stabilization points) !

        Sets lastr2stab, but not lastr2nostab !

        """

        out = si.splrep(self.datapoints.jds, self.datapoints.mags, w=1.0/self.datapoints.magerrs, xb=None, xe=None, k=self.k, task=-1, s=None, t=self.getintt(), full_output=1, per=0, quiet=1)
        # We check if it worked :
        if not out[2] <= 0: 
            raise RuntimeError("Problem with spline representation, message = %s" % (out[3]))

        self.c = out[0][1] # save the coeffs

        #import matplotlib.pyplot as plt
        #plt.plot(self.datapoints.jds, self.datapoints.magerrs)
        #plt.show()

        self.lastr2stab = out[1]
        return out[1]

def klgfbInterp(self, wvl, n, l):
        '''A Python version of the CHIANTI IDL procedure karzas_xs.

        Interpolates free-bound gaunt factor of Karzas and Latter, (1961, Astrophysical Journal
        Supplement Series, 6, 167) as a function of wavelength (wvl).
        '''
        try:
            klgfb = self.Klgfb
        except:
            self.Klgfb = ch_io.klgfbRead()
            klgfb = self.Klgfb
        # get log of photon energy relative to the ionization potential
        sclE = np.log(self.Ip/(wvl*ch_const.ev2ang))
        thisGf = klgfb['klgfb'][n-1, l]
        spl = splrep(klgfb['pe'], thisGf)
        gf = splev(sclE, spl)
        return gf

def ioneq_one(self, stage, **kwargs):
        """
        Calculate the equilibrium fractional ionization of the ion as a function of temperature.

        Uses the `ChiantiPy.core.ioneq` module and does a first-order spline interpolation to the data. An
        ionization equilibrium file can be passed as a keyword argument, `ioneqfile`. This can
        be passed through as a keyword argument to any of the functions that uses the
        ionization equilibrium.

        Parameters
        ----------
        stage : int
            Ionization stage, e.g. 25 for Fe XXV
        """
        tmp = ioneq(self.Z)
        tmp.load(ioneqName=kwargs.get('ioneqfile', None))
        ionization_equilibrium = splev(self.Temperature,
                                       splrep(tmp.Temperature, tmp.Ioneq[stage-1,:], k=1), ext=1)
        return np.where(ionization_equilibrium < 0., 0., ionization_equilibrium)

def setDriftData(self, alpha, CL, CDi, CM):
        self.alpha = alpha
        self.CL    = CL
        self.CDi   = CDi
        self.CM    = CM

        # Calculate coefficients for lift
        fitParams, fitCovariances = curve_fit(liftFunc, self.alpha, self.CL)

        self.CL_a1 = fitParams[0]
        self.CL_a2 = fitParams[1]

        # Calculate coefficients for induced drag
        fitParams, fitCovariances = curve_fit(inducedDragFunc, self.alpha, self.CDi)

        self.CDi_a2 = fitParams[0]

        # Spline interpolation for moment
        self.CMspl = interpolate.splrep(self.alpha, self.CM, s=1e-9)

def setStabilityData(self, heelAngles, GZ):
        self.heelAngles = heelAngles
        self.GZ         = GZ

        self.restoringMoment = self.GZ*self.Volume*self.rho*self.g

        self.restoringMomentSpl = interpolate.splrep(self.heelAngles, self.restoringMoment)

        # Calculate maximum heel angle
        nTest = 100
        heelAnglesTest = np.linspace(0, np.max(self.heelAngles), nTest)
        restoringMomentTest = interpolate.splev(heelAnglesTest, self.restoringMomentSpl)

        iMax = np.argmax(restoringMomentTest)

        self.maxHeelAngle = heelAnglesTest[iMax]

def alt_sky_model(self):
        fac = 10
        mn = 1e9
        mx = 0.
        for fiber in self.good_fibers:
            mn = np.min([mn, fiber.wavelength.min()])
            mx = np.max([mx, fiber.wavelength.max()])
        xs = np.linspace(0, 1, self.D*fac)
        A = np.zeros((len(xs), len(self.good_fibers)))
        for i, fiber in enumerate(self.good_fibers):
            y = fiber.spectrum / fiber.fiber_to_fiber
            xp = np.interp(fiber.wavelength, np.linspace(mn, mx, self.D*fac),
                           xs, left=0.0, right=0.0)
            tck = splrep(xp, y)
            A[:, i] = splev(xs, tck)
        ys = biweight_location(A, axis=(1,))
        self.masterwave = np.linspace(mn, mx, self.D*fac)
        B, c = bspline_x0(self.masterwave, nknots=self.D)
        sol = np.linalg.lstsq(c, ys)[0]
        self.mastersky = np.dot(c, sol)

def derivatives(xs,ys,ders=(1,)):
    '''
    Calculate the numerical derivatives of `ys` at points `xs` by use of the third-order spline interpolation.

    Parameters
    ----------
    xs : 1d ndarray
        The sample points of the argument.
    ys: 1d ndarray
        The corresponding sample points of the function.
    ders : tuple of integer
        The derivatives to calculate.

    Returns
    -------
    2d ndarray
        The derivatives at the sample points of the argument.
    '''
    assert len(xs)==len(ys)
    xs,ys=np.asarray(xs),np.asarray(ys)
    result=np.zeros((len(ders),len(xs)),dtype=ys.dtype)
    tck=ip.splrep(xs,ys,s=0)
    for i,der in enumerate(ders):
        result[i]=ip.splev(xs,tck,der=der)
    return result

def _setup_channels(self):
        """
        Setup channel, specifically the wavelength or temperature response functions.

        Notes
        -----
        This should be replaced once the response functions are available in SunPy. Probably should
        configure wavelength response function interpolators also.
        """
        aia_fn = pkg_resources.resource_filename('synthesizAR', 'instruments/data/sdo_aia.json')
        with open(aia_fn, 'r') as f:
            aia_info = json.load(f)

        for channel in self.channels:
            channel['name'] = str(channel['wavelength'].value).strip('.0')
            channel['instrument_label'] = '{}_{}'.format(self.fits_template['detector'],
                                                         channel['telescope_number'])
            channel['wavelength_range'] = None
            x = aia_info[channel['name']]['temperature_response_x']
            y = aia_info[channel['name']]['temperature_response_y']
            channel['temperature_response_spline'] = splrep(x, y)
            x = aia_info[channel['name']]['response_x']
            y = aia_info[channel['name']]['response_y']
            channel['wavelength_response_spline'] = splrep(x, y)

def lininterp2(x1, y1, x):
    """Linear interpolation at points x between numpy arrays (x1, y1).
    Only y1 is allowed to be two-dimensional.  The x1 values should be sorted
    from low to high.  Returns a numpy.array of y values corresponding to
    points x.
    """
    return splev(x, splrep(x1, y1, s=0, k=1))

def compute_alpha_splines(TT,minT,ibrav,splinesoptions):
    """
    This function calculates the thermal expansions alphaT at different temperatures
    as the previous function but using spline interpolation as implemented in
    scipy.interpolate.

    """
    alphaT = np.zeros(len(TT)*6)
    alphaT.shape = (len(TT),6)

    x = np.array(TT)
    y0 = np.array(minT[:,0])
    y1 = np.array(minT[:,1])
    y2 = np.array(minT[:,2])

    if (splinesoptions=={}):
        tck0 = interpolate.splrep(x, y0)
        tck1 = interpolate.splrep(x, y1)
        tck2 = interpolate.splrep(x, y2)  
    else:
        tck0 = interpolate.splrep(x, y0, k=splinesoptions['k0'], s=splinesoptions['s0'])
        tck1 = interpolate.splrep(x, y1, k=splinesoptions['k1'], s=splinesoptions['s1'])
        tck2 = interpolate.splrep(x, y2, k=splinesoptions['k2'], s=splinesoptions['s2'])        

    ynew0 = interpolate.splev(x, tck0, der=0)
    alphaT[:,0] = interpolate.splev(x, tck0, der=1)
    ynew1 = interpolate.splev(x, tck1, der=0)
    alphaT[:,1] = interpolate.splev(x, tck1, der=1)
    ynew2 = interpolate.splev(x, tck2, der=0)
    alphaT[:,2] = interpolate.splev(x, tck2, der=1)

    # now normalize the alphaTs properly. It must be different for different ibrav
    # to avoid a divide by 0 error (minT is zero for lattice parameters not defined
    # in the system)
    if ibrav==4:
        alphaT[:,0] = alphaT[:,0]/minT[:,0]
        alphaT[:,2] = alphaT[:,2]/minT[:,2]

    return alphaT

################################################################################

def compute_blackbody_response(self, Teffs=None):
        """
        Computes blackbody intensities across the entire range of
        effective temperatures.

        @Teffs: an array of effective temperatures. If None, a default
        array from ~300K to ~500000K with 97 steps is used. The default
        array is uniform in log10 scale.

        Returns: n/a
        """

        if Teffs == None:
            log10Teffs = np.linspace(2.5, 5.7, 97) # this corresponds to the 316K-501187K range.
            Teffs = 10**log10Teffs

        # Energy-weighted intensities:
        log10ints_energy = np.array([np.log10(self._bb_intensity(Teff, photon_weighted=False)) for Teff in Teffs])
        self._bb_func_energy = interpolate.splrep(Teffs, log10ints_energy, s=0)
        self._log10_Inorm_bb_energy = lambda Teff: interpolate.splev(Teff, self._bb_func_energy)

        # Photon-weighted intensities:
        log10ints_photon = np.array([np.log10(self._bb_intensity(Teff, photon_weighted=True)) for Teff in Teffs])
        self._bb_func_photon = interpolate.splrep(Teffs, log10ints_photon, s=0)
        self._log10_Inorm_bb_photon = lambda Teff: interpolate.splev(Teff, self._bb_func_photon)

        self.content.append('blackbody')
        self.atmlist.append('blackbody')

def burgess_tully_descale(x, y, energy_ratio, c, scaling_type):
        """
        Convert scaled Burgess-Tully parameters to physical quantities. For more details see
        [1]_.

        References
        ----------
        .. [1] Burgess, A. and Tully, J. A., 1992, A&A, `254, 436 <http://adsabs.harvard.edu/abs/1992A%26A...254..436B>`_ 
        """
        nots = splrep(x, y, s=0)
        if scaling_type == 1:
            x_new = 1.0 - np.log(c)/np.log(energy_ratio + c)
            upsilon = splev(x_new, nots, der=0)*np.log(energy_ratio + np.e)
        elif scaling_type == 2:
            x_new = energy_ratio/(energy_ratio + c)
            upsilon = splev(x_new, nots, der=0)
        elif scaling_type == 3:
            x_new = energy_ratio/(energy_ratio + c)
            upsilon = splev(x_new, nots, der=0)/(energy_ratio + 1.0)
        elif scaling_type == 4:
            x_new = 1.0 - np.log(c)/np.log(energy_ratio + c)
            upsilon = splev(x_new, nots, der=0)*np.log(energy_ratio + c)
        elif scaling_type == 5:
            # dielectronic
            x_new = energy_ratio/(energy_ratio + c)
            upsilon = splev(x_new, nots, der=0)/energy_ratio
        elif scaling_type == 6:
            # protons
            x_new = energy_ratio/(energy_ratio + c)
            upsilon = 10**splev(x_new, nots, der=0)
        else:
            raise ValueError('Unrecognized BT92 scaling option.')

        return upsilon

def from_tck(self, tck, points = None):
    """Change spline parameter and knot structure using a tck object returned by splprep or splrep

    Arguments:
        tck (tuple): t,c,k tuple returned by splprep
        points (int, array or None): optional sample points pecification
    """
    t,c,k = tck;
    self.degree = k;

    # set knots
    self._knots = t[k+1:-k-1];  
    self._knots_all = t;
    self._nparameter = self._knots.shape[0] + k + 1;

    #set parameter
    if isinstance(c, list):
      c = np.vstack(c);
    elif isinstance(c, np.ndarray) and c.ndim == 1:
      c = np.vstack([c]);
    c = c[:,:self._nparameter].T;

    self.ndim = c.shape[1];    
    self._parameter = c[:];

    #set points    
    if points is not None:
      self.set_points(points);

    # values and projection matrix will need to be updated after change of the knots
    self._values = None;
    self._projection = None;
    self._projection_inverse = None;  


  ############################################################################
  ### Spline getter

def tck_1d(self, dimension = 0):
    """Returns tck tuple for use with 1d spline functions

    Arguments:
      dimension (int): dimension for which to return tck tuple

    Returns:
      tuple: tck couple as retyurned by splrep
    """
    k = self.degree;
    p = np.pad(self._parameter[:,dimension], (0,k+1), 'constant');
    return (self._knots_all, p, k);

def set_values(self, values, points = None):
    """Calculate the bspline parameter for the given values and sample points 

    Arguments:
      values (array): values of data points
      points (array or None): sample points for the data values, if None use internal points or linspace(0,1,values.shape[0])
    """
    values = np.array(values, dtype = float);

    #set points
    if points is None:
      if self._points is None:
        self.set_points(values.shape[0]);
    else:
      self.set_points(points);

    if values.shape[0] != self._points.shape[0]:
      raise ValueError('number of values %d mismatch number of points %d' % (values.shape[0], self._points.shape[0]));

    #set parameter from values
    if self.with_projection and self.projection_inverse is not False:
      self._parameter = self._projection_inverse.dot(values); 
      self._values = values;  
    else:
      tck = splrep(self._points, values, t = self._knots, task = -1, k = self.degree);
      self._parameter = tck[1][:self._nparameter];
      # values changed
      self._values = None;
      #self._values = values; #fast but imprecise as values on spline due approximation might differ!

def from_tck(self, tck, points = None):
    """Set spline parameter and knot structure using a tck object returned by splrep

    Arguments:
        tck (tuple): t,c,k tuple as returned by splrep
        points (int, array or None): optional sample points pecification
    """
    t,c,k = tck;
    self.degree = k;

    # set knots
    self._knots = t[k+1:-k-1];  
    self._knots_all = t;
    self._nparameter = self._knots.shape[0] + k + 1;

    # set parameter
    self._parameter = c[:self._nparameter];

    #set points
    if points is not None:
      self.set_points(points);

    # values and projection matrix will need to be updated after change of the knots
    self._values = None;
    self._projection = None;
    self._projection_inverse = None;  


  ############################################################################
  ### Spline getter

def tck(self):
    """Returns a tck tuple for use with spline functions

    Note:
      This returns a tck tuple as splrep
    """  
    k = self.degree;
    p = np.pad(self._parameter, (0,k+1), 'constant');
    return (self._knots_all, p, k);

def set_values(self, values, points = None):
    """Calculate the bspline parameter for the given values and sample points 

    Arguments:
      values (array): values of data points
      points (array or None): sample points for the data values, if None use internal points or linspace(0,1,values.shape[0])
    """
    values = np.array(values, dtype = float);

    #set points
    if points is None:
      if self._points is None:
        self.set_points(values.shape[0]);
    else:
      self.set_points(points);

    if values.shape[0] != self._points.shape[0]:
      raise ValueError('number of values %d mismatch number of points %d' % (values.shape[0], self._points.shape[0]));

    #set parameter from values
    if self.with_projection and self.projection_inverse is not False:
      self._parameter = self._projection_inverse.dot(values); 
      self._values = values;  
    else:
      tck = splrep(self._points, values, t = self._knots, task = -1, k = self.degree);
      self._parameter = tck[1][:self._nparameter];
      # values changed
      self._values = None;
      #self._values = values; #fast but imprecise as values on spline due approximation might differ!

def from_tck(self, tck, points = None):
    """Set spline parameter and knot structure using a tck object returned by splrep

    Arguments:
        tck (tuple): t,c,k tuple as returned by splrep
        points (int, array or None): optional sample points pecification
    """
    t,c,k = tck;
    self.degree = k;

    # set knots
    self._knots = t[k+1:-k-1];  
    self._knots_all = t;
    self._nparameter = self._knots.shape[0] + k + 1;

    # set parameter
    self._parameter = c[:self._nparameter];

    #set points
    if points is not None:
      self.set_points(points);

    # values and projection matrix will need to be updated after change of the knots
    self._values = None;
    self._projection = None;
    self._projection_inverse = None;  


  ############################################################################
  ### Spline getter

def tck(self):
    """Returns a tck tuple for use with spline functions

    Note:
      This returns a tck tuple as splrep
    """  
    k = self.degree;
    p = np.pad(self._parameter, (0,k+1), 'constant');
    return (self._knots_all, p, k);

def convex_hull_removal(pixel, wvl):
    """
    Remove the convex-hull of the signal by hull quotient.

    Parameters:
        pixel: `list`
            1D HSI data (p), a pixel.
        wvl: `list`
            Wavelength of each band (p x 1).

    Results: `list`
        Data with convex hull removed (p).

    Reference:
        Clark, R.N. and T.L. Roush (1984) Reflectance Spectroscopy: Quantitative
        Analysis Techniques for Remote Sensing Applications, J. Geophys. Res., 89,
        6329-6340.
    """
    points = list(zip(wvl, pixel))
    # close the polygone
    poly = [(points[0][0],0)]+points+[(points[-1][0],0)]
    hull = _jarvis.convex_hull(poly)
    # the last two points are on the x axis, remove it
    hull = hull[:-2]
    x_hull = [u for u,v in hull]
    y_hull = [v for u,v in hull]

    tck = interpolate.splrep(x_hull, y_hull, k=1)
    iy_hull = interpolate.splev(wvl, tck, der=0)

    norm = []
    for ysig, yhull in zip(pixel, iy_hull):
        if yhull != 0:
            norm.append(ysig/yhull)
        else:
            norm.append(1)

    return norm, x_hull, y_hull

def smooth_plot(figure,X,Y,color1,color2,xlabel="",ylabel="",filled=False,n_points=400,smooth_factor=1.0,spline_order=3,linewidth=3,alpha=1.0,label=""):
    """
    """
    X_smooth = np.linspace(X.min(),X.max(),n_points)
    tck = splrep(X,Y,s=smooth_factor,k=spline_order)
    Y_smooth = splev(X_smooth,tck,der=0)

    font = fm.FontProperties(family = 'Trebuchet', weight ='light')
    #font = fm.FontProperties(family = 'CenturyGothic',fname = '/Library/Fonts/Microsoft/Century Gothic', weight ='light')
    figure.patch.set_facecolor('white')
    axes = figure.add_subplot(111)
    axes.plot(X,Y,linewidth=1,color=tuple(color2),alpha=0.2)
    if filled:
        axes.fill_between(X_smooth,Y_smooth,0,color=color2,alpha=0.1)
    for i in xrange(100):
        color = tuple(color1*(1.0-i/100.0) + color2*(i/100.0))
        if i == 0:
            axes.plot(X_smooth[(i*n_points/100):((i+1)*n_points)/100+1],Y_smooth[(i*n_points/100):((i+1)*n_points)/100+1],linewidth=linewidth,color=color,alpha=alpha,label=label)
        else:
            axes.plot(X_smooth[(i*n_points/100):((i+1)*n_points)/100+1],Y_smooth[(i*n_points/100):((i+1)*n_points)/100+1],linewidth=linewidth,color=color,alpha=alpha)
    axes.set_xlim(X.min(),X.max())
    axes.set_xlabel(xlabel,fontproperties=font, size=10, style='italic')
    axes.set_xticklabels(axes.get_xticks(),fontproperties=font, size=12)
    if '%' in ylabel:
        axes.set_ylim(0,np.minimum(2*Y.max(),100))
    axes.set_ylabel(ylabel, fontproperties=font, size=10, style='italic')
    axes.set_yticklabels(axes.get_yticks(),fontproperties=font, size=12)

def spline(self):
        """
        I return a new pycs.gen.spl.Spline object corresponding to the source.
        So this is a bit the inverse of the constructor.
        You can then put this spline object as ML of a lightcurve, as source spline, or whatever.

        Note that my output spline has LOTs of knots... it is an interpolating spline, not a regression spline !


        ..note:: This spline is a priori for display purposes only. To do an interpolation, it might be safer (and faster) to use
            the above linear interpolation eval() function.
            But making a plot, you'll see that this spline seems well accurate.
        """

        x = self.ijds.copy()
        y = self.imags.copy()
        magerrs = np.zeros(len(x))

        out = si.splrep(x, y, w=None, xb=None, xe=None, k=3, task=0, s=0.0, t=None, full_output=1, per=0, quiet=1)
        # s = 0.0 means interpolating spline !

        tck = out[0]

        # From this we want to build a real Spline object.
        datapoints = pycs.gen.spl.DataPoints(x, y, magerrs, splitup=False, sort=False, stab=False)
        outspline = pycs.gen.spl.Spline(datapoints, t = tck[0], c = tck[1], k = tck[2], plotcolour=self.plotcolour)

        # Conditions for the knots (no, everything is ok with them, no need to tweak anything).
        #intt = outspline.getintt()
        #outspline.setintt(intt)
        #print tck[2]
        #sys.exit()

        outspline.knottype = "MadeBySource"
        outspline.showknots = False # Usually we have a lot of them, slows down.
        return outspline

def getintt(self):
        """
        Returns the internal knots (i.e., not even the datapoints extrema)
        This is what you need to feed into splrep !
        There are nint - 1 such knots
        """
        return self.t[(self.k+1):-(self.k+1)].copy() # We cut the outer knots.

def r2(self, nostab=True, nosquare=False):
        """
        Evaluates the spline, compares it with the data points and returns a weighted sum of residuals r2.

        If nostab = False, stab points are included 
        This is precisely the same r2 as is used by splrep for the fit, and thus the same value as 
        returned by optc !

        This method can set lastr2nostab, so be sure to end any optimization with it.

        If nostab = True, we don't count the stab points
        """

        if nostab == True :
            splinemags = self.eval(nostab = True, jds = None)
            errs = self.datapoints.mags[self.datapoints.mask] - splinemags
            werrs = errs/self.datapoints.magerrs[self.datapoints.mask]
            if nosquare:
                r2 = np.sum(np.fabs(werrs))
            else:
                r2 = np.sum(werrs * werrs)
            self.lastr2nostab = r2
        else :
            splinemags = self.eval(nostab = False, jds = None)
            errs = self.datapoints.mags - splinemags
            werrs = errs/self.datapoints.magerrs
            if nosquare:
                r2 = np.sum(np.fabs(werrs))
            else:
                r2 = np.sum(werrs * werrs)
            self.lastr2stab = r2

        return r2

        #if red:
        #   return chi2/len(self.datapoints.jds)

def _interpolate(self, xa, xb, data):
        from scipy import interpolate   

        f = interpolate.splrep(xa, data)
        return interpolate.splev(xb, f)

def turning_points(mu, vv, Gv, Bv, dv=0.1):
    DD = np.sqrt(const.h/(8*mu*const.m_u*const.c*100))*1.0e10/np.pi
    # Gv spline
    gsp = splrep(vv, Gv, s=0)
    # Bv spline
    bsp = splrep(vv, Bv, s=0)

    # vibrational QN at which to evaluate turning points
    V = np.arange(dv-1/2, vv[-1], dv) 
    # add a point close to v=-0.5, the bottom of the well
    V = np.append([-1/2+0.0001], V)
    Rmin = []; Rmax = []; E = []
    # compute turning points using RKR method
    print(u"RKR: v   Rmin(A)  Rmax(A)  E(cm-1)")
    for vib in V:
        E.append(G(vib, gsp))    # energy of vibrational level
        ff = fg_integral(vib, gsp, bsp, lambda x, y: 1)
        gg = fg_integral(vib, gsp, bsp, B)
        fg = np.sqrt(ff/gg + ff**2)
        Rmin.append((fg - ff)*DD) # turning points
        Rmax.append((fg + ff)*DD)
        frac, integ = np.modf(vib)
        if frac > 0 and frac < dv: 
            print(u"     {:d}   {:6.4f}    {:6.4f}    {:6.2f}".
                  format(int(vib), Rmin[-1], Rmax[-1], np.float(E[-1])))

    return Rmin, Rmax, E

def formPEC(R, Rmin, Rmax, E, De, limb):
    evcm = const.e/(const.c*const.h*100) # converts cm-1 to eV

    # combine Rmin with Rmax to form PEC
    Re = (Rmin[0] + Rmax[0])/2
    print(u"RKR: Re = {:g}".format(Re))

    RTP = np.append(Rmin[::-1], Rmax, 0)  # radial positions of turning-points
    PTP = np.append(E[::-1], E, 0)  # potential energy at turning point

    # interpolate
    psp = splrep(RTP, PTP, s=0)
    # Interpolate RKR curve to this grid
    PEC = splev(R, psp, der=0)

    # extrapolate using analytical function
    inner_limb_Morse(R, PEC, RTP, PTP, Re, De)
    if limb=='L': 
        outer_limb_LeRoy(R, PEC, RTP, PTP, De)
    else:
        outer_limb_Morse(R, PEC, RTP, PTP, De)

    PTP /= evcm
    PEC /= evcm  # convert to eV

    return PEC, RTP, PTP

# analytical functions

def p2eRatio(self):
        """
        Calculates the proton density to electron density ratio using Eq. 7 of [1]_.

        Notes
        ------
        Uses the abundance and ionization equilibrium.

        References
        ----------
        .. [1] Young, P. R. et al., 2003, ApJS, `144, 135 <http://adsabs.harvard.edu/abs/2003ApJS..144..135Y>`_
        """
        if hasattr(self, 'Temperature'):
            temperature = self.Temperature
        else:
            temperature = self.IoneqAll['ioneqTemperature']
        if not hasattr(self, 'AbundanceName'):
            AbundanceName = self.Defaults['abundfile']
        else:
            AbundanceName = self.AbundanceName

        tmp_abundance = io.abundanceRead(abundancename=AbundanceName)
        abundance = tmp_abundance['abundance'][tmp_abundance['abundance']>0]
        denominator = np.zeros(len(self.IoneqAll['ioneqTemperature']))
        for i in range(len(abundance)):
            for z in range(1,i+2):
                denominator += z*self.IoneqAll['ioneqAll'][i,z,:]*abundance[i]

        p2eratio = abundance[0]*self.IoneqAll['ioneqAll'][0,1,:]/denominator
        nots = interpolate.splrep(np.log10(self.IoneqAll['ioneqTemperature']),p2eratio,s=0)
        self.ProtonDensityRatio = interpolate.splev(np.log10(temperature),nots,der=0,ext=1)

def ioneqOne(self):
        '''
        Provide the ionization equilibrium for the selected ion as a function of temperature.
        returned in self.IoneqOne
        '''
        #
        if hasattr(self, 'Temperature'):
            temperature = self.Temperature
        else:
            return
        #
        if hasattr(self, 'IoneqAll'):
            ioneqAll = self.IoneqAll
        else:
            ioneqAll = io.ioneqRead(ioneqname = self.Defaults['ioneqfile'])
            self.ioneqAll = self.IoneqAll
        #
        ioneqTemperature = ioneqAll['ioneqTemperature']
        Z = self.Z
        Ion = self.Ion
        Dielectronic = self.Dielectronic
        ioneqOne = np.zeros_like(temperature)
        #
        thisIoneq = ioneqAll['ioneqAll'][Z-1,Ion-1 + Dielectronic].squeeze()
        gioneq = thisIoneq > 0.
        goodt1 = self.Temperature >= ioneqTemperature[gioneq].min()
        goodt2 = self.Temperature <= ioneqTemperature[gioneq].max()
        goodt = np.logical_and(goodt1,goodt2)
        y2 = interpolate.splrep(np.log(ioneqTemperature[gioneq]),np.log(thisIoneq[gioneq]),s=0)
        #
        if goodt.sum() > 0:
            if self.Temperature.size > 1:
                gIoneq = interpolate.splev(np.log(self.Temperature[goodt]),y2)   #,der=0)
                ioneqOne[goodt] = np.exp(gIoneq)
            else:
                gIoneq = interpolate.splev(np.log(self.Temperature),y2)
                ioneqOne = np.exp(gIoneq)*np.ones(self.NTempDen, 'float64')
            self.IoneqOne = ioneqOne
        else:
            self.IoneqOne = np.zeros_like(self.Temperature)

def calculate_free_free_loss(self, **kwargs):
        """
        Calculate the free-free energy loss rate of an ion. The result is returned to the
        `free_free_loss` attribute.

        The free-free radiative loss rate is given by Eq. 5.15a of [1]_. Writing the numerical
        constant in terms of the fine structure constant :math:`\\alpha`,

        .. math::
           \\frac{dW}{dtdV} = \\frac{4\\alpha^3h^2}{3\pi^2m_e}\left(\\frac{2\pi k_B}{3m_e}\\right)^{1/2}Z^2T^{1/2}\\bar{g}_B

        where where :math:`Z` is the nuclear charge, :math:`T` is the electron temperature, and
        :math:`\\bar{g}_{B}` is the wavelength-averaged and velocity-averaged Gaunt factor. The
        Gaunt factor is calculated using the methods of [2]_. Note that this expression for the
        loss rate is just the integral over wavelength of Eq. 5.14a of [1]_, the free-free emission, and
        is expressed in units of erg :math:`\mathrm{cm}^3\,\mathrm{s}^{-1}`.

        References
        ----------
        .. [1] Rybicki and Lightman, 1979, Radiative Processes in Astrophysics,
            `(Wiley-VCH) <http://adsabs.harvard.edu/abs/1986rpa..book.....R>`_
        .. [2] Karzas and Latter, 1961, ApJSS, `6, 167
            <http://adsabs.harvard.edu/abs/1961ApJS....6..167K>`_
        """
        # interpolate wavelength-averaged K&L gaunt factors
        gf_kl_info = ch_io.gffintRead()
        gamma_squared = self.ionization_potential/ch_const.boltzmann/self.Temperature
        for i, atemp in enumerate(self.Temperature):
            print('%s T:  %10.2e gamma_squared  %10.2e'%(self.ion_string, atemp, gamma_squared[i]))
        gaunt_factor = splev(np.log(gamma_squared),
                             splrep(gf_kl_info['g2'],gf_kl_info['gffint']), ext=3)
        # calculate numerical constant
        prefactor = (4.*(ch_const.fine**3)*(ch_const.planck**2)/3./(np.pi**2)/ch_const.emass
                     * np.sqrt(2.*np.pi*ch_const.boltzmann/3./ch_const.emass))
        # include abundance and ionization equilibrium
        prefactor *= self.abundance*self.ioneq_one(self.stage, **kwargs)

        self.free_free_loss = prefactor*(self.Z**2)*np.sqrt(self.Temperature)*gaunt_factor

def freeFreeLoss(self, **kwargs):
        """
        Calculate the free-free energy loss rate of an ion. The result is returned to the
        `free_free_loss` attribute.

        The free-free radiative loss rate is given by Eq. 5.15a of [1]_. Writing the numerical
        constant in terms of the fine structure constant :math:`\\alpha`,

        .. math::
           \\frac{dW}{dtdV} = \\frac{4\\alpha^3h^2}{3\pi^2m_e}\left(\\frac{2\pi k_B}{3m_e}\\right)^{1/2}Z^2T^{1/2}\\bar{g}_B

        where where :math:`Z` is the nuclear charge, :math:`T` is the electron temperature, and
        :math:`\\bar{g}_{B}` is the wavelength-averaged and velocity-averaged Gaunt factor. The
        Gaunt factor is calculated using the methods of [2]_. Note that this expression for the
        loss rate is just the integral over wavelength of Eq. 5.14a of [1]_, the free-free emission, and
        is expressed in units of erg :math:`\mathrm{cm}^3\,\mathrm{s}^{-1}`.

        References
        ----------
        .. [1] Rybicki and Lightman, 1979, Radiative Processes in Astrophysics,
            `(Wiley-VCH) <http://adsabs.harvard.edu/abs/1986rpa..book.....R>`_
        .. [2] Karzas and Latter, 1961, ApJSS, `6, 167
            <http://adsabs.harvard.edu/abs/1961ApJS....6..167K>`_
        """
        # interpolate wavelength-averaged K&L gaunt factors
        gf_kl_info = ch_io.gffintRead()
        gamma_squared = self.IprErg/ch_const.boltzmann/self.Temperature
#        for i, atemp in enumerate(self.Temperature):
#            print('%s T:  %10.2e gamma_squared  %10.2e'%(self.ion_string, atemp, gamma_squared[i]))
        gaunt_factor = splev(np.log(gamma_squared),
                             splrep(gf_kl_info['g2'],gf_kl_info['gffint']), ext=3)
        # calculate numerical constant
        prefactor = (4.*(ch_const.fine**3)*(ch_const.planck**2)/3./(np.pi**2)/ch_const.emass
                     * np.sqrt(2.*np.pi*ch_const.boltzmann/3./ch_const.emass))
        # include abundance and ionization equilibrium
        prefactor *= self.abundance*self.IoneqOne

        self.FreeFreeLoss = {'rate':prefactor*(self.Z**2)*np.sqrt(self.Temperature)*gaunt_factor}

def ioneqOne(self):
        '''
        Provide the ionization equilibrium for the selected ion as a function of temperature.
        Similar to but not identical to ion.ioneqOne()
        returned in self.IoneqOne
        '''
        #
        if hasattr(self, 'Temperature'):
            temperature = self.Temperature
        else:
            return
        #
        if hasattr(self, 'IoneqAll'):
            ioneqAll = self.IoneqAll
        else:
            self.IoneqAll = ch_io.ioneqRead(ioneqname = self.Defaults['ioneqfile'])
            ioneqAll = self.IoneqAll
        #
        ioneqTemperature = ioneqAll['ioneqTemperature']
        Z = self.Z
        stage = self.stage
        ioneqOne = np.zeros_like(temperature)
        #
        thisIoneq = ioneqAll['ioneqAll'][Z-1,stage-1].squeeze()
        gioneq = thisIoneq > 0.
        goodt1 = self.Temperature >= ioneqTemperature[gioneq].min()
        goodt2 = self.Temperature <= ioneqTemperature[gioneq].max()
        goodt = np.logical_and(goodt1,goodt2)
        y2 = splrep(np.log(ioneqTemperature[gioneq]),np.log(thisIoneq[gioneq]),s=0)
        #
        if goodt.sum() > 0:
            if self.Temperature.size > 1:
                gIoneq = splev(np.log(self.Temperature[goodt]),y2)   #,der=0)
                ioneqOne[goodt] = np.exp(gIoneq)
            else:
                gIoneq = splev(np.log(self.Temperature),y2)
                ioneqOne = np.exp(gIoneq)
                ioneqOne = np.atleast_1d(ioneqOne)
            self.IoneqOne = ioneqOne
        else:
            self.IoneqOne = np.zeros_like(self.Temperature)

def __init__(self, t, x, y, omega, origin):
        self.t      = t
        self.x      = x
        self.y      = y
        self.omega  = omega
        self.origin = origin

        self.n = len(t)

        self.x_spl     = interpolate.splrep(self.t, self.x)
        self.y_spl     = interpolate.splrep(self.t, self.y)
        self.omega_spl = interpolate.splrep(self.t, self.omega)

def draw_spline(image, lane, ys, color=[0, 0, 255]):
    pts = [[x, ys[i]]for i, x in enumerate(lane) if not math.isnan(x)]
    if len(pts)-1 > 0:
        spline = interpolate.splrep([pt[1] for pt in pts], [pt[0] for pt in pts], k=len(pts)-1)
        for i in range(min([pt[1] for pt in pts]), max([pt[1] for pt in pts])):
            image[i, int(interpolate.splev(i, spline)), :] = color
    return image

def __init__(self, A, h, alpha, CL, CD, rho, smoothing=0, k_spline=3):
        self.A   = A
        self.h   = h
        self.Asp = 2*self.h**2/self.A

        self.rho = rho

        # Input lift and drag data
        self.n     = len(alpha)
        self.alpha = alpha
        self.CL    = CL
        self.CD    = CD

        self.k_spline = k_spline

        # -------- Spline interpolation -----------------------------
        if len(self.alpha.shape) == 1:
            self.interpolationMethod = 'spline'
            self.nrControlParameters = 1
            self.CL_spline = interpolate.splrep(self.alpha, self.CL, s=smoothing, k=self.k_spline)
            self.CD_spline = interpolate.splrep(self.alpha, self.CD, s=smoothing, k=self.k_spline)
        elif len(self.alpha.shape) == 2:
            self.interpolationMethod = 'griddata'
            self.nrControlParameters = 2
            self.CL_RbfModel = interpolate.Rbf(self.alpha[:, 0],self.alpha[:, 1], self.CL, smooth=smoothing)
            self.CD_RbfModel = interpolate.Rbf(self.alpha[:, 0],self.alpha[:, 1], self.CD, smooth=smoothing)

def setCalmWaterResistanceData(self, U, Fr, Re, CR, Cf):
        self.Fr = Fr
        self.U  = U
        self.Re = Re
        self.CR = CR
        self.Cf = Cf

        self.CRspl = interpolate.splrep(self.Fr, self.CR, s=1e-9)
        self.Cfspl = interpolate.splrep(self.Fr, self.Cf, s=1e-9)

def necessaryDriftAngle(self, U, Fy, scale=1):
        alpha_test = np.linspace(0, np.max(self.alpha), 5)

        Fx_test, Fy_test, Mz_test = self.driftInducedForces(U, alpha_test, scale=scale)

        alpha_spline = interpolate.splrep(Fy_test, alpha_test)

        return interpolate.splev(Fy, alpha_spline)

def reqRevolutions(self, T_req, U):
        J_min = 0.05

        n_min = U/(self.J_max*self.D)
        n_max = U/(J_min*self.D)

        n_test = 10
        n = np.linspace(n_min, n_max, n_test)

        T = self.Thrust(U, n)

        n_spl = interpolate.splrep(T, n)

        return float(interpolate.splev(T_req, n_spl))

def figure(self,mode,data,name,**options):
        '''
        Generate a figure to view the data.

        Parameters
        ----------
        mode : 'L','P'
            'L' for lines and 'P' for pcolor.
        data : ndarray
            The data to be viewed.
        name : str
            The name of the figure.
        options : dict
            Other options.
        '''
        assert mode in ('L','P')
        plt.title(os.path.basename(name))
        if mode=='L':
            if options.get('interpolate',False):
                plt.plot(data[:,0],data[:,1],'r.')
                X=np.linspace(data[:,0].min(),data[:,0].max(),10*data.shape[0])
                for i in xrange(1,data.shape[1]):
                    tck=interpolate.splrep(data[:,0],data[:,i],k=3)
                    Y=interpolate.splev(X,tck,der=0)
                    plt.plot(X,Y,label=options['legend'][i-1] if 'legend' in options else None)
                if 'legend' in options:
                    leg=plt.legend(fancybox=True,loc=options.get('legendloc',None))
                    leg.get_frame().set_alpha(0.5)
            else:
                plt.plot(data[:,0],data[:,1:])
                if 'legend' in options:
                    leg=plt.legend(options['legend'],fancybox=True,loc=options.get('legendloc',None))
                    leg.get_frame().set_alpha(0.5)
        elif mode=='P':
            plt.colorbar(plt.pcolormesh(data[:,:,0],data[:,:,1],data[:,:,2]))
            if 'axis' in options: plt.axis(options.get('axis','equal'))
        if self.show and self.suspend: plt.show()
        if self.show and not self.suspend: plt.pause(App.SUSPEND_TIME)
        if self.savefig: plt.savefig('%s.png'%name)
        plt.close()

def interpolate_to_mesh_indices(self, loop):
        """
        Return interpolated loop indices to the temperature and density meshes defined for
        the atomic data. For use with `~scipy.ndimage.map_coordinates`.
        """
        nots_itemperature = splrep(self.temperature.value, np.arange(self.temperature.shape[0]))
        nots_idensity = splrep(self.density.value, np.arange(self.density.shape[0]))
        itemperature = splev(np.ravel(loop.electron_temperature.value), nots_itemperature)
        idensity = splev(np.ravel(loop.density.value), nots_idensity)

        return itemperature, idensity

def process(self, spectrogram):
        mult = np.multiply(self.weights, np.ones((1, spectrogram.shape[1])))
        square_mag = np.power(spectrogram.magnitude(),2)
        square_mag = np.vstack((square_mag[1:,:],np.flipud(square_mag[1:,:])))
        square_mag_sum = np.sum(square_mag, 0)
        energy_threshold = np.percentile(square_mag_sum, 25)
        res = np.fft.rfft(square_mag.T, 2 * (self.n_bins - 1)).T
        resN = np.abs(res)
        resP = np.angle(res)
        yin = np.zeros((spectrogram.shape))
        yin[0,:] = 1
        tmp = np.zeros((spectrogram.shape[1],))
        for tau in range(1, self.n_bins):
            yin[tau,:] = square_mag_sum - resN[tau,:] * np.cos(resP[tau,:])
            tmp += yin[tau,:]
            yin[tau,:] *= tau/tmp

        tau = self.tau_min + np.argmin(yin[self.tau_min:self.tau_max,:],0)
        y_min = np.min(yin[self.tau_min:self.tau_max,:],0)
        tau = tau[:,np.newaxis]
        if self.interp:
            ranges = np.hstack((tau - 5, tau - 4, tau - 3, tau-2,tau-1,tau, 
                tau+1,tau + 2, tau + 3, tau + 4, tau + 5))
            new_tau = np.empty_like(tau)
            for frame in range(spectrogram.shape[1]):                    
                r = ranges[frame]
                r[0] = max(r[0], 0)
                r[1] = min(r[1], self.n_bins-1)
                val = yin[r.astype(int), frame]
                tck = interpolate.splrep(r, val, s=0, k = 2)
                xnew = np.arange(r[0],r[-1], (r[-1] - r[0]) /10)
                ynew = interpolate.splev(xnew, tck, der=0)
                new_tau[frame] = xnew[np.argmin(ynew)]
                y_min[frame] = np.min(ynew)            
            tau = new_tau
        tau = tau[:,0]
        pitch_confidence = 1- y_min
        pitch = np.zeros((spectrogram.shape[1],))
        pitch[tau!=0] = np.nan_to_num(self.sample_rate * 1.0 / (tau[tau!=0]))
        pitch[tau==0] = 0
        pitch_confidence[tau==0] = 0

        pitch[square_mag_sum < energy_threshold] = 0
        pitch_confidence[square_mag_sum < energy_threshold] = 0
        return (pitch, pitch_confidence)