Python scipy.constants.h() Examples

def get_electron_wavelength(accelerating_voltage):
    """Calculates the (relativistic) electron wavelength in Angstroms for a
    given accelerating voltage in kV.

    Parameters
    ----------
    accelerating_voltage : float or 'inf'
        The accelerating voltage in kV. Values `numpy.inf` and 'inf' are
        also accepted.

    Returns
    -------
    wavelength : float
        The relativistic electron wavelength in Angstroms.

    """
    if accelerating_voltage in (np.inf, "inf"):
        return 0
    E = accelerating_voltage * 1e3
    wavelength = (
        h / math.sqrt(2 * m_e * e * E * (1 + (e / (2 * m_e * c * c)) * E)) * 1e10
    )
    return wavelength 

def plotSpectre(transitions, eneval, spectre):
    """ plot the UV-visible spectrum using matplotlib. Absissa are converted in nm. """

    # lambda in nm
    lambdaval = [cst.h * cst.c / (val * cst.e) * 1.e9 for val in eneval]

    # plot gaussian spectra
    plt.plot(lambdaval, spectre, "r-", label = "spectre")

    # plot transitions
    plt.vlines([val[1] for val in transitions], \
               0., \
               [val[2] for val in transitions], \
               color = "blue", \
               label = "transitions" )

    plt.xlabel("lambda   /   nm")
    plt.ylabel("Arbitrary unit")
    plt.title("UV-visible spectra")
    plt.grid()
    plt.legend(fancybox = True, shadow = True)
    plt.show() 

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 Responsivity(wavelength, quantumEffic):
    """
    Responsivity quantifies the amount of output seen per watt of radiant
    optical power input [1]. But, for this application it is interesting to
    define spectral responsivity that is the output per watt of monochromatic
    radiation.

    The model used here is based on Equations 7.114 in Dereniak's book.

    Args:
        | wavelength: spectral variable [m]
        | quantumEffic: spectral quantum efficiency

    Returns:
        | responsivity in [A/W]
    """

    return (const.e * wavelength * quantumEffic) / (const.h * const.c)


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

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 get_unique_families(hkls):
    """Returns unique families of Miller indices, which must be permutations of
    each other.

    Parameters
    ----------
    hkls : list
        List of Miller indices ([h, k, l])

    Returns
    -------
    pretty_unique : dictionary
        A dict with unique hkl and multiplicity {hkl: multiplicity}.
    """

    def is_perm(hkl1, hkl2):
        h1 = np.abs(hkl1)
        h2 = np.abs(hkl2)
        return all([i == j for i, j in zip(sorted(h1), sorted(h2))])

    unique = collections.defaultdict(list)
    for hkl1 in hkls:
        found = False
        for hkl2 in unique.keys():
            if is_perm(hkl1, hkl2):
                found = True
                unique[hkl2].append(hkl1)
                break
        if not found:
            unique[hkl1].append(hkl1)

    pretty_unique = {}
    for k, v in unique.items():
        pretty_unique[sorted(v)[-1]] = len(v)

    return pretty_unique 

def __init__(self, atom1, state1, atom2, state2):

        self.atom1 = atom1
        if (issubclass(type(self.atom1), DivalentAtom)
            and (len(state1) != 5 or (state1[4] != 0 and state1[4] != 1))
                ):
            raise ValueError("For divalent atoms state specification has to "
                             "include total spin angular momentum s as the last "
                             "number in the state specification [n,l,j,m_j,s].")
        self.state1 = state1
        # add exlicitly total spin of the state for Alkaline atoms
        if (len(self.state1) == 4): self.state1.append(0.5)

        self.atom2 = atom2
        if (issubclass(type(self.atom2), DivalentAtom)
                and (len(state1) != 5 or (state1[4] != 0 and state1[4] != 1))
                    ):
            raise ValueError("For divalent atoms state specification has to "
                             "include total spin angular momentum s as the last "
                             "numbre in the state specification [n,l,j,m_j,s].")
        self.state2 = state2
        # add exlicitly total spin of the state for Alkaline atoms
        if (len(self.state2) == 4):
            self.state2.append(0.5)

        self.pairStateEnergy = (
            atom1.getEnergy(self.state1[0],
                            self.state1[1],
                            self.state1[2],
                            s=self.state1[4])
            + atom2.getEnergy(self.state2[0],
                              self.state2[1],
                              self.state2[2],
                              s=self.state2[4])
            ) * C_e / C_h * 1e-9 

def getRecoilEnergy(self):
        """
        Recoil energy for atoms in given optical lattice

        Returns:
            float: recoil energy in units of J
        """
        latticeConstant = self.trapWavenegth / 2
        Er = C_h**2 / (8 * self.atom.mass * latticeConstant**2)
        return Er 

