Python scipy.constants.pi() Examples

def __plot__(self, orbit, plotter, next_operation=None):
        radius = radius_from_altitude(self.apocenter_altitude, orbit.body)
        if orbit.apocenter_radius > radius:
            label = 'Lowered apocenter'
        else:
            label = 'Raised apocenter'
        self.__apply__(orbit)

        with saved_state(orbit):
            if (next_operation is not None and
                    isinstance(next_operation, TimeOperation)):
                orbit.apply_maneuver(next_operation)
                f2 = orbit.f
                if f2 == 0:
                    f2 = 2 * pi
            else:
                f2 = 2 * pi

        plotter._plot_orbit(orbit, f1=0, f2=f2, label=label) 

def __init__(self, axes=None, num_points=100):
        if axes:
            self.fig = axes.get_figure()
        else:
            self.fig = plt.figure()
            axes = self.fig.add_subplot(111, projection='3d')
        self.axes = axes
        self.axes.set_xlabel("$x$ [km]")
        self.axes.set_ylabel("$y$ [km]")
        self.axes.set_zlabel("$z$ [km]")

        # These are used to fix aspect ratio of final plot.
        # See Plotter3D._force_aspect()
        self._coords_x = np.array(0)
        self._coords_y = np.array(0)
        self._coords_z = np.array(0)

        self.points_per_rad = num_points / (2 * pi) 

def from_state_vector(cls, r, v, body, ref_epoch=J2000):
        """Create orbit from given state vector."""
        elements = elements_from_state_vector(r, v, body.mu)

        self = cls(
            a=elements.a,
            e=elements.e,
            i=elements.i,
            raan=elements.raan,
            arg_pe=elements.arg_pe,
            M0=mean_anomaly_from_true(elements.e, elements.f),
            body=body,
            ref_epoch=ref_epoch)

        # Fix mean anomaly at epoch for new orbit and position.
        oldM0 = self.M0
        self.M0 = ou.mod(self.M - self.n * self.t, 2 * pi)
        assert self.M0 == oldM0

        return self 

def entr_mixed(x, s, shift, delta_0, act_s1, fit_param_fe):
    """
    Returns a fit function for entropies based on the arctan function and the dilute species model fit of SrFeOx
    (see docstring in Atanfit.entropy)
    :param x:               absolute delta
    :param s:               slope, measure for the preference of B species reduction over A species reduction
    :param shift:           constant entropy shift
    :param delta_0:         shift from absolute delta
    :return:                dS of solid solution at delta = x with delta_0
    """
    efe = entr_fe(x + delta_0, fit_param_fe)
    return (
        ((act_s1 * efe) / pi) * (pd.np.arctan((x - delta_0) * s) + pi / 2)
        + (1 - act_s1) * efe
        + shift
    ) 

def test_set_apside_altitudes(self):
        set_apocenter = SetApocenterAltitudeTo(1000 * kilo)
        self.LEO.propagate_anomaly_to(M=0)
        self.LEO.apply_maneuver(set_apocenter)
        self.assertAlmostEqual(
            self.LEO.apocenter_radius,
            self.LEO.body.mean_radius + set_apocenter.apocenter_altitude)
        self.assertTrue(self.LEO.e > 0)

        set_pericenter = SetPericenterAltitudeTo(150 * kilo)
        self.LEO.propagate_anomaly_to(M=pi)
        self.LEO.apply_maneuver(set_pericenter)
        self.assertAlmostEqual(
            self.LEO.pericenter_radius,
            self.LEO.body.mean_radius + set_pericenter.pericenter_altitude)
        self.assertAlmostEqual(
            self.LEO.apocenter_radius,
            self.LEO.body.mean_radius + set_apocenter.apocenter_altitude) 

