def __startBarionicAccretionRate(self):

        np = 1000
        deltaz = self._zmax / float(np)

        z = [self._zmax - i * deltaz for i in range(np)]
        z.append(0)
        z = array(z)
        fbt2 = array([self.fbstruc(zi) for zi in z])
        ascale = array([1.0 / (1.0 + zi) for zi in z])
        self._ascale = ascale

        tck = spint.splrep(ascale, fbt2)
        ab3 = spint.splev(ascale, tck, der=1)

        def a5(z, i):
            a = 1.0 / (1.0 + z)
            a2 = a * a
            a3 = -1.0 * ab3[i] * a2
            a4 = a3
            a5 = self._cosmology.getRobr0() * abs(a4) \
                 / self._cosmology.dt_dz(z)
            return a5

        self._abt2 = array([a5(z[i], i) for i in range(z.size)])
        self._tck_ab = spint.splrep(self._ascale, self._abt2)

    def RealParts(self, om):
        """ calculate functions on real axis
        """
        Func=[]
        derivs=[]
        for i in range(len(self.pos)):
            En = self.pos[i]

            (x0, dh0) = swing_make_mesh(300, 1e-2*abs(En), 2*abs(En), 7*abs(En))

            wb = array([self.F0(x+En, En) for x in x0])
            
            F0j = array([self.F0(x, En) for x in om])

            Frc = kramskron(om, F0j, wb, x0, dh0, En)

            spl = interpolate.splrep(om, real(Frc), k=3, s=0.0)            
            dersr = [interpolate.splev(0.0, spl, der=0),
                     interpolate.splev(0.0, spl, der=1),
                     interpolate.splev(0.0, spl, der=2)]

            spl = interpolate.splrep(om, imag(Frc), k=3, s=0.0)
            dersi = [interpolate.splev(0.0, spl, der=0),
                     interpolate.splev(0.0, spl, der=1),
                     interpolate.splev(0.0, spl, der=2)]

            Func.append(Frc)
            ders = array(dersr + dersi)
            derivs.append(array(ders))
        return (Func, derivs)

	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 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 setGrid(self):
     from scipy import interpolate
     x0 = numpy.logspace(-3,1,81)
     etas = numpy.linspace(0.,2.,21)
     qs = numpy.linspace(0.2,1.,17)
     grid1 = numpy.empty((x0.size,x0.size,etas.size,qs.size))
     grid2 = numpy.empty(grid1.shape)
     for i in range(qs.size):
         q = qs[i]
         q2 = q**2
         b = 1-q2
         for j in range(etas.size):
             eta = etas[j]
             g = 0.5*eta-1. # g = -1*gamma
             for k in range(x0.size):
                 x = x0[k]
                 for l in range(x0.size):
                     y = x0[l]
                     qb = ((2*x*y)/b)**g  # q_bar
                     qt = q*(x*y)**0.5/b  # q_tilde
                     sb = 0.5*(x/y - y/x) + s**2*b/(2*x*y)
                     nu1 = s**2*b/(2*x*y)
                     nu2 = nu1+ 0.5*b*(x/y + y/(x*q2))
                     nu = numpy.logspace(nu1,nu2,1001)
                     mu = nu-sb
                     t = (1+mu**2)**0.5
                     f1 = (t-mu)**0.5/t
                     f2 = (t+mu)**0.5/t
                     ng = nu**g
                     I1 = interpolate.splrep(nu,f1*ng)
                     I2 = interpolate.splrep(nu,f2*ng)
                     grid1[k,l,i,j] = qt*interpolate.splint(nu1,nu2,I1)
                     grid2[k,l,i,j] = qt*interpolate.splint(nu1,nu2,I2)
             pylab.imshow(grid1[:,:,i,j])
             pylab.show()