def getTransitionWavelength(self, n1, l1, j1, n2, l2, j2, s=0.5, s2=None):
        """
            Calculated transition wavelength (in vacuum) in m.

            Returned values is given relative to the centre of gravity of the
            hyperfine-split states.

            Args:
                n1 (int): principal quantum number of the state
                    **from** which we are going
                l1 (int): orbital angular momentum of the state
                    **from** which we are going
                j1 (float): total angular momentum of the state
                    **from** which we are going
                n2 (int): principal quantum number of the state
                    **to** which we are going
                l2 (int): orbital angular momentum of the state
                    **to** which we are going
                j2 (float): total angular momentum of the state
                    **to** which we are going
                s (float): optional, spin of the intial state
                    (for Alkali this is fixed to 0.5)
                s2 (float): optional, spin of the final state.
                    If not set, defaults to same value as :obj:`s`

            Returns:
                float:
                    transition wavelength (in m). If the returned value is
                    negative, level from which we are going is **above**
                    the level to which we are going.
        """
        if s2 is None:
            s2 = s
        return (C_h * C_c) / ((self.getEnergy(n2, l2, j2, s=s2)
                               - self.getEnergy(n1, l1, j1, s=s)) * C_e) 

def calcLuxEquivalent(wavelength,rad_min,rad_dynrange,units):
    """Calc the interpolation function between lux and photon rate radiance

    Assuming single wavelength colour, the specified wavelength value is used to 
    calculate the lux equivalent lux image for the radiance input range.

    Args:
        | wavelength (np.array): wavelength vector
        | sysresp (np.array): system response spectral vector
        | rad_min (float): minimum photon rate radiance lookup table 
        | rad_dynrange (float): maximum photon rate radiance in lookup table 
        | units (string): input radiance units q/s or W

    Returns:
        | interpolation function

    Raises:
        | No exception is raised.

    Author: CJ Willers
    """
    if 'q' in units:
        conversion = wavelength / (const.h * const.c)
    else:
        conversion = 1.

    Wm2tolux = 683 * 1.019 * np.exp(-285.51 * (wavelength*1e6 - 0.5591)**2)
    # convert from q/s to W
    rad_minW = rad_min / conversion
    rad_dynrangeW = rad_dynrange / conversion
    radW = np.linspace(0.99*rad_minW, 1.01*(rad_minW+rad_dynrangeW), 1000) 
    lux =  Wm2tolux * radW
    # convert from W back to q/s when setting up the function    
    fintp = interpolate.interp1d((radW*wavelength / (const.h * const.c)).reshape(-1), lux.reshape(-1))
    return fintp


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

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 deeStarPeak(wavelength,temperature,eta,halfApexAngle):
        i = 0
        for wlc in wavelength:
            wl =  np.linspace(wlc/100, wlc, 1000).reshape(-1, 1)
            LbackLambda = ryplanck.planck(wl,temperature, type='ql')/np.pi
            Lback = np.trapz(LbackLambda.reshape(-1, 1),wl, axis=0)[0]
            Eback = Lback * np.pi * (np.sin(halfApexAngle)) ** 2
            # funny construct is to prevent divide by zero
            tempvar = np.sqrt(eta/(Eback+(Eback==0))) * (Eback!=0) + 0 * (Eback==0)
            dstarwlc[i] = 1e-6 * wlc * tempvar/(const.h * const.c * np.sqrt(2))
            #print(Eback)
            i = i + 1
        return dstarwlc * 100. # to get cm units 

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 calc_abscoeff(tran, itera, alt=None):
    _up, _low = Tran_tag[tran][1:].astype(int)
    _v0 = Freqi[_up, _low]*1.e9
    _const = C.h*_v0/4./N.pi
    if alt is None:
        lowerPop = N.array(Ite_pop)[alt][itera][_low][0]
        upperPop = N.array(Ite_pop)[alt][itera][_up][0]
    else:
        lowerPop = N.array(Ite_pop)[:, itera, _low, 0]
        upperPop = N.array(Ite_pop)[:, itera, _up, 0]
    BLU = Blu[_up, _low]
    BUL = Bul[_up, _low]
    # AUL = Ai[_up, _low]
    return (lowerPop*BLU - upperPop*BUL)*_const 

def basic(LowerPop, UpperPop, Blu, Bul, Freq):
    u"""
    calculate absorption coefficient.\n
    Freq should be Hz
    """
    _const = h*Freq/4./np.pi
    return (LowerPop*Blu - UpperPop*Bul)*_const 

