Python scipy.constants.e() Examples

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 eVtoJoule(EeV):
    """
    Convert energy in eV to Joule.

    Args:
        | E: Energy in eV

    Returns:
        | EJ: Energy in J
    """

    return EeV * const.e


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

def JouleTeEv(EJ):
    """
    Convert energy in Joule to eV.

    Args:
        | EJ: Energy in J

    Returns:
        | EeV: Energy in eV
    """

    return EJ / const.e




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

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 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 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 _eFieldCouplingDivE(self, n1, l1, j1, mj1, n2, l2, j2, mj2, s=0.5):
        # eFied coupling devided with E (witout actuall multiplication to getE)
        # delta(mj1,mj2') delta(l1,l2+-1)
        if ((abs(mj1 - mj2) > 0.1) or (abs(l1 - l2) != 1)):
            return 0

        # matrix element
        result = self.atom.getRadialMatrixElement(n1, l1, j1,
                                                  n2, l2, j2,
                                                  s=s) *\
            physical_constants["Bohr radius"][0] * C_e

        sumPart = self.eFieldCouplingSaved.getAngular(l1, j1, mj1,
                                                      l2, j2, mj2,
                                                      s=s)
        return result * sumPart 

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 shotnoise(sensor_signal_in):
    r"""This routine adds photon shot noise to the signal of the photosensor that is in photons.

    The photon shot noise is due to the random arrival of photons and can be
    described by a Poisson process. Therefore, for each :math:`(i,j)`-th element of
    the matrix :math:`\Phi_{q}` that contains the number of collected photons, a photon
    shot noise  is simulated as a Poisson process :math:`\mathcal{P}` with mean
    :math:`\Lambda`:

    :math:`\Phi_{ph.shot}=\mathcal{P}(\Lambda), \,\,\,\,\mbox{ where   } \Lambda = \Phi_{q}.`

    We use the `ryutils.poissonarray` function that generates Poisson random numbers
    with mean :math:`\Lambda`.  That is, the number of collected photons in 
    :math:`(i,j)`-th pixel of the simulated photosensor in the matrix :math:`\Phi_{q}` is
    used as the mean :math:`\Lambda` for the generation of Poisson random numbers to
    simulate the photon shot noise. The input of the `ryutils.poissonarray` function will
    be the matrix :math:`\Phi_{q}` that contains the number of collected photons. The
    output will be the matrix :math:`\Phi_{ph.shot} \rightarrow \Phi_{q}`, i.e., the signal
    with added photon shot noise.  The matrix :math:`\Phi_{ph.shot}` is recalculated
    each time the simulations are started, which corresponds to the temporal nature
    of the photon shot noise.

    Args:
        | sensor_signal_in (np.array[N,M]): photon irradiance in, in photons
 
    Returns:
        | sensor_signal_out (np.array[N,M]): photon signal out, in photons

    Raises:
        | No exception is raised.

    Author: Mikhail V. Konnik, revised/ported by CJ Willers

    Original source: http://arxiv.org/pdf/1412.4031.pdf
    """

    # Since this is a shot noise, it must be random every time we apply it, unlike the Fixed Pattern Noise (PRNU or DSNU).
    # If seed is omitted or None, current system time is used
    return ryutils.poissonarray(sensor_signal_in, seedval=0)

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

