Python scipy.constants.k() Examples

def population(self, Mole, Temp):
        """Calculate population for each level at given temperature"""
        RoTemp = np.reshape(Temp*1., (Temp.size, 1))
        Qr = self.weighti*np.exp(-(self.e_cm1*100*c*h)/(k*RoTemp))
        # RoPr = Qr/Ntotal  # This is for all transitions
        RoPr = Qr/(Qr.sum(axis=1).reshape(RoTemp.size, 1))  # only given trans.
        linet = []
        for xx in range(self.nt):
            gdu, gdl = self.weighti[self.tran_tag[xx][1:].astype(int)]
            _up = int(self.tran_tag[xx][1])
            _low = int(self.tran_tag[xx][2])
            Aei = self.ai[_up, _low]
            line_const = (c*10**2)**2*Aei*(gdu/gdl)*1.e-6*1.e14 /\
                         (8.*np.pi*(self.freq_array[xx]*1.e9)**2)
            # Hz->MHz,cm^2 ->nm^2
            # W = C.h*C.c*E_cm1[_low]*100.  # energy level above ground state
            "This is the function of calculating H2O intensity"
            line = (1.-np.exp(-h*(self.freq_array[xx]*1.e9) /
                              k/RoTemp))*line_const
            linet.append(line[:, 0]*RoPr[:, _low])  # line intensity non-LTE
        Ni_LTE = Mole.reshape((Mole.size, 1))*RoPr  # *0.75  # orth para ratio
        Ite_pop = [[Ni_LTE[i].reshape((self.ni, 1))] for i in range(Mole.size)]
        return Ite_pop

#import numba 

def test_tlcorrection(self):
        # testing long correction function
        x_for_long = np.array([20.,21.,22.,23.,24.,25.])
        y_for_long = np.array([1.0,1.0,1.0,1.0,1.0,1.0])

        nu0 = 1.0/(514.532)*1e9 #laser wavenumber at 514.532
        nu = 100.0*x_for_long # cm-1 to m-1
        T = 23.0+273.15 # the temperature in K

        x_long,long_res,eselong = rp.tlcorrection(x_for_long, y_for_long,23.0,514.532,correction = 'long',normalisation='area') # using the function
        t0 = nu0**3.0*nu/((nu0-nu)**4)
        t1= 1.0 - np.exp(-h*c*nu/(k*T)) # c in m/s  : t1 dimensionless
        long_calc= y_for_long*t0*t1 # pour les y
        long_calc = long_calc/np.trapz(long_calc,x_for_long) # area normalisation

        np.testing.assert_equal(long_res,long_calc)
        np.testing.assert_equal(x_for_long,x_long)

        x_long,long_res,eselong = rp.tlcorrection(x_for_long, y_for_long,23.0,514.532,correction = 'long',normalisation='no') # using the function
        t0 = nu0**3.0*nu/((nu0-nu)**4)
        t1= 1.0 - np.exp(-h*c*nu/(k*T)) # c in m/s  : t1 dimensionless
        long_calc= y_for_long*t0*t1 # pour les y

        np.testing.assert_equal(long_res,long_calc)
        np.testing.assert_equal(x_for_long,x_long) 

def kTCnoiseCsn(temptr, sensecapacity):
    """

        Args:
            | temptr (scalar): temperature in K
            | sensecapacity (): sense node capacitance F
 
        Returns:
            | n (scalar): noise as number of electrons 

        Raises:
            | No exception is raised.
    """
    return np.sqrt(const.k * temptr * sensecapacity) / const.e

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

def FermiDirac(Ef, EJ, T):
    """
    Returns the Fermi-Dirac probability distribution, given the crystal's
    Fermi energy, the temperature and the energy where the distribution values
    is required.

    Args:
        | Ef: Fermi energy in J
        | EJ: Energy in J
        | T : Temperature in K

    Returns:
        | fermiD : the Fermi-Dirac distribution
    """
    #prevent divide by zero
    den = (1 + np.exp( ( EJ - Ef ) / ( T * const.k) ) )
    return 1 / den



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

