qutip/utilities.py
"""
This module contains utility functions that are commonly needed in other
qutip modules.
"""
__all__ = ['n_thermal', 'clebsch', 'convert_unit']
import numpy as np
def n_thermal(w, w_th):
"""
Return the number of photons in thermal equilibrium for an harmonic
oscillator mode with frequency 'w', at the temperature described by
'w_th' where :math:`\\omega_{\\rm th} = k_BT/\\hbar`.
Parameters
----------
w : float or ndarray
Frequency of the oscillator.
w_th : float
The temperature in units of frequency (or the same units as `w`).
Returns
-------
n_avg : float or array
Return the number of average photons in thermal equilibrium for a
an oscillator with the given frequency and temperature.
"""
if isinstance(w, np.ndarray):
return 1.0 / (np.exp(w / w_th) - 1.0)
else:
if (w_th > 0) and np.exp(w / w_th) != 1.0:
return 1.0 / (np.exp(w / w_th) - 1.0)
else:
return 0.0
def _factorial_prod(N, arr):
arr[:int(N)] += 1
def _factorial_div(N, arr):
arr[:int(N)] -= 1
def _to_long(arr):
prod = 1
for i, v in enumerate(arr):
prod *= (i+1)**int(v)
return prod
def clebsch(j1, j2, j3, m1, m2, m3):
"""Calculates the Clebsch-Gordon coefficient
for coupling (j1,m1) and (j2,m2) to give (j3,m3).
Parameters
----------
j1 : float
Total angular momentum 1.
j2 : float
Total angular momentum 2.
j3 : float
Total angular momentum 3.
m1 : float
z-component of angular momentum 1.
m2 : float
z-component of angular momentum 2.
m3 : float
z-component of angular momentum 3.
Returns
-------
cg_coeff : float
Requested Clebsch-Gordan coefficient.
"""
if m3 != m1 + m2:
return 0
vmin = int(np.max([-j1 + j2 + m3, -j1 + m1, 0]))
vmax = int(np.min([j2 + j3 + m1, j3 - j1 + j2, j3 + m3]))
c_factor = np.zeros((int(j1 + j2 + j3 + 1)), np.int32)
_factorial_prod(j3 + j1 - j2, c_factor)
_factorial_prod(j3 - j1 + j2, c_factor)
_factorial_prod(j1 + j2 - j3, c_factor)
_factorial_prod(j3 + m3, c_factor)
_factorial_prod(j3 - m3, c_factor)
_factorial_div(j1 + j2 + j3 + 1, c_factor)
_factorial_div(j1 - m1, c_factor)
_factorial_div(j1 + m1, c_factor)
_factorial_div(j2 - m2, c_factor)
_factorial_div(j2 + m2, c_factor)
C = np.sqrt((2.0 * j3 + 1.0)*_to_long(c_factor))
s_factors = np.zeros(((vmax + 1 - vmin), (int(j1 + j2 + j3))), np.int32)
sign = (-1) ** (vmin + j2 + m2)
for i,v in enumerate(range(vmin, vmax + 1)):
factor = s_factors[i,:]
_factorial_prod(j2 + j3 + m1 - v, factor)
_factorial_prod(j1 - m1 + v, factor)
_factorial_div(j3 - j1 + j2 - v, factor)
_factorial_div(j3 + m3 - v, factor)
_factorial_div(v + j1 - j2 - m3, factor)
_factorial_div(v, factor)
common_denominator = -np.min(s_factors, axis=0)
numerators = s_factors + common_denominator
S = sum([(-1)**i * _to_long(vec) for i,vec in enumerate(numerators)]) * \
sign / _to_long(common_denominator)
return C * S
# -----------------------------------------------------------------------------
# Functions for unit conversions
#
_e = 1.602176565e-19 # C
_kB = 1.3806488e-23 # J/K
_h = 6.62606957e-34 # Js
_unit_factor_tbl = {
# "unit": "factor that convert argument from unit 'unit' to Joule"
"J": 1.0,
"eV": _e,
"meV": 1.0e-3 * _e,
"GHz": 1.0e9 * _h,
"mK": 1.0e-3 * _kB,
}
def convert_unit(value, orig="meV", to="GHz"):
"""
Convert an energy from unit `orig` to unit `to`.
Parameters
----------
value : float / array
The energy in the old unit.
orig : str, {"J", "eV", "meV", "GHz", "mK"}, default: "meV"
The name of the original unit.
to : str, {"J", "eV", "meV", "GHz", "mK"}, default: "GHz"
The name of the new unit.
Returns
-------
value_new_unit : float / array
The energy in the new unit.
"""
if orig not in _unit_factor_tbl:
raise TypeError("Unsupported unit %s" % orig)
if to not in _unit_factor_tbl:
raise TypeError("Unsupported unit %s" % to)
return value * (_unit_factor_tbl[orig] / _unit_factor_tbl[to])
def convert_GHz_to_meV(w):
"""
Convert an energy from unit GHz to unit meV.
Parameters
----------
w : float / array
The energy in the old unit.
Returns
-------
w_new_unit : float / array
The energy in the new unit.
"""
# 1 GHz = 4.1357e-6 eV = 4.1357e-3 meV
w_meV = w * 4.1357e-3
return w_meV
def convert_meV_to_GHz(w):
"""
Convert an energy from unit meV to unit GHz.
Parameters
----------
w : float / array
The energy in the old unit.
Returns
-------
w_new_unit : float / array
The energy in the new unit.
"""
# 1 meV = 1.0/4.1357e-3 GHz
w_GHz = w / 4.1357e-3
return w_GHz
def convert_J_to_meV(w):
"""
Convert an energy from unit J to unit meV.
Parameters
----------
w : float / array
The energy in the old unit.
Returns
-------
w_new_unit : float / array
The energy in the new unit.
"""
# 1 eV = 1.602e-19 J
w_meV = 1000.0 * w / _e
return w_meV
def convert_meV_to_J(w):
"""
Convert an energy from unit meV to unit J.
Parameters
----------
w : float / array
The energy in the old unit.
Returns
-------
w_new_unit : float / array
The energy in the new unit.
"""
# 1 eV = 1.602e-19 J
w_J = 0.001 * w * _e
return w_J
def convert_meV_to_mK(w):
"""
Convert an energy from unit meV to unit mK.
Parameters
----------
w : float / array
The energy in the old unit.
Returns
-------
w_new_unit : float / array
The energy in the new unit.
"""
# 1 mK = 0.0000861740 meV
w_mK = w / 0.0000861740
return w_mK
def convert_mK_to_meV(w):
"""
Convert an energy from unit mK to unit meV.
Parameters
----------
w : float / array
The energy in the old unit.
Returns
-------
w_new_unit : float / array
The energy in the new unit.
"""
# 1 mK = 0.0000861740 meV
w_meV = w * 0.0000861740
return w_meV
def convert_GHz_to_mK(w):
"""
Convert an energy from unit GHz to unit mK.
Parameters
----------
w : float / array
The energy in the old unit.
Returns
-------
w_new_unit : float / array
The energy in the new unit.
"""
# h v [Hz] = kB T [K]
# h 1e9 v [GHz] = kB 1e-3 T [mK]
# T [mK] = 1e12 * (h/kB) * v [GHz]
w_mK = w * 1.0e12 * (_h / _kB)
return w_mK
def convert_mK_to_GHz(w):
"""
Convert an energy from unit mK to unit GHz.
Parameters
----------
w : float / array
The energy in the old unit.
Returns
-------
w_new_unit : float / array
The energy in the new unit.
"""
w_GHz = w * 1.0e-12 * (_kB / _h)
return w_GHz
def _version2int(version_string):
str_list = version_string.split(
"-dev")[0].split("rc")[0].split("a")[0].split("b")[0].split(
"post")[0].split('.')
return sum([int(d if len(d) > 0 else 0) * (100 ** (3 - n))
for n, d in enumerate(str_list[:3])])