def __plot__(self, orbit, plotter, next_operation=None):
        if orbit.apocenter_radius > self.apocenter_radius:
            label = 'Lowered apocenter'
        else:
            label = 'Raised apocenter'
        self.__apply__(orbit)

        with saved_state(orbit):
            if (next_operation is not None and
                    isinstance(next_operation, TimeOperation)):
                orbit.apply_maneuver(next_operation)
                f2 = orbit.f
                if f2 == 0:
                    f2 = 2 * pi
            else:
                f2 = 2 * pi

        plotter._plot_orbit(orbit, f1=0, f2=f2, label=label) 

def __plot__(self, orbit, plotter, next_operation=None):
        if orbit.pericenter_radius > self.pericenter_radius:
            label = 'Lowered pericenter'
        else:
            label = 'Raised pericenter'
        self.__apply__(orbit)

        with saved_state(orbit):
            if (next_operation is not None and
                    isinstance(next_operation, TimeOperation)):
                orbit.apply_maneuver(next_operation)
                f2 = orbit.f
                if almost_equal(f2, pi):
                    f2 = pi + 2 * pi
            else:
                f2 = pi + 2 * pi

        plotter._plot_orbit(orbit, f1=pi, f2=f2, label=label) 

def __plot__(self, orbit, plotter, next_operation=None):
        radius = radius_from_altitude(self.pericenter_altitude, orbit.body)
        if orbit.pericenter_radius > radius:
            label = 'Lowered pericenter'
        else:
            label = 'Raised pericenter'
        self.__apply__(orbit)

        with saved_state(orbit):
            if (next_operation is not None and
                    isinstance(next_operation, TimeOperation)):
                orbit.apply_maneuver(next_operation)
                f2 = orbit.f
                if almost_equal(f2, pi):
                    f2 = pi + 2 * pi
            else:
                f2 = pi + 2 * pi

        plotter._plot_orbit(orbit, f1=0, f2=f2, label=label) 

def __plot__(self, orbit, plotter, next_operation=None):
        if self.delta < 0:
            label = 'Lowered pericenter'
        else:
            label = 'Raised pericenter'
        self.__apply__(orbit)

        with saved_state(orbit):
            if (next_operation is not None and
                    isinstance(next_operation, TimeOperation)):
                orbit.apply_maneuver(next_operation)
                f2 = orbit.f
                if almost_equal(f2, pi):
                    f2 = pi + 2 * pi
            else:
                f2 = pi + 2 * pi

        plotter._plot_orbit(orbit, f1=0, f2=f2, label=label) 

def velocity_delta(self, orbit):
        with saved_state(orbit):
            if self.raise_pericenter:
                orbit.propagate_anomaly_to(M=pi)
                radius = orbit.apocenter_radius
            else:
                orbit.propagate_anomaly_to(M=0)
                radius = orbit.pericenter_radius

            old_velocity = orbit.v

            a, e = elements_for_apsides(radius, radius)
            orbit.a = a
            orbit.e = e

            new_velocity = orbit.v

            return new_velocity - old_velocity 

def shiftmp(freq, s11, shiftplanes):#{{{
    """ Adjusts the reflection phase like if the monitor planes were not centered.

    For symmetric metamaterial cell, this function is not needed. The symmetry requires that
    the monitor planes in front of and behind the mm cell are centered.

    However, for an asymmetric metamaterial, the correct position has to be found. Otherwise
    the Fresnel inversion gives negative imaginary part of N and/or negative real part of Z, which
    is quite spooky for passive medium. 
    
    Even such metamaterials, however, may be properly homogenized if we define the 
    position of monitor planes as a function of frequency. We can assume that:
    1) This optimum shift shall hold for all simulations with one or more unit cellnumber.
    2) When the wave direction is reversed (to get s21, s22 parameters), the shift should be negated.
    These rules should enable us to homogenize any asymmetric non-chiral metamaterial. 

    Note that this shifting is still an experimental technique and has to be tested out thoroughly. 
    """
    return np.array(s11) * np.exp(1j*np.array(shiftplanes)/(c/freq) * 2*pi * 2)
