Python integrate.simps方法代码示例

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import simps [as 别名]
def bunching(p_array, lambda_mod, smooth_sigma=None):
    """
    Function calculates bunching factor for wavelength lambda_mod

    $b(\lambda) = \frac{1}{N_0}\left| \langle e^{- i\frac{ 2 \pi}{\lambda} s} N(s)\rangle \right|$

    :param p_array: ParticleArray
    :param lambda_mod: wavelength
    :param smooth_sigma: smoothing parameter
    :return: bunching factor
    """
    if smooth_sigma is None:
        smooth_sigma = min(np.std(p_array.tau())*0.01, lambda_mod*0.1)

    B = s_to_cur(p_array.tau(), sigma=smooth_sigma, q0=np.sum(p_array.q_array), v=speed_of_light)

    b = np.abs(simps(B[:, 1] / speed_of_light * np.exp(-1j * 2 * np.pi / lambda_mod * B[:, 0]), B[:, 0])) / np.sum(
        p_array.q_array)
    return b 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import simps [as 别名]
def slice_bunching(tau, charge, lambda_mod, smooth_sigma=None):
    """
    Function calculates bunching factor for wavelength lambda_mod

    $b(\lambda) = \frac{1}{N_0}\left| \langle e^{- i\frac{ 2 \pi}{\lambda} s} N(s)\rangle \right|$

    :param p_array: ParticleArray
    :param lambda_mod: wavelength
    :param smooth_sigma: smoothing parameter
    :return: bunching factor
    """
    if smooth_sigma is None:
        smooth_sigma = lambda_mod/10.

    B = s_to_cur(tau, sigma=smooth_sigma, q0=charge, v=speed_of_light)

    b = np.abs(simps(B[:, 1] / speed_of_light * np.exp(-1j * 2 * np.pi / lambda_mod * B[:, 0]), B[:, 0])) / charge
    return b 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import simps [as 别名]
def p_quad_values(self, mode, xvec, pvec):

        r"""Calculates the discretized p-quadrature probability distribution of the specified mode.

        Args:
            mode (int): the mode to calculate the p-quadrature probability values of
            xvec (array): array of discretized :math:`x` quadrature values
            pvec (array): array of discretized :math:`p` quadrature values

        Returns:
            array: 1D array of size len(pvec), containing reduced p-quadrature
            probability values for a specified range of x and p.
        """

        W = self.wigner(mode, xvec, pvec)
        y = []
        for i in range(0, len(pvec)):
            res = simps(W[i, : len(xvec)], xvec)
            y.append(res)
        return np.array(y) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import simps [as 别名]
def x_quad_values(self, mode, xvec, pvec):

        r"""Calculates the discretized x-quadrature probability distribution of the specified mode.

        Args:
            mode (int): the mode to calculate the x-quadrature probability values of
            xvec (array): array of discretized :math:`x` quadrature values
            pvec (array): array of discretized :math:`p` quadrature values

        Returns:
            array: 1D array of size len(xvec), containing reduced x-quadrature
            probability values for a specified range of x and p.
        """

        W = self.wigner(mode, xvec, pvec)
        y = []
        for i in range(0, len(xvec)):
            res = simps(W[: len(pvec), i], pvec)
            y.append(res)
        return np.array(y) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import simps [as 别名]
def _integrate(self, y, f):
        """
        Integrate `y` along axis=1, i.e. over freq axis for all T.

        Parameters
        ----------
        y : 2d array (nT, ndos) where nT = len(self.T), ndos = len(self.dos)
        f : self.f, (len(self.dos),)

        Returns
        -------
        array (nT,)
        """
        mask = np.isnan(y)
        if mask.any():
            self._printwarn("HarmonicThermo._integrate: warning: "
                            " %i NaNs found in y!" %len(mask))
            if self.fixnan:
                self._printwarn("HarmonicThermo._integrate: warning: "
                                "fixing %i NaNs in y!" %len(mask))
                y[mask] = self.nanfill
        # this call signature works for scipy.integrate,{trapz,simps}
        return self.integrator(y, x=f, axis=1) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import simps [as 别名]
def expected(data, airFrame):
    alpha, V, lift_to_drag = data

    pdf = airFrame.pdf.score_samples(np.vstack([alpha.ravel(), V.ravel()]).T)
    pdf = np.exp(pdf.reshape(lift_to_drag.shape))
    expected_value = 0
    numerator_list = []
    denominator_list = []
    for i in range(len(lift_to_drag)):
        numerator = simps(lift_to_drag[i]*pdf[i], alpha[i])
        denominator = simps(pdf[i], alpha[i])
        numerator_list.append(numerator)
        denominator_list.append(denominator)
    numerator = simps(numerator_list, V[:, 0])
    denominator = simps(denominator_list, V[:, 0])
    expected_value = numerator/denominator
    return(expected_value) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import simps [as 别名]