def VegaFilterMagnitude(filter,spectrum,redshift):
    """
    Determines the Vega magnitude (up to a constant) given an input filter,
        SED, and redshift.
    """
    from scipy.interpolate import splev,splint,splrep
    from scipy.integrate import simps
    from math import log10

    wave = spectrum[0].copy()
    data = spectrum[1].copy()

    # Redshift the spectrum and determine the valid range of wavelengths
    wave *= (1.+redshift)
    data /= (1.+redshift)
    wmin,wmax = filter[0][0],filter[0][-1]
    cond = (wave>=wmin)&(wave<=wmax)

    # Evaluate the filter at the wavelengths of the spectrum
    response = splev(wave[cond],filter)

    # Determine the total observed flux (without the bandpass correction)
    observed = splrep(wave[cond],(response*data[cond]),s=0,k=1)
    flux = splint(wmin,wmax,observed)

    # Determine the magnitude of Vega through the filter
    vwave,vdata = getSED('Vega')
    cond = (vwave>=wmin)&(vwave<=wmax)
    response = splev(vwave[cond],filter)
    vega = splrep(vwave[cond],response*vdata[cond],s=0,k=1)
    vegacorr = splint(wmin,wmax,vega)

    return -2.5*log10(flux/vegacorr)#+2.5*log10(1.+redshift)

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 GetSourceSize(self,kpc=False):
     self.z=source_redshifts[self.name]
     self.Da = astCalc.da(self.z)
     self.scale = self.Da*1e3*np.pi/180./3600.
     if len(self.srcs) == 1 or self.name == 'J0837':
         self.Re_v = self.Ddic['Source 1 re']*0.05
         self.Re_i = self.Re_v.copy()
         self.Re_lower = self.Ldic['Source 1 re']*0.05
         self.Re_upper = self.Udic['Source 1 re']*0.05
     elif len(self.srcs) == 2 and self.name != 'J0837':
         print 'test this out...!'
         Xgrid = np.logspace(-4,5,1501)
         Res = []
         for i in range(len(self.imgs)):
             #if self.name == 'J1605':
             #    source = 
             source = self.fits[i][-3]*self.srcs[0].eval(Xgrid) + self.fits[i][-2]*self.srcs[1].eval(Xgrid)
             R = Xgrid.copy()
             light = source*2.*np.pi*R
             mod = splrep(R,light,t=np.logspace(-3.8,4.8,1301))
             intlight = np.zeros(len(R))
             for i in range(len(R)):
                 intlight[i] = splint(0,R[i],mod)
             model = splrep(intlight[:-300],R[:-300])
             
             if len(model[1][np.where(np.isnan(model[1])==True)]>0):
                 print "arrays need to be increasing monotonically! But don't worry about it"
                 model = splrep(intlight[:-450],R[:-450])
             reff = splev(0.5*intlight[-1],model)
             Res.append(reff*0.05)
         self.Re_v,self.Re_i = Res
     if kpc:
         return [self.Re_v*self.scale, self.Re_i*self.scale]
     return [self.Re_v, self.Re_i]

def Interpo (wave, flux, variance) :
    wave_min = 1500
    wave_max = 12000
    dw = 2

    #wavelength = np.linspace(wave_min,wave_max,(wave_max-wave_min)/pix+1)
    wavelength = np.arange(math.ceil(wave_min), math.floor(wave_max), dtype=int, step=dw) #creates N equally spaced wavelength values
    inter_flux = []
    inter_var  = []
    output     = []

    lower = wave[0] # Find the area where interpolation is valid
    upper = wave[-1]

    good_data = np.where((wave >= lower) & (wave <= upper))	#creates an array of wavelength values between minimum and maximum wavelengths from new spectrum

    influx = inter.splrep(wave[good_data], flux[good_data])	#creates b-spline from new spectrum

    invar  = inter.splrep(wave[good_data], variance[good_data]) # doing the same with the errors

    inter_flux = inter.splev(wavelength, influx)	#fits b-spline over wavelength range
    inter_var  = inter.splev(wavelength, invar)   # doing the same with errors

    missing_data = np.where((wavelength < lower) | (wavelength > upper))
    inter_flux[missing_data] = float('NaN')  # set the bad values to NaN !!!
    inter_var[missing_data] =  float('NaN')

    output = np.array([wavelength, inter_flux, inter_var]) # put the interpolated data into the new table

    return output # return new table