#}}} 

def PerdidaPresion(self):
        entrada = self.kwargs["entrada"]
        rhoS = entrada.solido.rho
        rhoG = entrada.Gas.rho
        QS = entrada.solido.caudal
        Q = entrada.caudalmasico

        if self.kwargs["modelo_DeltaP"] == 0:
            # Cinetic loss
            DeltaP = Pressure(self.kf/2.*rhoG*self.V**2)
        elif self.kwargs["modelo_DeltaP"] == 1:
            # Casal-Martinez-Benet
            DeltaP = Pressure(0.003*self.N*rhoG*self.V**2, "inH2O")
        elif self.kwargs["modelo_DeltaP"] == 2:
            # Leith-Licht
            DeltaP = Pressure(0.003*rhoG*16*(entrada.Q/self.num_ciclones)**2 /
                              self.De**2/self.Bc/self.Hc, "mmHg")
        else:
            # Sheferd, Lapple y Ter Linden
            Ae = self.Bc*self.Hc
            As = pi/4*self.De**2
            rhom = rhoG+QS/rhoG/(Q+QS/rhoG)*(rhoS-rhoG)
            DeltaP = Pressure(1.078*(Ae/As)**1.21*rhom*self.V**2, "mmH2O")
        return DeltaP 

def velocidad_f_presion(self):
        entrada = self.kwargs["entrada"]
        rhoS = entrada.solido.rho
        rhoG = entrada.Gas.rho
        dP = self.DeltaPAdmisible
        QS = entrada.solido.caudal
        Q = entrada.caudalmasico

        if self.kwargs["modelo_DeltaP"] == 0:
            velocidad = sqrt(dP*2/self.kf/rhoG)
        elif self.kwargs["modelo_DeltaP"] == 1:
            velocidad = sqrt(dP.inH2O/self.N/0.003/rhoG)
        elif self.kwargs["modelo_DeltaP"] == 2:
            velocidad = sqrt(dP*self.De**2*self.Bc*self.Hc/0.003/rhoG/16)
        else:
            Ae = self.Bc*self.Hc
            As = pi/4*self.De**2
            rhom = rhoG+QS/rhoG/(Q+QS/rhoG)*(rhoS-rhoG)
            velocidad = sqrt(dP.mmH2O/1.078/(Ae/As)**1.21/rhom)
        return velocidad 

def _visco0(self, rho, T, fase=None):
        a = [17.67484, -2.78751, 311498.7, -48826500, 3938774000, -1.654629e11,
             2.86561e12]
        Tr = T/0.29944
        y = 0.68321*(a[0] + a[1]*log10(Tr) + a[2]/Tr**2 + a[3]/Tr**3 +
                     a[4]/Tr**4 + a[5]/Tr**5 + a[6]/Tr**6)
        nt = 266.93*(T*self.M)**0.5/y
        om = rho/1673.0
        c = [1.03010, -0.99175, 2.47127, -3.11864, 1.57066]
        b = [0.48148, -1.18732, 2.80277, -5.41058, 7.04779, -3.76608]
        sum1 = sum([ci*om**i for i, ci in enumerate(c)])
        sum2 = sum([bi*om**i for i, bi in enumerate(b)])
        sigma = 3.05e-10*(sum1-sum2*log10(T/122.1))
        br = 2.0/3.0*pi*Avogadro*sigma**3
        brho = rho/self.M*1000*br
        d = [1, 0.27676, 0.014355, 2.6480, -1.9643, 0.89161]
        nd = sum([di*brho**i for i, di in enumerate(d)])
        return unidades.Viscosity(nd*nt/100, "muPas") 

def _areaFzFull(self):
        """
        area of z-faces prior to deflation
        """
        if self.isSymmetric:
            return np.kron(
                    np.ones_like(self.vectorNz), pi*(
                        self.vectorNx**2 -
                        np.r_[0, self.vectorNx[:-1]]**2
                    )
            )
        return np.kron(
            np.ones(self._ntNz), np.kron(
                self.hy,
                0.5 * (self.vectorNx[1:]**2 - self.vectorNx[:-1]**2)
            )
        ) 