def simps_approx():
    # Compute the left hand side analytically
    loglhs = psi*a - b * np.log1p(np.exp(psi))
    lhs = np.exp(loglhs)

    # Compute the right hand side with quadrature
    from scipy.integrate import simps
    # Lay down a grid of omegas
    # TODO: How should we choose the bins?
    omegas = np.linspace(1e-15, 5, 1000)
    # integrand = lambda om: 2**-b \
    #                        * np.exp((a-b/2.)*psi 0 0.5*om*psi**2) \
    #                        * polya_gamma_density(om, b, 0)
    # y = map(integrand, omegas)
    # rhs = simps(integrand, y)
    logy = -b * np.log(2) + (a-b/2.) * psi -0.5*omegas*psi**2
    logy += log_polya_gamma_density(omegas, b, 0, trunc=21)
    y = np.exp(logy)
    rhs = simps(y, omegas)

    print("Numerical Quadrature")
    print("log LHS: ", loglhs)
    print("log RHS: ", np.log(rhs)) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import simps [as 别名]
def plot_density():
    omegas = np.linspace(1e-16, 5, 1000)
    logpomega = log_polya_gamma_density(omegas, b, 0, trunc=1000)
    pomega = np.exp(logpomega).real

    import matplotlib.pyplot as plt
    plt.ion()
    plt.figure()
    plt.plot(omegas, pomega)

    y = -b * np.log(2) + (a-b/2.) * psi -0.5*omegas*psi**2


    from scipy.integrate import simps
    Z = simps(pomega, omegas)
    print("Z: ", Z) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import simps [as 别名]
def get_dband_center(self, d_cols):
        """
        Get d-band center of the DosX object.

        Parameters:
        -----------
        d_cols: The column number range for d orbitals, int or tuple of int.

        Examples:
        ---------
        # The 5 - 9 columns are state density for d orbitals.
        >>> dos.get_dband_center(d_cols=(5, 10))
        """
        d_cols = (d_cols, d_cols+1) if type(d_cols) is int else d_cols

        # 合并d轨道DOS
        start, end = d_cols
        yd = np.sum(self.data[:, start:end], axis=1)

        #获取feimi能级索引
        for idx, E in enumerate(self.data[:, 0]):
            if E >= 0:
                nfermi = idx
                break
        E = self.data[: nfermi+1, 0]  # negative inf to Fermi
        dos = yd[: nfermi+1]          # y values from negative inf to Fermi
        # Use Simpson integration to get d-electron number
        nelectro = simps(dos, E)
        # Get total energy of dband
        tot_E = simps(E*dos, E)
        dband_center = tot_E/nelectro
        self.dband_center = dband_center

        return dband_center 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import simps [as 别名]
def flux_total(self):
        C_fi = 3.9614e19 #ph/(sec * rad * GeV * A)
        mrad = 1e-3 # transform rad to mrad
        S = lambda w: 9.*sqrt(3)/8./pi*w*simps(kv(5./3.,linspace(w, 20, num=200)))
        F = lambda eph: mrad*C_fi*self.energy*self.I*eph/self.eph_c*S(eph/self.eph_c)
        return F 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import simps [as 别名]
def radiation_integrals(lattice, twiss_0, nsuperperiod = 1):
    #TODO: add I4 for rectangular magnets I4 = Integrate(2 Dx(z)*k(z)*h(z), Z)
    
    n_points_element = 20
    
    tws_elem = twiss_0
    (I1, I2, I3, I4, I5) = (0., 0., 0., 0., 0.)
    h = 0.
    for elem in lattice.sequence:
        if elem.__class__ in (SBend, RBend, Bend) and elem.l != 0:
            Dx = []
            Hinvariant = []
            Z = []
            h = elem.angle/elem.l

            for z in np.linspace(0, elem.l, num=n_points_element, endpoint=True):
                tws_z = elem.transfer_map(z)*tws_elem
                Dx.append(tws_z.Dx)
                Z.append(z)
                Hx = (tws_z.gamma_x*tws_z.Dx*tws_z.Dx + 2.*tws_z.alpha_x*tws_z.Dxp*tws_z.Dx
                                        + tws_z.beta_x*tws_z.Dxp*tws_z.Dxp)
                Hinvariant.append(Hx)
            #H = array(h)
            H2 = h*h
            H3 = np.abs(h*h*h)
            I1 += h*simps(np.array(Dx), Z)
            I2 += H2*elem.l  #simps(H2, Z)*nsuperperiod
            I3 += H3*elem.l  #simps(H3, Z)*nsuperperiod
            I4 += h*(2*elem.k1 + H2)*simps(np.array(Dx), Z)
            I5 += H3*simps(np.array(Hinvariant), Z)
        tws_elem = elem.transfer_map*tws_elem
    #if abs(tws_elem.beta_x - twiss_0.beta_x)>1e-7 or abs(tws_elem.beta_y - twiss_0.beta_y)>1e-7:
    #    print( "WARNING! Results may be wrong! radiation_integral() -> beta functions are not matching. ")
        #return None
    return (I1*nsuperperiod, I2*nsuperperiod, I3*nsuperperiod, I4*nsuperperiod, I5*nsuperperiod) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import simps [as 别名]