def Bv_T(Freq, T):
    # brbr = 1  # .e7*1.e-4
    Bv_out = 2.*h*Freq**3/c**2/(np.exp(h*Freq/k/T)-1.)
    return Bv_out 

def _snr_constant(self):
        temp = 290  # noise reference temperature (room temperature kelvin)
        # convert from dB
        noise_figure = 10 ** (self.receiver_noise / 10)
        loss = 10 ** (self.loss / 10)
        # calculate part of snr that is independent of:
        #   rcs, transmitted power, gain and range
        return (const.c ** 2 * self.number_pulses * self.duty_cycle) / \
               (64 * np.pi ** 3 * const.k * temp * self.band_width *
                noise_figure * self.frequency ** 2 * loss) 

def getAverageSpeed(self, temperature):
        """
            Average (mean) speed at a given temperature

            Args:
                temperature (float): temperature (K)

            Returns:
                float: mean speed (m/s)
        """
        return sqrt(8. * C_k * temperature / (pi * self.mass)) 

def getNumberDensity(self, temperature):
        """ Atom number density at given temperature

            See `calculation of basic properties example snippet`_.

            .. _`calculation of basic properties example snippet`:
                ./Rydberg_atoms_a_primer.html#General-atomic-properties

            Args:
                temperature (float): temperature in K
            Returns:
                float: atom concentration in :math:`1/m^3`
        """
        return self.getPressure(temperature) / (C_k * temperature) 

def kTCnoiseGv(temptr, gv):
    """

        Args:
            | temptr (scalar): temperature in K
            | gv (scalar): sense node gain V/e
 
 
        Returns:
            | n (scalar): noise as number of electrons 

        Raises:
            | No exception is raised.
    """
    return np.sqrt(const.k * temptr / (const.e * gv))

############################################################
##
#def 
    """

        Args:
            | (): 
            | (): 
            | (): 
            | (): 
            | (): 
            | (): 
            | (): 
 
        Returns:
            | n (scalar): noise as number of electrons 

        Raises:
            | No exception is raised.
    """


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

def __init__(self):
        """ Precalculate the Planck function constants.

        Reference: http://www.spectralcalc.com/blackbody/appendixC.html
        """

        import scipy.optimize

        self.c1em = 2 * np.pi * const.h * const.c * const.c
        self.c1el = self.c1em * (1.0e6)**(5-1) # 5 for lambda power and -1 for density
        self.c1en = self.c1em * (100)**3 * 100 # 3 for wavenumber, 1 for density
        self.c1ef = 2 * np.pi * const.h / (const.c * const.c)

        self.c1qm = 2 * np.pi * const.c
        self.c1ql = self.c1qm * (1.0e6)**(4-1) # 5 for lambda power and -1 for density
        self.c1qn = self.c1qm * (100)**2 * 100 # 2 for wavenumber, 1 for density
        self.c1nf = 2 * np.pi  / (const.c * const.c)

        self.c2m = const.h * const.c / const.k
        self.c2l = self.c2m * 1.0e6 # 1 for wavelength density
        self.c2n = self.c2m * 1.0e2 # 1 for cm-1 density
        self.c2f = const.h / const.k

        self.sigmae = const.sigma
        self.zeta3 = 1.2020569031595942853
        self.sigmaq = 4 * np.pi * self.zeta3 * const.k ** 3 \
               / (const.h ** 3 * const.c ** 2)

        self.a2 = scipy.optimize.brentq(self.an, 1, 2, (2) )
        self.a3 = scipy.optimize.brentq(self.an, 2, 3, (3) )
        self.a4 = scipy.optimize.brentq(self.an, 3.5, 4, (4) )
        self.a5 = scipy.optimize.brentq(self.an, 4.5, 5, (5) )

        self.wel = 1e6 * const.h * const.c /(const.k * self.a5)
        self.wql = 1e6 * const.h * const.c /(const.k * self.a4)
        self.wen = self.a3 * const.k /(100 * const.h * const.c )
        self.wqn = self.a2 * const.k /(100 * const.h * const.c )
        self.wef = self.a3 * const.k /(const.h )
        self.wqf = self.a2 * const.k /(const.h ) 