def get_interaction_constant(accelerating_voltage):
    """Calculates the interaction constant, sigma, for a given
    acelerating voltage.

    Parameters
    ----------
    accelerating_voltage : float
        The accelerating voltage in V.

    Returns
    -------
    sigma : float
        The relativistic electron wavelength in m.

    """
    E = accelerating_voltage
    wavelength = get_electron_wavelength(accelerating_voltage)
    sigma = 2 * pi * (m_e + e * E)

    return sigma 

def _getRadialDipoleSemiClassical(self, n1, l1, j1, n2, l2, j2,
                                      s=0.5):

        # get the effective principal number of both states
        nu = np.sqrt( - self.scaledRydbergConstant
                     / self.getEnergy(n1, l1, j1, s=s))
        nu1 = np.sqrt( - self.scaledRydbergConstant
                      / self.getEnergy(n2, l2, j2, s=s))

        # get the parameters required to calculate the sum
        l_c = (l1 + l2 + 1.) / 2.
        nu_c = sqrt(nu * nu1)

        delta_nu = nu - nu1
        delta_l = l2 - l1

        # I am not sure if this correct

        gamma = (delta_l * l_c) / nu_c

        if delta_nu == 0:
            g0 = 1
            g1 = 0
            g2 = 0
            g3 = 0
        else:

            g0 = (1. / (3. * delta_nu)) * (
                angerj(delta_nu - 1., - delta_nu)
                - angerj(delta_nu + 1, - delta_nu))
            g1 = -(1. / (3. * delta_nu)) * (
                angerj(delta_nu - 1., - delta_nu)
                + angerj(delta_nu + 1, -delta_nu))
            g2 = g0 - np.sin(np.pi * delta_nu) / (np.pi * delta_nu)
            g3 = (delta_nu / 2.) * g0 + g1

        radial_ME = (3 / 2) * nu_c**2 * (1 - (l_c / nu_c)**(2))**0.5 * \
            (g0 + gamma * g1 + gamma**2 * g2 + gamma**3 * g3)
        return float(radial_ME) 

def getC3term(self, n, l, j, n1, l1, j1, n2, l2, j2, s=0.5):
        """
            C3 interaction term for the given two pair-states

            Calculates :math:`C_3` intaraction term for
                :math:`|n,l,j,n,l,j\\rangle \
                 \\leftrightarrow |n_1,l_1,j_1,n_2,l_2,j_2\\rangle`

            Args:
                n (int): principal quantum number
                l (int): orbital angular momenutum
                j (float): total angular momentum
                n1 (int): principal quantum number
                l1 (int): orbital angular momentum
                j1 (float): total angular momentum
                n2 (int): principal quantum number
                l2 (int): orbital angular momentum
                j2 (float): total angular momentum
                s (float): optional, total spin angular momentum of state.
                    By default 0.5 for Alkali atoms.

            Returns:
                float:  :math:`C_3 = \\frac{\\langle n,l,j |er\
                |n_1,l_1,j_1\\rangle \
                \\langle n,l,j |er|n_2,l_2,j_2\\rangle}{4\\pi\\varepsilon_0}`
                (:math:`h` Hz m :math:`{}^3`).
        """
        d1 = self.getRadialMatrixElement(n, l, j, n1, l1, j1, s=s)
        d2 = self.getRadialMatrixElement(n, l, j, n2, l2, j2, s=s)
        d1d2 = 1 / (4.0 * pi * epsilon_0) * d1 * d2 * C_e**2 *\
            (physical_constants["Bohr radius"][0])**2
        return d1d2 

