Python scipy.constants.epsilon_0() Examples

def permittivity2conductivity(complex_eps, freq):#{{{
    """
    Enables to use the same dispersive materials for time- and frequency-domain simulation

    Complex permittivity can express also the conductivity of the sample (in the same
    manner as dielectric losses) with the corresponding relation:
        complex_eps = real_eps - 1j conductivity / (frequency * 2*pi * epsilon_0)
    Therefore it should be inverted for classic D-conductivity:
        conductivity = -
    In order to simulate any lossy medium with the freq-domain solver, we invert this relation
    to obtain a (nondispersive) conductivity for one frequency. But it does not give the same results
    as time-domain simulation.

        What we know: 
           (we know that     c**.5 * np.pi  = 54395 ) 
        function of c, f, 2pi, eps0, eps.im/eps.r
        should give dimension 1  to feed unitless meep
        should be proportional to eps.im/eps.r
        should give ca. 50000 for omega = 2pi * 800 GHz and eps.im/eps.r=0.02
           => should give 2.5e6 for (eps.im/eps.r=1)
        should be proportional to frequency omega
            => should give 5e-7 for omega = 1 and  (eps.im/eps.r=1)
        already was pre-divided for meep by c = 3e8  (which counts here)
            => in real life it shall be  3e-6 * 3e8 = 148
        should be proportional to epsilon0 [C/Vsm], which is similar in units to conductivity
            => if epsilon0 was 1, it would be 1.7e13 -> roughly c**2
    """
    ## TODO resolve the Bulgarian constant when running freq-domain simulation
    # return complex_eps.imag * freq * 2*np.pi * epsilon_0 * complex_eps.real  ## orig. idea
    # return complex_eps.imag * freq * 2*np.pi * epsilon_0 * complex_eps.real ## also wrong
    # return complex_eps.imag / complex_eps.real * 2*np.pi * c
    #return complex_eps.imag / complex_eps.real * 6.28*freq * 8.85e-12 * c
    magic_constant = 1.65e13       ## A. K. A. bulgarian constant...
    return complex_eps.imag / complex_eps.real * 2 * np.pi * freq * epsilon_0 * c**.5 * np.pi
#}}} 

def __init__(self, where=None, lfconductivity=15e3, f_c=1e14):
        """
        At microwave frequencies (where omega << gamma), the most interesting property of the metal 
        is its conductivity sigma(omega), which may be approximated by a nearly constant real value:
            
            conductivity = epsilon/1j * frequency * eps0 * 2pi

        Therefore, for convenience, we allow the user to specify the low-frequency conductivity sigma0 and compute
        the scattering frequency as:
        """
        #self.gamma = omega_p**2 * epsilon_0 / lfconductivity
        """
        Values similar to those of aluminium are used by default:
            dielectric part of permittivity: 1.0, 
            plasma frequency = 3640 THz
            low-frequency conductivity: 40e6 [S/m]
        """
        ## Design a new metallic model
        ## We need to put the scattering frequency below fc, so that Re(eps) is not constant at f_c
        self.gamma = .5 * f_c           

        ## The virtual plasma frequency is now determined by lfconductivity
        f_p = (self.gamma * lfconductivity / (2*pi) / epsilon_0)**.5       
        self.eps = (f_p/f_c)**2  ## add such an epsilon value, that just shifts the permittivity to be positive at f_c
        #self.eps = 1 ##XXX
        print "F_C", f_c
        print "GAMMA", self.gamma
        print "F_P", f_p
        print "LFC", lfconductivity
        print "eps", self.eps
        print


        ## Feed MEEP with a (fake) Drude model
        f_0 = 1e5           ## arbitrary low frequency that makes Lorentz model behave as Drude model
        self.pol = [
                {'omega': f_0, 'gamma': self.gamma, 'sigma': f_p**2 / f_0**2}, # (Lorentz) model
                ]
        self.name = "Drude metal for <%.2g Hz" % f_c
        self.where = where
#}}} 

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 __init__(self):
        self.n = 2
        f = 1000000      # f is a scalar, it's the trap frequency
        self.w = 2 * pi * f
        self.C = (4 * pi * constants.epsilon_0) ** (-1) * constants.e ** 2
        # C is a scalar, it's the I
        self.m = 39.96 * 1.66e-27 