def limitzero(a, thr=0.6):
    r"""Performs an asymetric clipping to prevent negative values.
    The lower-end values are clumped up towards the lower positive values, while 
    upper-end values are not affected.  

    This function is used to prevent negative random variables for wide sigma and low
    mean value, e.g., N(1,.5).  If the random variables are passed through this function
    The resulting distribution is not normal any more, and has no known analytical form.

    A threshold value of around 0.6 was found to work well for N(1,small) up to N(1,.5).

    Before you use this function, first check the results using the code below in the main
    body of this file.

    Args:
        | a (np.array): an array of floats, 

    Returns:
        | scaling factors.

    Raises:
        | No exception is raised.

    Author: CJ Willers
    """
    ashape = a.shape
    a = a.flatten()
    a = np.where(a<thr, thr * np.exp(( a - thr) / (1 * thr)), 0) + np.where(a >= thr, a, 0)

    return   a.reshape(ashape)


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

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 getTransitionFrequency(self, n1, l1, j1, n2, l2, j2, s=0.5, s2=None):
        """
            Calculated transition frequency in Hz

            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 the same value as :obj:`s`

            Returns:
                float:
                    transition frequency (in Hz). 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 (self.getEnergy(n2, l2, j2, s=s2)
                - self.getEnergy(n1, l1, j1, s=s)) * C_e / C_h 

def getRabiFrequency2(self,
                          n1, l1, j1, mj1,
                          n2, l2, j2, q,
                          electricFieldAmplitude,
                          s=0.5):
        """
            Returns a Rabi frequency for resonant excitation with a given
            electric field amplitude

            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)
                electricFieldAmplitude : amplitude of electric field
                    driving (V/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`
        """
        mj2 = mj1 + q
        if abs(mj2) - 0.1 > j2:
            return 0
        dipole = self.getDipoleMatrixElement(n1, l1, j1, mj1,
                                             n2, l2, j2, mj2, q,
                                             s=s) *\
            C_e * physical_constants["Bohr radius"][0]
        freq = electricFieldAmplitude * abs(dipole) / hbar
        return freq 

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 getEnergyDefect(self, n, l, j, n1, l1, j1, n2, l2, j2,
                        s=0.5):
        """
            Energy defect for the given two pair-states (one of the state has
            two atoms in the same state)

            Energy difference between the states
            :math:`E(n_1,l_1,j_1,n_2,l_2,j_2) - E(n,l,j,n,l,j)`

            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. Spin angular momentum
                    (default 0.5 for Alkali)

            Returns:
                float:  energy defect (SI units: J)
        """
        return C_e * (self.getEnergy(n1, l1, j1, s=s)
                      + self.getEnergy(n2, l2, j2, s=s)
                      - 2 * self.getEnergy(n, l, j, s=s)) 

def __getEnergyDefect(self,
                          n, l, j,
                          nn, ll, jj,
                          n1, l1, j1,
                          n2, l2, j2):
        """
        Energy defect between |n,l,j>x|nn,ll,jj> state and |n1,l1,j1>x|n1,l1,j1>
        state of atom1 and atom2 in respective spins states s1 and s2

        Takes spin vales s1 and s2 as the one defined when defining calculation.

        Args:
            n (int): principal quantum number
            l (int): orbital angular momenutum
            j (float): total angular momentum
            nn (int): principal quantum number
            ll (int): orbital angular momenutum
            jj (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

        Returns:
            float:  energy defect (SI units: J)
        """
        return C_e * (self.atom1.getEnergy(n1, l1, j1, s=self.s1)
                      + self.atom2.getEnergy(n2, l2, j2, s=self.s2)
                      - self.atom1.getEnergy(n, l, j, s=self.s1)
                      - self.atom2.getEnergy(nn, ll, jj, s=self.s2)) 

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 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 convert_to_electrons(strh5):
    r"""This routine converts photon rate irradiance to electron rate irradiance
     
    Args:
        | strh5 (hdf5 file): hdf5 file that defines all simulation parameters
 
    Returns:
        | in strh5: (hdf5 file) updated data fields

    Raises:
        | No exception is raised.

    Author: Mikhail V. Konnik, revised/ported by CJ Willers

    Original source: http://arxiv.org/pdf/1412.4031.pdf
    """

    # Converting the signal from Photons to Electrons
    # Quantum Efficiency = Quantum Efficiency Interaction X Quantum Yield Gain.
    # diagram node 4 quantum efficiency stored in rystare/quantumEfficiency
    strh5['rystare/quantumEfficiency'][...] = strh5['rystare/photondetector/externalquantumeff'][()] * strh5['rystare/photondetector/quantumyield'][()]

    # number of electrons [e] generated in detector
    strh5['rystare/signal/electronRateIrradiance'][...] = strh5['rystare/signal/photonRateIrradianceNU'][()] * strh5['rystare/quantumEfficiency'][()] 
    # diagram node 4 photon rate x mean value of the quantum efficiency stored in rystare/signal/electronRateIrradiance

    return strh5

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

def multiply_integration_time(strh5):
    r"""This routine multiplies with integration time
     
    The input to the model of the photosensor is assumed to be a matrix :math:`E_{q}\in R^{N\times M}` 
    that has been converted to electrons, corresponding to electron rate  [e/s].  
    The electron rate is converted to electron count into the pixel by accounting for 
    detector integration time:

    :math:`\Phi_{q}  =  \textrm{round} \left(  E_{q} \cdot t_I  \right),`

    where :math:`t_{I}` is integration     (exposure) time.

     Args:
        | strh5 (hdf5 file): hdf5 file that defines all simulation parameters
 
    Returns:
         | in strh5: (hdf5 file) updated data fields

    Raises:
        | No exception is raised.

    Author: Mikhail V. Konnik, revised/ported by CJ Willers

    Original source: http://arxiv.org/pdf/1412.4031.pdf
    """

    #the number of electrons accumulated during the integration time  (rounded).
    strh5['rystare/signal/lightelectronsnoshotnoise'][...] = np.round(strh5['rystare/signal/electronRate'][()] * strh5['rystare/photondetector/integrationtime'][()])

    return strh5

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

def multiply_detector_area(strh5):
    r"""This routine multiplies  detector area

    The input to the model of the photosensor is assumed to be a matrix :math:`E_{q}\in R^{N\times M}` 
    that has been converted to electronrate irradiance, corresponding to electron rate  [e/(m2.s)].  
    The electron rate irriance is converted to electron rate into the pixel by accounting for 
    detector area:

    :math:`\Phi_{q}  =  \textrm{round} \left(  E_{q} \cdot P_A    \right),`

    where :math:`P_A` is the area of a pixel [m2].
     
     Args:
        | strh5 (hdf5 file): hdf5 file that defines all simulation parameters
 
    Returns:
         | in strh5: (hdf5 file) updated data fields

    Raises:
        | No exception is raised.

    Author: Mikhail V. Konnik, revised/ported by CJ Willers

    Original source: http://arxiv.org/pdf/1412.4031.pdf
    """


    # calculate radiant flux [W] from irradiance [W/m^2] and area
    strh5['rystare/signal/electronRate'][...]  = strh5['rystare/detectorArea'][()] * strh5['rystare/signal/electronRateIrradiance'][()] 

    return strh5

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

def dplanck(spectral, temperature, type='el'):
    """Temperature derivative of Planck law exitance.

    Calculates the temperature derivative for Planck law spectral exitance
    from a surface at the stated temperature. dM/dT can be given in radiant or
    photon rate units, depending on user input in type. Temperature can be a 
    scalar, a list or an array. 

    Args:
        | spectral (scalar, np.array (N,) or (N,1)):  spectral vector in  [micrometer], [cm-1] or [Hz].
        | temperature (scalar, list[M], np.array (M,), (M,1) or (1,M)):  Temperature in [K]
        | type (string):
        |  'e' signifies Radiant values in [W/(m^2.K)].
        |  'q' signifies photon rate values  [quanta/(s.m^2.K)].
        |  'l' signifies wavelength spectral vector  [micrometer].
        |  'n' signifies wavenumber spectral vector [cm-1].
        |  'f' signifies frequency spectral vecor [Hz].

    Returns:
        | (scalar, np.array[N,M]):  spectral radiant exitance (not radiance) in units selected.
        | For type = 'el' units will be [W/(m2.um.K)]
        | For type = 'qf' units will be [q/(s.m2.Hz.K)]
        | Other return types are similarly defined as above.
        | Returns None on error.

    Raises:
        | No exception is raised, returns None on error.
    """

    if type in list(dplancktype.keys()):
        #select the appropriate fn as requested by user
        exitance = dplancktype[type](spectral, temperature)
    else:
        # return all zeros if illegal type
        exitance = - np.ones(spectral.shape)

    return exitance


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

def planck(spectral, temperature, type='el'):
    """Planck law spectral exitance.

    Calculates the Planck law spectral exitance from a surface at the stated 
    temperature. Temperature can be a scalar, a list or an array. Exitance can 
    be given in radiant or photon rate units, depending on user input in type.

    Args:
        | spectral (scalar, np.array (N,) or (N,1)):  spectral vector.
        | temperature (scalar, list[M], np.array (M,), (M,1) or (1,M)):  Temperature in [K]
        | type (string):
        |  'e' signifies Radiant values in [W/m^2.*].
        |  'q' signifies photon rate values  [quanta/(s.m^2.*)].
        |  'l' signifies wavelength spectral vector  [micrometer].
        |  'n' signifies wavenumber spectral vector [cm-1].
        |  'f' signifies frequency spectral vecor [Hz].

    Returns:
        | (scalar, np.array[N,M]):  spectral radiant exitance (not radiance) in units selected.
        | For type = 'el' units will be [W/(m^2.um)].
        | For type = 'qf' units will be [q/(s.m^2.Hz)].
        | Other return types are similarly defined as above.
        | Returns None on error.

    Raises:
        | No exception is raised, returns None on error.
    """
    if type in list(plancktype.keys()):
        #select the appropriate fn as requested by user
        exitance = plancktype[type](spectral, temperature)
    else:
        # return all minus one if illegal type
        exitance = None

    return exitance


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

def stefanboltzman(temperature, type='e'):
    """Stefan-Boltzman wideband integrated exitance.

    Calculates the total Planck law exitance, integrated over all wavelengths,
    from a surface at the stated temperature. Exitance can be given in radiant or
    photon rate units, depending on user input in type.

    Args:
        | (scalar, list[M], np.array (M,), (M,1) or (1,M)):  Temperature in [K]
        | type (string):  'e' for radiant or 'q' for photon rate exitance.

    Returns:
        | (float): integrated radiant exitance in  [W/m^2] or [q/(s.m^2)].
        | Returns a -1 if the type is not 'e' or 'q'

    Raises:
        | No exception is raised.
    """

    #confirm that only vector is used, break with warning if so.
    if isinstance(temperature, np.ndarray):
        if len(temperature.flat) != max(temperature.shape):
            print('ryplanck.stefanboltzman: temperature must be of shape (M,), (M,1) or (1,M)')
            return -1

    tempr = np.asarray(temperature).astype(float)
    #use dictionary to switch between options, lambda fn to calculate, default -1
    rtnval = {
              'e': lambda temp: pconst.sigmae * np.power(temp, 4) ,
              'q': lambda temp: pconst.sigmaq * np.power(temp, 3)
              }.get(type, lambda temp: -1)(tempr)
    return rtnval



################################################################
##
# dictionaries to allow case-like statements in generic planck functions, below. 

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 __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 responsivity_FPN_light(strh5):
    r"""Multiploying the photon signal with the PRNU.

    The Photo Response Non-Uniformity (PRNU) is the spatial variation in pixel
    conversion gain (from photons to electrons). When viewing a uniform scene the pixel signals will 
    differ because of the PRNU,  mainly due to variations in the individual pixel's characteristics
    such as detector area and spectral response. These variations occur during the manufacture 
    of the substrate and the detector device.
 
    The PRNU is signal-dependent (proportional to the input signal) and is
    fixed-pattern (time-invariant). For visual (silicon) sensors the PRNU factor is typically
    :math:`0.01\dots 0.05`\ , but for HgCdTe sensors it can be as large as  :math:`0.02\dots 0.25`\ .
    It varies from sensor to sensor, even within the same manufacturing batch. 

    The photo response non-uniformity (PRNU) is considered as a temporally-fixed light 
    signal non-uniformity. The PRNU is modelled using a Gaussian distribution for each
    :math:`(i,j)`-th pixel of the matrix :math:`I_{e^-}`, 
    as :math:`I_{PRNU.e^-}=I_{e^-}(1+\mathcal{N}(0,\sigma_{PRNU}^2))`
    where :math:`\sigma_{PRNU}` is the PRNU factor value. 

    Args:
        | strh5 (hdf5 file): hdf5 file that defines all simulation parameters
 
    Returns:
        | in strh5: (hdf5 file) updated data fields

    Raises:
        | No exception is raised.

    Author: Mikhail V. Konnik, revised/ported by CJ Willers

    Original source: http://arxiv.org/pdf/1412.4031.pdf
    """

    #the random generator seed is fixed with value from input
    np.random.seed(strh5['rystare/photondetector/lightPRNU/seed'][()])

    #matrix for the PRNU
    normalisedVariation = FPN_models(
        strh5['rystare/imageSizePixels'][()][0], strh5['rystare/imageSizePixels'][()][1],
        'pixel', strh5['rystare/photondetector/lightPRNU/model'][()], strh5['rystare/photondetector/lightPRNU/sigma'][()])

    # diagram node 2  NU stored in 'rystare/photondetector/lightPRNU/value'

    #np.random.randn has mean=0, variance = 1, so we multiply with variance and add to mean
    strh5['rystare/photondetector/lightPRNU/value'][...] = (1 + normalisedVariation)

    # diagram node 3 photon rate multiplied with PRNU stored in 'rystare/signal/photonRateIrradianceNU'
    #apply the PRNU noise to the light signal of the photosensor.
    strh5['rystare/signal/photonRateIrradianceNU'][...] = strh5['rystare/signal/photonRateIrradiance'][()] * strh5['rystare/photondetector/lightPRNU/value'][()]
        
    return strh5


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

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) ================== 

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 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 )


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