def qvRealtion(voltage):

    global chargeVoltage, chargeQuantity
	
    vPoints = chargeVoltage
    qPoints = chargeQuantity
    SOCPoints = chargeSOC
	
    qvCurveTck = interpolate.splrep(vPoints, qPoints, s = 0.002)
    batteryQuantity = interpolate.splev(voltage, qvCurveTck, der = 0)
    SOCVCurveTck = interpolate.splrep(SOCPoints, vPoints,s = 0.002)
    batterySOC = interpolate.splev(voltage, SOCVCurveTck, der = 0)
	
    newX = np.arange(vPoints[0], vPoints[-1], 0.1)
    newY = interpolate.splev(newX, qvCurveTck, der = 0)
	
    plt.plot(vPoints, qPoints, 'o', newX, newY, '-')
    plt.title("V-Q curve")
	
    newSOC = np.arange(SOCPoints[0], SOCPoints[-1], 0.01)
    newVoltage = interpolate.splev(newSOC, SOCVCurveTck, der = 0)
	
    plt.plot(SOCPoints, vPoints, 'o', newSOC, newVoltage, '-')
    plt.title("V-SOC curve")
    plt.show()
	
    return batteryQuantity

def conv2(lam,flu,ac,bc,rot):
	
	lux = 299792.458
	
	WAV = np.array([3530,4860,6260,7670,8350])
		
	"""Determining Limb Darkening Coefficients"""
	tck = interpolate.splrep(WAV,ac,k=3,s=0)
	A = interpolate.splev(lam,tck,der=0)
	tck = interpolate.splrep(WAV,bc,k=3,s=0)
	B = interpolate.splev(lam,tck,der=0)
		
	F = flu.copy()
	
	MIN = lam * ( 1. - float(rot) / lux )
	MAX = lam * ( 1. + float(rot) / lux )
	
	lar = len(F)
	
	if rot > 0.0:
		H = Cfunctions.Conv(lam.astype("double"),flu.astype("double"),F.astype("double"), \
		    A.astype("double"),B.astype("double"),MIN.astype("double"),MAX.astype("double"),rot,lar)
	else:
		H = flu
	
	Fnu = np.array(H)

	return Fnu

def envelope(min_arr, max_arr, N, periodic=0):
    #Cubic Spline by default
    order_max = 3
    order_min = 3

    min_arr = np.asarray(min_arr)
    max_arr = np.asarray(max_arr)

    if min_arr.shape[1] < 4:
        order_min = 1 #Do linear interpolation if not enough points
    elif min_arr.shape[1] < 5:
        order_min = 2 #Do quad interpolation if not enough points
    else:
        order_min = 3

    if max_arr.shape[1] < 4:
        order_max = 1  #Do linear interpolation if not enough points
    elif max_arr.shape[1] < 5:
        order_max = 2 #Do quad interpolation if not enough points
    else:
        order_max = 3
    # Mirror Method requires per flag = 1
    # No extrapolation requires per flag = 0
    t = interpolate.splrep(*min_arr, k=order_min, per=periodic)
    top = interpolate.splev(np.arange(N), t)

    b = interpolate.splrep(*max_arr, k=order_max, per=periodic)
    bot = interpolate.splev(np.arange(N), b)
    mean = (top + bot)/2

    return mean

def get_optimal_splines(events, optimise_bin_edges, k=1):

    cut_events = {}
    cut_energies, ga_cuts, xi_cuts = [], [], []
    for elow, ehigh, emid in zip(optimise_bin_edges[:-1],
                                 optimise_bin_edges[1:],
                                 np.sqrt(optimise_bin_edges[:-1] *
                                         optimise_bin_edges[1:])):

        for key in events:
            cut_events[key] = events[key][
                (events[key]["MC_Energy"] > elow) &
                (events[key]["MC_Energy"] < ehigh)]

        res = optimize.differential_evolution(
            cut_and_sensitivity,
            bounds=[(.5, 1), (0, 0.5)],
            maxiter=1000, popsize=10,
            args=(cut_events,
                  np.array([elow / energy_unit,
                            ehigh / energy_unit]) * energy_unit,
                  alpha)
        )

        if res.success:
            cut_energies.append(emid.value)
            ga_cuts.append(res.x[0])
            xi_cuts.append(res.x[1])

    spline_ga = interpolate.splrep(cut_energies, ga_cuts, k=k)
    spline_xi = interpolate.splrep(cut_energies, xi_cuts, k=k)

    return (spline_ga, ga_cuts), (spline_xi, xi_cuts)