def natural_chromaticity(lattice, tws_0, nsuperperiod = 1):
    #edge_ksi_x, edge_ksi_y = edge_chromaticity(lattice, tws_0)
    edge_ksi_x, edge_ksi_y = 0, 0
    tws_elem = tws_0
    #M = TransferMap()
    integr_x = 0.
    integr_y = 0.
    for elem in lattice.sequence:
        if elem.__class__ in [SBend, RBend, Bend, Quadrupole]:
            bx = []
            by = []
            k = []
            h = []
            Z = []
            for z in np.linspace(0, elem.l,num = 5, endpoint=True):
                twiss_z = elem.transfer_map(z)*tws_elem
                bx.append(twiss_z.beta_x)
                by.append(twiss_z.beta_y)
                k.append(elem.k1)
                #print elem.l
                if elem.__class__ != Quadrupole and elem.l != 0:
                    h.append(elem.angle/elem.l)
                else:
                    h.append(0.)
                Z.append(z)

            H2 = np.array(h)*np.array(h)
            X = np.array(bx)*(np.array(k)+ H2)
            Y = -np.array(by)*np.array(k)
            integr_x += simps(X, Z)
            integr_y += simps(Y, Z)
        elif elem.__class__ == Multipole:
            twiss_z = elem.transfer_map*tws_elem
            integr_x += twiss_z.beta_x*elem.kn[1]
            integr_y -= twiss_z.beta_y*elem.kn[1]
        tws_elem = elem.transfer_map*tws_elem
    ksi_x = -(integr_x - edge_ksi_x)/(4*pi)
    ksi_y = -(integr_y - edge_ksi_y)/(4*pi)
    return np.array([ksi_x*nsuperperiod, ksi_y*nsuperperiod]) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import simps [as 别名]
def LossShape(bunch, wake):
    """
    loss, spread, peak
    dimensions:
                wake - m , Volt/pC
                out - V/pC;
    """
    w = wake[1]
    bi2 = bunch[1]
    s = wake[0]
    loss = simps(-bi2*w, s)
    spread = np.sqrt(simps(bi2*(w + loss)**2, s))
    peak = max(abs(w))
    return loss, spread, peak 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import simps [as 别名]
def pipe_wake(z, current, tube_radius, tube_len, conductivity, tau, roughness, d_oxid):

    Q = simps(current, z)/speed_of_light
    print ("Charge = ", Q*1e12, "pC")
    
    xb = -z[::-1]
    yb = current[::-1]/(speed_of_light*Q)
    Q = Q*1e12 #C->pC

    ds=xb[3]-xb[0]
    xb = np.append(xb, np.arange(1, 100001)*ds + xb[-1])
    yb = np.append(yb, np.arange(1, 100001)*0)

    # roughness and axid layer are here
    eps_r = 2.
    Ind = mu_0*((eps_r-1.)/eps_r*d_oxid + 0.01035*roughness)

    # the result is in V
    W = ResistiveZaZb(xb, yb, tube_radius, conductivity, tau, Ind)#*Q*L

    W = W.real*Q*tube_len
    n = len(current)
    bunch = [xb[:n], yb[:n]]
    wake = [xb[:n], W[:n]]
    # postprocessing
    L, S, P = LossShape(bunch, wake)
    print ('Loss [V]:  ', L)
    print ('Spread [V]:', S)
    print ('Peak [V]:  ', P)
    return bunch, wake 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import simps [as 别名]
def sample_normal_expectations(self, gaussian_state):
        """Returns the expectation value E(f) and the variance var(f)
        for some normal distribution X~N(mu, cov).

        Args:
            mu (array): means vector
            cov (array): covariance matrix
            func (function): function acting on the random variables X, P, XP,
                returning a second order polynomial

        Returns:
            tuple: tuple of expectation value and variance.
        """

        def _sample(func, correction=0, mu=None, cov=None):
            """wrapped function"""
            if mu is None:
                mu = gaussian_state[0]

            if cov is None:
                cov = gaussian_state[1]

            X, P = np.mgrid[-7:7:0.01, -7:7:0.01]
            grid = np.dstack((X, P))
            XP = np.prod(grid, axis=2)

            poly = func(X, P, XP)
            PDF = multivariate_normal.pdf(grid, mu, cov)

            Ex = simps(simps(poly * PDF, P[0]), X.T[0])
            ExSq = simps(simps(poly ** 2 * PDF, P[0]), X.T[0])

            var = ExSq - Ex ** 2 + correction

            return Ex, var

        return _sample