def PopuSource(LowerPop, UpperPop, LowerDeg, UpperDeg, Freq_c):
    sij = (2.*h*Freq_c**3)/(c**2*(LowerPop*UpperDeg/UpperPop/LowerDeg-1.))
    return sij 

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 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 noise_profile(self, df):
        """noise_profile(bw) computes amplifier ASE (W) in signal bandwidth (Hz)

        Noise is calculated at amplifier input

        :bw: signal bandwidth = baud rate in Hz
        :type bw: float

        :return: the asepower in W in the signal bandwidth bw for 96 channels
        :return type: numpy array of float

        ASE power using per channel gain profile inputs:

            NF_dB - Noise figure in dB, vector of length number of channels or
                    spectral slices
            G_dB  - Actual gain calculated for the EDFA, vector of length number of
                    channels or spectral slices
            ffs     - Center frequency grid of the channels or spectral slices in
                    THz, vector of length number of channels or spectral slices
            dF    - width of each channel or spectral slice in THz,
                    vector of length number of channels or spectral slices

        OUTPUT:

            ase_dBm - ase in dBm per channel or spectral slice

        NOTE:

            The output is the total ASE in the channel or spectral slice. For
            50GHz channels the ASE BW is effectively 0.4nm. To get to noise power
            in 0.1nm, subtract 6dB.

        ONSR is usually quoted as channel power divided by
        the ASE power in 0.1nm RBW, regardless of the width of the actual
        channel.  This is a historical convention from the days when optical
        signals were much smaller (155Mbps, 2.5Gbps, ... 10Gbps) than the
        resolution of the OSAs that were used to measure spectral power which
        were set to 0.1nm resolution for convenience.  Moving forward into
        flexible grid and high baud rate signals, it may be convenient to begin
        quoting power spectral density in the same BW for both signal and ASE,
        e.g. 12.5GHz."""

        ase = h * df * self.channel_freq * db2lin(self.nf)  # W
        return ase  # in W at amplifier input 

def getStateC3(self, n, l, j, coupledStatesList, s=0.5, debugOutput=False):
        r"""
        Van der Waals atom-surface interaction coefficient for
        a given state (:math:`C_3` in units of
        :math:`\mathrm{J}\cdot\mathrm{m}^3` )

        Args:
            n (int): principal quantum number of the state
            l (int): orbital angular momentum of the state
            j (int): total angular momentum of state
            coupledStatesList (array): array of states that are strongly
                dipole-coupled to the initial state, whose contribution
                to :math:`C_3` will be take into account. Format
                `[[n1,l1,j1],...]`
            s (float, optional): total spin angular momentum for the considered
                state. Default value of 0.5 is correct for `AlkaliAtoms`, but
                it has to be explicitly specifiied for `DivalentAtom`.
            debugOutput (bool, optional): prints additional output information,
                False by default.

        Returns:
            float, float:
                :math:`C_3` (in units of :math:`{\rm J}\cdot {\rm m}^3` ),
                estimated error :math:`\delta C_3`
        """

        if debugOutput:
            print("%s ->\tC3 contr. (kHz mum^3) \tlambda (mum)\tn"
                  % (printStateString(n, l, j, s=s))
                  )

        totalShift = 0
        sumSqError = 0

        for state in coupledStatesList:
            c3, err, refIndex = self.getC3contribution(
                n, l, j,
                state[0], state[1], state[2],
                s=s)
            if debugOutput:
                print(
                    "-> %s\t%.3f +- %.3f    \t%.3f\t\t%.3f\n" %
                    (printStateString(state[0], state[1], state[2], s=s),
                     c3/C_h*(1e6)**3*1e-3,
                     err/C_h*(1e6)**3*1e-3,
                     self.atom.getTransitionWavelength(
                         n, l, j,
                         state[0], state[1], state[2], s=s, s2=s) * 1e6,
                     refIndex)
                    )
            totalShift += c3
            sumSqError += err**2
        error = np.sqrt(sumSqError)
        if debugOutput:
            print("= = = = = = \tTotal shift of %s\t= %.3f+-%.4f kHz mum^3\n" %
                  (printStateString(n, l, j, s=s),
                   totalShift/C_h * (1e6)**3 * 1e-3,
                   error/C_h * (1e6)**3 * 1e-3))

        return totalShift, error  # in J m^3 

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 