def _BlochFunction(self, x, stateVector, q, k=1.):
        r"""
        Bloch wavefunctions

        Args:
            x (x): position (in units \2 pi/k, for default value of laser
                wavevector unit k=1, one full wavelength is 2\pi)
            stateVector: eigen vector obtained by diagonalisation of
                interaction Hamiltonian in a subspace given by the selected
                quasimomentum
            q (float): quasimomentum (in units of driving laser k)
            k (float): driving laser wavevector, define units for momentum and
                distance;
                if k==1 (default value), reciprocal lattice momentum is 2,
                and the full range of quasimomentum is from -1 to +1;
                one full wavelength is the 2\pi.

        Retruns:
            float:
        """

        index = len(stateVector) // 2 + 2  # Align Bloch functions in phase
        angle = np.angle(stateVector[index])
        sign = np.exp(-1j*angle)
        temp = 0 + 0j
        for l in np.arange(-self.lLimit, self.lLimit + 1, 1):
            temp += sign * stateVector[l + self.lLimit] \
                * np.exp(1.j * (2. * k * l + q) * x)
        return temp 

def _thermo0(self, rho, T, fase):
        # λ1 in Eq 3 is always 0

        GT = [-2.903423528e5, 4.680624952e5, -1.8954783215e5, -4.8262235392e3,
              2.243409372e4, -6.6206354818e3, 8.9937717078e2, -6.0559143718e1,
              1.6370306422]
        lo = 0
        for i in range(-3, 6):
            lo += GT[i+3]*T**(i/3.)

        l2, lc = 0, 0
        if rho:
            tita = (rho-221)/221
            k = [-1.304503323e1, 1.8214616599e1, -9.903022496e3, 7.420521631e2,
                 -3.0083271933e-1, 9.6456068829e1, 1.350256962e4]
            l2 = exp(k[0]+k[3]/T) * (
                exp(rho.gcc**0.1*(k[1]+k[2]/T**1.5) +
                    tita*rho.gcc**0.5*(k[4]+k[5]/T+k[6]/T**2))-1)

            # Critical enhancement
            deltarho = (rho-221)/221
            deltaT = (T-282.34)/282.34

            xt = rho**2*fase.kappa*5.039/221**2
            B = abs(deltarho)/abs(deltaT)**1.19                         # Eq 11
            Gamma = xt*abs(deltaT)**1.19                                # Eq 12
            xi = 0.69/(B**2*5.039/Gamma/Boltzmann/282.34)**0.5          # Eq 14

            # Eq 19
            F = exp(-18.66*deltaT**2) * exp(-4.25*deltarho**4)

            # Eq 18
            c = (self.M/rho.gcc/Avogadro/Boltzmann/T)**0.5
            d = Boltzmann*T**2/6/pi/fase.mu.muPas/xi
            lc = c*d*fase.dpdT_rho**2*fase.kappa**0.5*F

        return unidades.ThermalConductivity(lo+l2+lc, "mWmK") 

def calcularRendimientos_parciales(self, A):
        """Calculate the separation efficiency per diameter"""
        entrada = self.kwargs["entrada"]
        rendimiento_parcial = []
        for dp in entrada.solido.diametros:
            if dp <= 1e-6:
                q = dp*e*1e8
            else:
                q = pi*epsilon_0*self.potencialCarga*dp**2 * \
                    (1+2*(self.epsilon-1.)/(self.epsilon+2.))
            U = q*self.potencialDescarga/(3*pi*dp*entrada.Gas.mu)
            rendimiento_parcial.append(Dimensionless(1-exp(-U*A/entrada.Q)))
        return rendimiento_parcial 

def epoch(self, value):
        """Set epoch, adjusting current mean anomaly (from which
        other anomalies are calculated).
        """
        t = (value - self.ref_epoch).sec
        self._M = self.M0 + self.n * t
        self._M = ou.mod(self._M, 2 * pi)
        self._t = t 