def IXV(V, IVbeta, tDetec, iPhoto,I0):
    """
    This function provides the diode curve for a given photocurrent.

    The same function is also used to calculate the dark current, using
    IVbeta=1 and iPhoto=0

    Args:
        | V: bias [V]
        | IVbeta: diode equation non linearity factor;
        | tDetec: detector's temperature [K]
        | iPhoto: photo-induced current, added to diode curve [A]
        | I0: reverse sat current [A]

    Returns:
        | current from detector [A]
    """

    # diode equation from Dereniak's book eq. 7.23
    return I0 * (np.exp(const.e * V / (IVbeta * const.k * tDetec)) - 1) - iPhoto




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

def Absorption(wavelength, Eg, tempDet, a0, a0p):
    """
    Calculate the spectral absorption coefficient
    for a semiconductor material with given material values.

    The model used here is based on Equations 3.5, 3.6 in Dereniaks book.

    Args:
        | wavelength: spectral variable [m]
        | Eg: bandgap energy [Ev]
        | tempDet: detector's temperature in [K]
        | a0: absorption coefficient [m-1] (Dereniak Eq 3.5 & 3.6)
        | a0p:  absorption coefficient in [m-1] (Dereniak Eq 3.5 & 3.6)

    Returns:
        | absorption: spectral absorption coefficient in [m-1]
    """

    #frequency/wavelength expressed as energy in Ev
    E = const.h * const.c / (wavelength * const.e )

    # the np.abs() in the following code is to prevent nan and inf values
    # the effect of the abs() is corrected further down when we select
    # only the appropriate values based on E >= Eg and E < Eg

    # Absorption coef - eq. 3.5- Dereniak
    a35 = (a0 * np.sqrt(np.abs(E - Eg))) + a0p
    # Absorption coef - eq. 3.6- Dereniak
    a36 = a0p * np.exp((- np.abs(E - Eg)) / (const.k * tempDet))
    absorption = a35 * (E >= Eg) + a36 * (E < Eg)

    return absorption

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

def planck(lam, T):
    """ Returns the spectral radiance of a black body at temperature T.

    Returns the spectral radiance, B(lam, T), in W.sr-1.m-2 of a black body
    at temperature T (in K) at a wavelength lam (in nm), using Planck's law.

    """

    lam_m = lam / 1.e9
    fac = h*c/lam_m/k/T
    B = 2*h*c**2/lam_m**5 / (np.exp(fac) - 1)
    return B 

def test_twomode(self):
        """Test if function returns two-mode squeezing parameters that correctly reconstruct the
        input normal mode frequencies."""
        w = -k * self.T / (0.5 * h * c * 100) * np.log(np.tanh(self.t))
        assert np.allclose(w, self.w) 