def __isCoupled(self, n, l, j, nn, ll, jj, n1, l1, j1, n2, l2, j2, limit):
        if ((abs(self.__getEnergyDefect(n, l, j,
                                        nn, ll, jj,
                                        n1, l1, j1,
                                        n2, l2, j2)
                 ) / C_h < limit)
            and not (n == n1 and nn == n2
                     and l == l1 and ll == l2
                     and j == j1 and jj == j2)
                and not ((abs(l1 - l) != 1
                          and( (abs(j - 0.5) < 0.1
                                and abs(j1 - 0.5) < 0.1) # j = 1/2 and j'=1/2 forbidden
                                or
                                (abs(j) < 0.1
                                and abs(j1 - 1) < 0.1)  # j = 0 and j'=1 forbidden
                                or
                                (abs(j-1) < 0.1
                                and abs(j1) < 0.1)  # j = 1 and j'=0 forbidden
                                )
                          )
                         or (abs(l2 - ll) != 1
                            and( (abs(jj - 0.5) < 0.1
                                and abs(j2 - 0.5) < 0.1) # j = 1/2 and j'=1/2 forbidden
                                or
                                (abs(jj) < 0.1
                                and abs(j2 - 1) < 0.1)  # j = 0 and j'=1 forbidden
                                or
                                (abs(jj-1) < 0.1
                                and abs(j2) < 0.1)  # j = 1 and j'=0 forbidden
                                )
                          )
                        )
                and not(abs(j)<0.1 and abs(j1)<0.1)  # j = 0 and j'=0 forbiden
                and not (abs(jj)<0.1 and abs(j2)<0.1)
                and not (abs(l)<0.1 and abs(l1)<0.1) # l = 0 and l' = 0 is forbiden
                and not (abs(ll)<0.1 and abs(l2)<0.1)
                ):
            # determine coupling
            dl = abs(l - l1)
            dj = abs(j - j1)
            c1 = 0
            if dl == 1 and (dj < 1.1):
                c1 = 1  # dipole coupling
            elif (dl == 0 or dl == 2 or dl == 1)and (dj < 2.1) and \
                    (2 <= self.interactionsUpTo):
                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) and \
                    (2 <= self.interactionsUpTo):
                c2 = 2  # quadrupole coupling
            else:
                return False
            return c1 + c2
        else:
            return False 

def __init__(self, molecules_state_consts):
        self.filename = molecules_state_consts
        fn = open(molecules_state_consts + 'oH2O16.lev_7levels', 'rb')
        self.fr = fn.readlines()
        fn.close()
        self.ni = int(self.fr[3])  # NUMBER OF ENERGY LEVELS
        self.e_cm1 = np.array([float(self.fr[xx].split()[1])
                               for xx in range(5, 5 + self.ni)])
        self.weighti = np.array([float(self.fr[xx].split()[2])
                                 for xx in range(5, 5 + self.ni)])

        Ai = np.zeros((self.ni, self.ni))
        Freqi = np.zeros((self.ni, self.ni))
        E_Ki = np.zeros((self.ni, self.ni))
        B_place = np.ones((self.ni, self.ni))*-1.
        Tran_tag = []
        Ai_tag = []
        Nt = int(self.fr[3+self.ni+3])  # int(9)
        for xx in range(8+self.ni, 8+self.ni+Nt):  # Nt Number of transition ,24,twolevel ->16):
            """up: temp[1]-1, low: temp[2]-1"""
            temp = np.array(self.fr[xx].split()).astype(np.float64)
            Tran_tag.append(temp[:3]-1)  # for 1... -> 0....
            Ai[int(temp[1]-1), int(temp[2])-1] = temp[3]
            Ai_tag.append(temp[3])
            # Freqi[int(temp[1]-1), int(temp[2])-1] = temp[4]
            _deltaEcm1 = (self.e_cm1[int(temp[1])-1] -
                          self.e_cm1[int(temp[2])-1])
            Freqi[int(temp[1]-1), int(temp[2])-1] = (_deltaEcm1*c *
                                                     100.*1.e-9)  # GHz
            E_Ki[int(temp[1]-1), int(temp[2])-1] = temp[5]
            B_place[int(temp[1]-1), int(temp[2])-1] = int(temp[0]-1)  # induced
            B_place[int(temp[2])-1, int(temp[1]-1)] = int(temp[0]-1)  # abs.ed.
        Tran_tag = np.array(Tran_tag)*1
        Ai_tag = np.array(Ai_tag)*1.
        Nt = len(Tran_tag)  # int(9)
        self.ai = Ai
        self.freqi = Freqi
        self.tran_tag = Tran_tag
        self.b_place = B_place
        self.ai_tag = Ai_tag
        self.nt = Nt
        Freq_array = np.array([])
        for xx in range(Nt):
            _up = int(Tran_tag[xx][1])
            _low = int(Tran_tag[xx][2])
            Freq_array = np.hstack((Freq_array, Freqi[_up, _low]))
        self.freq_array = Freq_array
        Bul = self.ai*c**2/(2*h*(self.freqi*1.e9)**3)
        Blu = Bul*self.weighti.reshape((self.ni, 1))/self.weighti
        Bul[Bul != Bul] = 0  # [i,j]: i->j
        Blu[Blu != Blu] = 0  # [i,j]: j->i
        self.blu = Blu
        self.bul = Bul