def t(self, value):
        """Set time since ref_epoch, adjusting current mean anomaly (from which
        other anomalies are calculated).
        """
        self._M = self.M0 + self.n * value
        self._M = ou.mod(self._M, 2 * pi)
        self._t = value 

def nz2rt(freq, N, Z, d):#{{{
    """ Returns the complex reflection and transmission parameters for a metamaterial slab.

    Useful for reverse calculation of eps and mu (to check results)
    
    Accepts: 
    * frequency array,
    * effective refractive index N, 
    * effective impedance Z, 
    * vacuum wave vector k and 
    * thickness d of the layer.
    """
    ## Direct derivation from infinite sum of internal reflections
    k = 2*pi * freq/c          # the wave vector
    t1 = 2 / (1+Z)              # transmission of front interface
    t2 = 2*Z / (Z+1)            # transmission of back interface
    t1prime = Z*t1       
    r1=(Z-1)/(Z+1)              # reflection of front interface
    r2=(1-Z)/(1+Z)              # reflection of back interface
    s12 = t1*t2*np.exp(1j*k*N*d) / (1 + r1*r2*np.exp(2j*k*N*d))
    s11 = r1 + t1prime*t1*r2*np.exp(2j*k*N*d)/(1+r1*r2*np.exp(2j*k*N*d))
    return s11, s12
    """
    Note: these results may be also re-expressed using goniometric functions.
    Equations from Smith2002 or Cai-Shalaev, mathematically equivalent to those above 
    (only Smith's s11 has negative sign convention).
    s12new = 1/(np.cos(N*k*d) - .5j*(Z+1/Z)*np.sin(N*k*d)) 
    s11new = -s12new * .5j*(Z-1/Z)*np.sin(N*k*d)      

    TODO: implement also for other surrounding media than vacuum.
    """
#}}}

## --- Preparation of data ------------------------------------ # {{{
## Get reflection and transmission data, prepare its parameters 

def M(self, value):
        warnings.warn('Setting anomaly does not set time, use KeplerianElements'
                      '.propagate_anomaly_to() instead.', OrbitalWarning)
        self._M = ou.mod(value, 2 * pi) 

def _plot_body(self, orbit):
        color = '#EBEBEB'
        if orbit.body.plot_color is not None:
            color = orbit.body.plot_color

        u = np.linspace(0, 2 * np.pi, 100)
        v = np.linspace(0, np.pi, 50)
        cx = orbit.body.mean_radius * np.outer(np.cos(u), np.sin(v))
        cy = orbit.body.mean_radius * np.outer(np.sin(u), np.sin(v))
        cz = orbit.body.mean_radius * np.outer(np.ones(np.size(u)), np.cos(v))
        cx, cy, cz = cx / kilo, cy / kilo, cz / kilo
        self.axes.plot_surface(cx, cy, cz, rstride=5, cstride=5, color=color,
                               edgecolors='#ADADAD', shade=False) 

def naive_hilbert_transform(x, y, new_x): ## or, just a discrete convolution with the 1/t function
    old_x_grid, new_x_grid = np.meshgrid(x, new_x)
    sharpness = 5000         # with ideally dense sampled data, this should converge to infinity; reduce it to avoid ringing 
    return -1j * np.sum(y * np.arctan(1/(new_x_grid - old_x_grid)/sharpness)*sharpness, axis=1) / len(x) / (2*pi) 

def top_tick_function(p): return p/(cellsize/2/np.pi) / 100 

def eccentric_anomaly_from_true(e, f):
    """Convert true anomaly to eccentric anomaly."""
    E = atan2(sqrt(1 - e ** 2) * sin(f), e + cos(f))
    E = mod(E, 2 * pi)
    return E 

def _edgeEyFull(self):
        """
        full vector of y-edge lengths (prior to deflating)
        """
        if self.isSymmetric:
            return 2*pi*self.gridN[:, 0]
        return np.kron(
            np.ones(self._ntNz),
            np.kron(self.hy, self.vectorNx)
        )