def getRabiFrequency(self,
                         n1, l1, j1, mj1,
                         n2, l2, j2, q,
                         laserPower, laserWaist,
                         s=0.5):
        """
            Returns a Rabi frequency for resonantly driven atom in a
            center of TEM00 mode of a driving field

            Args:
                n1,l1,j1,mj1 : state from which we are driving transition
                n2,l2,j2 : state to which we are driving transition
                q : laser polarization (-1,0,1 correspond to :math:`\\sigma^-`,
                    :math:`\\pi` and :math:`\\sigma^+` respectively)
                laserPower : laser power in units of W
                laserWaist : laser :math:`1/e^2` waist (radius) in units of m
                s (float): optional, total spin angular momentum of state.
                    By default 0.5 for Alkali atoms.

            Returns:
                float:
                    Frequency in rad :math:`^{-1}`. If you want frequency
                    in Hz, divide by returned value by :math:`2\\pi`
        """
        maxIntensity = 2 * laserPower / (pi * laserWaist**2)
        electricField = sqrt(2. * maxIntensity / (C_c * epsilon_0))
        return self.getRabiFrequency2(n1, l1, j1, mj1,
                                      n2, l2, j2, q,
                                      electricField,
                                      s=s) 

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 __init__(self, where=None, lfconductivity=15e3, f_c=1e14, gamma_factor = .5, epsplus=0):
        """
        At microwave frequencies (where omega << gamma), the most interesting property of the metal 
        is its conductivity sigma(omega), which may be approximated by a nearly constant real value:
            
            conductivity = epsilon/1j * frequency * eps0 * 2pi

        Therefore, for convenience, we allow the user to specify the low-frequency conductivity sigma0 and compute
        the scattering frequency as:

        Values similar to those of aluminium are used by default:
            dielectric part of permittivity: 1.0, 
            plasma frequency = 3640 THz
            low-frequency conductivity: 40e6 [S/m]

        The gamma_factor is required for numerical stability. It has always to be lower than the timestepping frequency f_c, 
        but sometimes the simulation is unstable even though and then it has to be reduced e.g. to 1e-2 or to even smaller value.
        """

        ## Design a new metallic model. We use 'omega_p' and 'gamma' in angular units, as is common in textbooks
        ## We need to put the scattering frequency below fc, so that Re(eps) is not constant at f_c (it would hinder
        ## our trick that ensures FDTD stability)
        self.gamma      = f_c * gamma_factor * 2*np.pi
        ## TODO: test when unstable -> reduce self.gamma by 100 -> test again -> finally write automatic selection of value

        ## Knowing the scattering frequency gamma, the virtual plasma frequency is now determined by lfconductivity
        omega_p = np.sqrt(self.gamma * lfconductivity / epsilon_0)

        ## The following step is required for FDTD stability: 
        ## Add such an high-frequency-epsilon value that shifts the permittivity to be positive at f_c
        ## If self.eps>1, the frequency where permittivity goes positive will generally be less than f_p
        self.eps = max((omega_p/2/np.pi / f_c)**2,   1.) + epsplus
        #print 'self.eps', self.eps
        #print "gamma =", self.gamma
        #print 'omega_p = ', omega_p
        #print 'LFC = %.3e' % (omega_p**2 * epsilon_0 / self.gamma)

        ## Feed MEEP with a Lorentz oscillator of arbitrarily low frequency f_0 so that it behaves as the Drude model
        omega_0 = .1           
        self.pol = [
                {'omega': omega_0/(2*np.pi), 'gamma': self.gamma/(2*np.pi), 'sigma': (omega_p/omega_0)**2}, # (Lorentz) model
                ## Note: meep also uses keywords 'omega' and 'gamma', but they are non-angular units
                ]
        self.name = "Drude metal for up to %.2g Hz" % f_c
        self.where = where
#}}}

## -- Prepared for FDTD simulations in the terahertz range (obsoleted) -- 

def _atomLightAtomCoupling(n, l, j, nn, ll, jj, n1, l1, j1, n2, l2, j2,
                           atom1, atom2=None, s=0.5, s2=None):
    """
        Calculates radial part of atom-light coupling

        This function might seem redundant, since similar function exist for
        each of the atoms. Function that is not connected to specific
        atomic species is provided in order to provides route to implement
        inter-species coupling.
    """
    if atom2 is None:
        # if not explicitly inter-species, assume it's the same species
        atom2 = atom1
    if s2 is None:
        s2 = s

    # determine coupling
    dl = abs(l - l1)
    dj = abs(j - j1)
    c1 = 0
    if dl == 1 and (dj < 1.1):
        c1 = 1  # dipole couplings1
    elif (dl == 0 or dl == 2 or dl == 1) and(dj < 2.1):
        c1 = 2  # quadrupole coupling
    else:
        return False
    dl = abs(ll - l2)
    dj = abs(jj - j2)
    c2 = 0
    if dl == 1 and (dj < 1.1):
        c2 = 1  # dipole coupling
    elif (dl == 0 or dl == 2 or dl == 1) and(dj < 2.1):
        c2 = 2  # quadrupole coupling
    else:
        return False

    radial1 = atom1.getRadialCoupling(n, l, j, n1, l1, j1, s=s)
    radial2 = atom2.getRadialCoupling(nn, ll, jj, n2, l2, j2, s=s2)

    # TO-DO: check exponent of the Boht radius (from where it comes?!)

    coupling = C_e**2 / (4.0 * pi * epsilon_0) * radial1 * radial2 *\
        (physical_constants["Bohr radius"][0])**(c1 + c2)
    return coupling


# ================== Saving and loading calculations (START) ==================