def energies(samples: list, w: np.ndarray, wp: np.ndarray) -> Union[list, float]:
    r"""Computes the energy :math:`E = \sum_{k=1}^{N}m_k\omega'_k - \sum_{k=N+1}^{2N}n_k\omega_k`
    of each GBS sample in units of :math:`\text{cm}^{-1}`.

    **Example usage:**

    >>> samples = [[1, 1, 0, 0, 0, 0], [1, 2, 0, 0, 1, 1]]
    >>> w  = np.array([300.0, 200.0, 100.0])
    >>> wp = np.array([700.0, 600.0, 500.0])
    >>> energies(samples, w, wp)
    [1300.0, 1600.0]

    Args:
        samples (list[list[int]] or list[int]): a list of samples from GBS, or alternatively a
            single sample
        w (array): normal mode frequencies of initial state in units of :math:`\text{cm}^{-1}`
        wp (array): normal mode frequencies of final state in units of :math:`\text{cm}^{-1}`

    Returns:
        list[float] or float: list of GBS sample energies in units of :math:`\text{cm}^{-1}`, or
        a single sample energy if only one sample is input
    """
    if not isinstance(samples[0], list):
        return np.dot(samples[: len(samples) // 2], wp) - np.dot(samples[len(samples) // 2 :], w)

    return [np.dot(s[: len(s) // 2], wp) - np.dot(s[len(s) // 2 :], w) for s in samples] 

def DopplerWind(Temp, FreqGrid, Para, wind_v, shift_direction='red'):
    u"""#doppler width
    #Para[transient Freq[Hz], relative molecular mass[g/mol]]"""
    # step1 = Para[0]/c*(2.*R*gct/(Para[1]*1.e-3))**0.5
    # outy = np.exp(-(Freq-Para[0])**2/step1**2) / (step1*(np.pi**0.5))
    #wind_v = speed[:,10] 
    #Temp=temp[10]
    #FreqGrid = Fre_range_i[0]
    wind = wind_v.reshape(wind_v.size, 1)
    FreqGrid = FreqGrid.reshape(1, FreqGrid.size)
    deltav = Para[0]*wind/c
    if shift_direction.lower() == 'red':
        D_effect = (deltav)
    elif shift_direction.lower() == 'blue':
        D_effect = (-deltav)
    else:
        raise ValueError('Set shift direction to "red" or "blue".')

#    step1 = Para[0]/c*(2.*R*Temp*np.log(2.)/(Para[1]*1.e-3))**0.5  # HWHM
#    outy = np.exp(-np.log(2.)*(FreqGrid-Para[0])**2/step1**2) *\
#                 (np.log(2.)/np.pi)**0.5/step1
#    outy_d = np.exp(-np.log(2.)*(FreqGrid+D_effect-Para[0])**2/step1**2) *\
#                   (np.log(2.)/np.pi)**0.5/step1
    GD = np.sqrt(2*k*ac/Para[1]*Temp)/c*Para[0]
    step1 = GD
    outy_d = wofz((FreqGrid+D_effect-Para[0])/GD).real / np.sqrt(np.pi) / GD
    #plot(FreqGrid, outy)
    #plot(FreqGrid, outy_d[:,0])
    return outy_d 

def Isaturation(mobE, tauE, mobH, tauH, me, mh, na, nd, Eg, tDetec, areaDet):
    """
    This function calculates the reverse saturation current, by
    Equation 7.22 in Dereniak's book

    Args:
        | mobE: electron mobility [m2/V.s]
        | tauE: electron lifetime [s]
        | mobH: hole mobility [m2/V.s]
        | tauH: hole lifetime [s]
        | me: electron effective mass [kg]
        | mh: hole effective mass [kg]
        | na: acceptor concentration [m-3]
        | nd: donor concentration [m-3]
        | Eg: energy bandgap in [Ev]
        | tDetec: detector's temperature in [K]
        | areaDet: detector's area [m2]

    Returns:
        | I0: reverse sat current [A]
    """

    # diffusion length [m] Dereniak Eq7.19
    Le=np.sqrt(const.k * tDetec * mobE * tauE / const.e)
    Lh=np.sqrt(const.k * tDetec * mobH * tauH / const.e)
    # intrinsic carrier concentration - Dereniak`s book eq. 7.1 - m-3
    # Eg here in eV units, multiply with q
    ni = (np.sqrt(4 * (2 * np.pi * const.k * tDetec / const.h ** 2) ** 3 *\
        np.exp( - (Eg * const.e) / (const.k * tDetec)) * (me * mh) ** 1.5))

    # reverse saturation current - Dereniak eq. 7.22 - Ampère
    if na == 0 or tauE == 0:
        elec = 0
    else:
        elec = Le / (na * tauE)
    if nd == 0 or tauH == 0:
        hole = 0
    else:
        hole = Lh / ( nd * tauH )

    I0 = areaDet * const.e * (ni ** 2) *( elec + hole )

    return (I0,ni)


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

def Noise(tempDet, IVbeta, Isat, iPhoto, vBias=0):
    """
    This function calculates the noise power spectral density produced in the
    diode: shot noise and thermal noise. The assumption is that all noise
    sources are white noise PSD.

    Eq 5.143 plus thermal noise, see Eq 5.148

    Args:
        | tempDet: detector's temperature [K]
        | IVbeta: detector nonideal factor [-]
        | Isat: reverse saturation current [A]
        | iPhoto: photo current [A]
        | vBias: bias voltage on the detector [V]

    Returns:
        | detector noise power spectral density [A/Hz1/2]
        | R0: dynamic resistance at zero bias.
        | Johnson noise only noise power spectral density [A/Hz1/2]
        | Shot noise only noise power spectral density [A/Hz1/2]
    """

    R0 = IVbeta * const.k * tempDet / (Isat * const.e)

    # johnson noise
    iJohnson = 4 * const.k * tempDet / R0

    # shot noise for thermal component Isat
    iShot1 = 2 * const.e * Isat

    # shot noise for thermal component Isat
    iShot2 = 2 * const.e * Isat *np.exp(const.e * vBias /(const.k * tempDet * IVbeta))

    # shot noise for photocurrent
    iShot3 = 2 * const.e * iPhoto

    # total noise
    noise = np.sqrt(iJohnson + iShot1 + iShot2 + iShot3 )

    return noise, R0, np.sqrt(iJohnson), np.sqrt(iShot1 + iShot2 + iShot3 )


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

def printConstants(self):
        """Print Planck function constants.

        Args:
            | None

        Returns:
            | Print to stdout

        Raises:
            | No exception is raised.
        """
        print('h = {:.14e} Js'.format(const.h))
        print('c = {:.14e} m/s'.format(const.c))
        print('k = {:.14e} J/K'.format(const.k))
        print('q = {:.14e} C'.format(const.e))

        print(' ')
        print('pi = {:.14e}'.format(const.pi))
        print('e = {:.14e}'.format(np.exp(1)))
        print('zeta(3) = {:.14e}'.format(self.zeta3 ))
        print('a2 = {:.14e}, root of 2(1-exp(-x))-x'.format(self.a2 ))
        print('a3 = {:.14e}, root of 3(1-exp(-x))-x'.format(self.a3 ))
        print('a4 = {:.14e}, root of 4(1-exp(-x))-x'.format(self.a4 ))
        print('a5 = {:.14e}, root of 5(1-exp(-x))-x'.format(self.a5 ))

        print(' ')
        print('sigmae = {:.14e} W/(m^2 K^4)'.format(self.sigmae))
        print('sigmaq = {:.14e} q/(s m^2 K^3)'.format(self.sigmaq))
        print(' ')
        print('c1em = {:.14e} with wavelenth in m'.format(self.c1em))
        print('c1qm = {:.14e} with wavelenth in m'.format(self.c1qm))
        print('c2m = {:.14e} with wavelenth in m'.format(self.c2m))
        print(' ')
        print('c1el = {:.14e} with wavelenth in $\mu$m'.format(self.c1el))
        print('c1ql = {:.14e} with wavelenth in $\mu$m'.format(self.c1ql))
        print('c2l = {:.14e} with wavelenth in $\mu$m'.format(self.c2l))
        print(' ')
        print('c1en = {:.14e} with wavenumber in cm$^{{-1}}$'.format(self.c1en))
        print('c1qn = {:.14e} with wavenumber in cm$^{{-1}}$'.format(self.c1qn))
        print('c2n = {:.14e} with wavenumber in cm$^{{-1}}$'.format(self.c2n))
        print(' ')
        print('c1ef = {:.14e} with frequency in Hz'.format(self.c1ef))
        print('c1nf = {:.14e} with frequency in Hz'.format(self.c1nf))
        print('c2f = {:.14e} with frequency in Hz'.format(self.c2f))
        print(' ')
        print('wel = {:.14e} um.K  Wien for radiant and wavelength'.format(self.wel))
        print('wql = {:.14e} um.K  Wien for photon rate and wavelength'.format(self.wql))
        print('wen = {:.14e} cm-1/K  Wien for radiant and wavenumber'.format(self.wen))
        print('wqn = {:.14e} cm-1/K  Wien for photon rate and wavenumber'.format(self.wqn))
        print('wef = {:.14e} Hz/K  Wien for radiant and frequency'.format(self.wef))
        print('wqf = {:.14e} Hz/K  Wien for photon rate and frequency'.format(self.wqf))
        print(' ') 

def planck_integral(wavenum, temperature, radiant,niter=512):
    """Integrated Planck law spectral exitance by summation from wavenum to infty.

    http://www.spectralcalc.com/blackbody/inband_radiance.html
    http://www.spectralcalc.com/blackbody/appendixA.html
    Integral of spectral radiance from wavenum (cm-1) to infinity. 
    follows Widger and Woodall, Bulletin of Am Met Soc, Vol. 57, No. 10, pp. 1217

    Args:
        | wavenum (scalar):  spectral limit.
        | temperature (scalar):  Temperature in [K]
        | radiant (bool): output units required, W/m2 or q/(s.m2)
 
    Returns:
        | (scalar):  exitance (not radiance) in units selected.
        | For radiant units will be [W/(m^2)].
        | For not radiant units will be [q/(s.m^2)].

    Raises:
        | No exception is raised, returns None on error.
    """
    # compute powers of x, the dimensionless spectral coordinate
    c1 = const.h * const.c / const.k
    x = c1 * 100. * wavenum / temperature 
    x2 = x * x  
    x3 = x * x2 

    # decide how many terms of sum are needed
    iter = int(2.0 + 20.0 / x )
    iter = iter if iter<niter else niter

    # add up terms of sum 
    sum = 0. 
    for n in range(1,iter):
        dn = 1.0 / n 
        if radiant:
            sum += np.exp(-n * x) * (x3 + (3.0 * x2 + 6.0 * (x + dn) * dn) * dn) * dn
        else:
            sum += np.exp(-n * x) * (x2 + 2.0 * (x + dn) * dn) * dn 

    if radiant:
        # in units of W/m2
        c2 = 2.0 * const.h * const.c * const.c
        rtnval = c2 * (temperature/c1)**4. * sum  
    else:
        # in units of photons/(s.m2)
        kTohc =  const.k * temperature / (const.h * const.c) 
        rtnval =  2.0 * const.c * kTohc**3.  * sum 

    # print('wn={} x={} T={} E={} n={}'.format(wavenum,x,temperature, rtnval/np.pi,iter))

    return rtnval * np.pi

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

def gbs_params(
    w: np.ndarray, wp: np.ndarray, Ud: np.ndarray, delta: np.ndarray, T: float = 0
) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    r"""Converts molecular information to GBS gate parameters.

    **Example usage:**

    >>> formic = data.Formic()
    >>> w = formic.w  # ground state frequencies
    >>> wp = formic.wp  # excited state frequencies
    >>> Ud = formic.Ud  # Duschinsky matrix
    >>> delta = formic.delta  # displacement vector
    >>> T = 0  # temperature
    >>> p = gbs_params(w, wp, Ud, delta, T)

    Args:
        w (array): normal mode frequencies of the initial electronic state in units of
            :math:`\mbox{cm}^{-1}`
        wp (array): normal mode frequencies of the final electronic state in units of
            :math:`\mbox{cm}^{-1}`
        Ud (array): Duschinsky matrix
        delta (array): Displacement vector, with entries :math:`\delta_i=\sqrt{\omega'_i/\hbar}d_i`,
            and :math:`d_i` is the Duschinsky displacement
        T (float): temperature (Kelvin)

    Returns:
        tuple[array, array, array, array, array]: the two-mode squeezing parameters :math:`t`,
        the first interferometer unitary matrix :math:`U_{1}`, the squeezing parameters :math:`r`,
        the second interferometer unitary matrix :math:`U_{2}`, and the displacement
        parameters :math:`\alpha`
    """
    if T < 0:
        raise ValueError("Temperature must be zero or positive")
    if T > 0:
        t = np.arctanh(np.exp(-0.5 * h * (w * c * 100) / (k * T)))
    else:
        t = np.zeros(len(w))

    U2, s, U1 = np.linalg.svd(np.diag(wp ** 0.5) @ Ud @ np.diag(w ** -0.5))
    alpha = delta / np.sqrt(2)

    return t, U1, np.log(s), U2, alpha