def fit_fwhm_enclosed_direct(peak, rad, flux):

    # We use splines to interpolate and derivate
    spl = splrep(rad, flux)
    # First derivative
    vald1 = splev(rad, spl, der=1)

    splinter = splrep(rad, vald1 - math.pi * peak * rad)

    roots = sproot(splinter)
    nroots = len(roots)
    if peak < 0:
        msg = "The method doesn't converge, peak is negative"
        fwhm = -99
    else:
        if nroots == 0:
            msg = "The method doesn't converge, no roots"
            fwhm = -99
        elif nroots == 1:
            r12 = roots[0]
            fwhm = 2 * r12
            msg = "The method converges, one root"
        else:
            msg = "The method doesn't converge, multiple roots"
            r12 = roots[0]
            fwhm = 2 * r12

    return peak, fwhm, msg

    def calc_spline_conv_coefs(self):
        """
        Calculate cubic spline coefficients for a set of curves.

        INPUT:
          None

        RETURN:
          None

        For evaluation/interpolation: 
         Make sure that splev is imported from scipy.interpolate.
         Call like:      ynew = splev(xnew,tck,der=0)
         or for derivs   yder = splev(xnew,tck,der=1).
        """

        # do some checks
        if(self.y.shape[0] != self.x.shape[0]):
            raise Exception, "".join(["dimensions do not agree for: y[n,m]",
                                      "and x[n]"])

        self.tck=[]
        if(len(self.y.shape) == 1):
            self.tck.append(interpolate.splrep(self.x, self.y, s=0))
        else:
            n = self.y.shape[1]
            for i in range(n):
                self.tck.append(interpolate.splrep(self.x, self.y[:,i], s=0))

        self.__kind_computed__ |= self.CONV_SPLINE
        self.importQuantity('tck', self.tck)
        self.importQuantity('__kind_computed__', self.__kind_computed__)

        return

	def pgrid(self,tname):

		self._hasTeff_ = True
		self._hasLogg_ = True
		self._hasMet_ = False
		
		grid = Table.read(tname,format='ascii.no_header')

		len_teff = len(np.unique(grid['col2']))
		len_logg = len(np.unique(grid['col3']))

		# Reads Teff

		self.gteff = np.zeros(len_teff)

		for i,ii in enumerate(np.arange(len_teff)*len_logg):
			self.gteff[i] = np.float(os.path.basename(grid['col1'][ii])[2:].split('_')[0])

		# Reads logg
		
		self.glogg = np.zeros(len_logg)

		for i,ii in enumerate(np.arange(len_logg)):
			self.glogg[i] = np.float(os.path.basename(grid['col1'][ii]).split('_')[1].split('.')[0])

		# evalute grid interpolation

		self._teffc = itrp.splrep(np.arange(len_teff),self.gteff)
		self._loggc = itrp.splrep(np.arange(len_logg),self.glogg)

    def calc_spline_sqrt_coefs(self):
        """
        Calculate cubic spline coefficients for a set of curves.

        INPUT:
          None

        RETURN:
          None

        For further usage compute spline coefficients for
        rmnc_b(s), zmns_b(s), pmns_b(s).

        Divide odd m terms by sqrt(s), 
        keep even m terms as they are.
        Use 3rd order qubic spline.

        ATTENTION: For interpolation keep in mind to 're-do' sqrt(s)!

        For evaluation/interpolation: 
         Use splev_sqrt(). This routine calls splev() which is
         to be found in from scipy.interpolate.
         Make sure that splev is imported from scipy.interpolate.
         Call like:      ynew = splev_sqrt(xnew,sqrt_tck,der=0)
         or for derivs   yder = splev_sqrt(xnew,sqrt_tck,der=1).
        """

        # as we shall not divide by zero...
        eps = 1.e-30
        self.x = numpy.where(self.x >= eps, self.x, eps)

        # do some checks
        if(self.y.shape[0] != self.x.shape[0]):
            raise Exception, "".join(["dimensions do not agree for: y[n,m]",
                                      "and x[n]"])
        if(self.y.shape[1] != self.ixm.shape[0]):
            raise Exception, "".join(["dimensions do not agree for: y[n,m]",
                                      "and ixm[m]"])

        self.sqrt_tck=[]
        if(len(self.y.shape) == 1):
            self.sqrt_tck.append(interpolate.splrep(self.x, self.y, s=0))
        else:
            sqx = numpy.sqrt(self.x)
            n = self.y.shape[1]
            for i in range(n):
                if(self.ixm[i]%2 == 0):
                    self.sqrt_tck.append(interpolate.splrep(self.x, 
                                                            self.y[:,i], 
                                                            s=0))
                else:
                    self.sqrt_tck.append(interpolate.splrep(self.x, 
                                                            self.y[:,i]/sqx, 
                                                            s=0))

        self.__kind_computed__ |= self.SQRT_SPLINE
        self.importQuantity('sqrt_tck', self.sqrt_tck)
        self.importQuantity('__kind_computed__', self.__kind_computed__)

        return

    def redshiftfunctions(self):   
        D=self.D
        zbins=self.zbins
        z2bins=self.z2bins
        Dabins=zbins*0.0
        Dmodbins=zbins*0.0
        Da2bins=numpy.zeros((z2bins.size,z2bins.size))
        volumebins=zbins*0.0
        for i in range(zbins.size):
            Dabins[i]=D.Da(zbins[i])
            Dmodbins[i]=D.distance_modulus(zbins[i])
            volumebins[i]=D.volume(zbins[i])
        for i in range(z2bins.size):
            for j in range(z2bins.size):
                if j>i:
                    Da2bins[i,j]=D.Da(z2bins[i],z2bins[j])

        self.Da_spline=interpolate.splrep(zbins,Dabins)
        self.Dmod_spline=interpolate.splrep(zbins,Dmodbins)

        self.volume_spline=interpolate.splrep(zbins,volumebins)

        z2d=iT.coords((z2bins.size,z2bins.size))*self.dz2
        self.Da_bispline=interpolate.RectBivariateSpline(z2bins,z2bins,Da2bins)

        #pickle the splines
        splinedump=open("redshiftsplines.pkl","wb")
        cPickle.dump([self.Da_spline,self.Dmod_spline,self.volume_spline,self.Da_bispline],splinedump,2)

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 bspline(pts, degree, smoothness):
    # cycle check
    periodic = False
    if pts[0] == pts[-1]:
        # pts.pop()
        periodic = True
        pts = pts[:-1]

    if periodic: pts = pts + pts[0 : degree+1]
    else: pts = [pts[0]] + pts + [pts[-1],pts[-1]]

    pts = array(pts)
    n_points = len(pts)
    x, y = pts[:,0], pts[:,1]

    t = range(len(x))
    ipl_t = linspace(1.0, len(pts) - degree, smoothness)

    x_tup = si.splrep(t, x, k=degree, per=periodic)
    y_tup = si.splrep(t, y, k=degree, per=periodic)
    x_list = list(x_tup)
    xl = x.tolist()

    y_list = list(y_tup)
    yl = y.tolist()

    if periodic:
        x_list[1] = [0.0] + xl + [0.0, 0.0, 0.0, 0.0]
        y_list[1] = [0.0] + yl + [0.0, 0.0, 0.0, 0.0]

    x_i = si.splev(ipl_t, x_list)
    y_i = si.splev(ipl_t, y_list)

    return zip(x_i, y_i)