''' Unit handling using external "quantities" module.
.. codeauthor:: Yoshifumi Futaana
.. note::
This module uses quantities package.
I appreciate the author, Darren Dale, for providing the software.
https://github.com/python-quantities/python-quantities
Unitq handles the unit using quantities package.
See http://pythonhosted.org/quantities/ for detailed documentation.
Additional functionalies are
* Shorter names of handling the qunatities. (For example, use :meth:`r` instead of :meth:`rescale`)
* Simple conversion of the string formatted unit to the unit object.
(Use the function :func:`nit`, or :func:`unit`)
* Some new units commonly used in plasma physics. (For example, :attr:`cu_flux`.)
* Dynamic definition of un-defined unit using SI prefix. (For example, 'pT' is automatically
defined as pico-tesla if you use it in :func:`nit` function.)
* You can set default unit system. (:func:`set_default_units`,
or wrapped function :func:`use` or specific functions
:func:`use_mksa` and :func:`use_plasma`.)
* Simple pretty-print (:func:`pprint`, or simply :func:`pp`)
* Helper function for implementing formula, simple extention for unitful
formula is possible. See :func:`make_unitful_function`. (Mainly for developer. It is
used in :mod:`irfpy.util.physics` intensively.
*Usage*
Basically, you may start with
>>> import irfpy.util.unitq as u
to use the unit system.
Simple example to create a unitful quantity is:
>>> print(8.57 * u.m)
8.57 m
>>> u.pprint(8.57 * u.m / (2.48 * u.s), fmt='%.3f')
3.456 m/s
Or you can make an unitful object using :func:`nit` function.
>>> u.pprint(u.nit('300 MPa')) # 300 Mega pascal
300.000 MPa
>>> u.pprint(u.nit('3.6 km/h').simplified)
1.000 m/s
You may get the value only as ``np.array`` object as a matter of fact.
>>> val = [8.57, 2.55] * u.m
>>> val_arr = val.n # .n will return the numpy array
>>> print(val_arr, isinstance(val_arr, np.ndarray)) # doctest: +NORMALIZE_WHITESPACE
[8.57 2.55] True
>>> val_arr[0] = 1.55
>>> print(val) # The change in the new array will not affect the original # doctest: +NORMALIZE_WHITESPACE
[8.57 2.55] m
Sample with angle and solid angle
>>> print(3.5 * u.sr)
3.5 sr
>>> print((1.5 * u.sr)**2)
2.25 sr**2
>>> print(((0.3 * u.sr)**2).s) # .s is to simplify. sr is unitless.
0.09 dimensionless
Electron volt temperature would be helpful.
>>> ev = 1 * eV_asT # A new unit eV_asT has been defined.
>>> print(ev.r('K')) # .r is to rescale the unit
11604.5052897 K
*Interactive use (``n``, ``u``, ``s``, ``r``, ``ns``, ``nr``, ``us``, ``ur``)*
The module is useful also for the interactive calculation.
Aliases are prepared for easy access to some methods/properties.
* ``n`` (numpy) returns the ``numpy.array`` expression of the quantity (``5.0`` or ``[2.0, 3.5]``).
* ``u`` (unit) returns the string expression of the units (``kg`` or ``cm/s``).
* ``s`` (simplified) is an alias to :meth:`simplified`.
* ``r`` (rescale) is an alias to :meth:`rescale`.
* ``ns`` (numpy simplified) returns the ``numpy.array`` after being simplified.
* ``nr`` (numpy rescaled) returns the ``numpy.array`` after being rescaled.
* ``us`` (unit simplified) returns the string expression of the unit after being simplified.
* ``ur`` (unit rescaled) returns the string expression of the unit after begin rescaled.
A series of sample is as follows.
Prepare the Planck constant. Internal expression is using
the unit of "Planck constant".
>>> h = k.Planck_constant # Planck constant
>>> print(h)
1 h (Planck_constant)
>>> print(h.n) # numpy.array of internal expression.
1.0
>>> print(h.u) # the unit of internal expression.
h
You may want to get the quantity in MKSA.
Usually you use ``simplified`` property to get it.
You may use ``s`` for a short expression.
To separate to the value and the unit, ``ns`` and ``us`` is prepared.
>>> print(h.simplified) # Usually, you use simplified for getting MKSA expression.
6.62606896e-34 kg*m**2/s
>>> print(h.s) # An alias is prepared.
6.62606896e-34 kg*m**2/s
>>> print(h.ns) # This simplifies the quantity, return as numpy array.
6.62606896e-34
>>> print(h.us) # This returns the unit of simplified quantity.
kg*m**2/s
Rescaling to another unit system is also done via ``r``, ``nr`` and ``ur``.
>>> pprint(h.r('g*cm^2/s'), fmt='%.3e') # You may rescale with ``r``.
6.626e-27 g*cm**2/s
>>> print('%.3e' % h.nr('g*cm^2/s')) # Return the rescaled value (float) only.
6.626e-27
>>> print(h.ur('g*cm^2/s')) # Return the rescaled unit
g*cm**2/s
'''
import logging as _logging
import re as _re
import quantities as _pq
import quantities.constants as k
from quantities import *
_pq.Quantity.s = _pq.Quantity.simplified
''' Alias to simplified.
'''
_pq.Quantity.n = lambda self: np.array(self)
_pq.Quantity.n = property(_pq.Quantity.n)
_pq.Quantity.u = lambda self: self.dimensionality.string
_pq.Quantity.u = property(_pq.Quantity.u)
_pq.Quantity.r = _pq.Quantity.rescale
''' Alias to rescale.
'''
_pq.Quantity.ns = lambda self: np.array(self.simplified)
_pq.Quantity.ns = property(_pq.Quantity.ns)
''' Get numpy array with simplied unit.
'''
_pq.Quantity.nr = lambda self, unit: np.array(self.rescale(unit))
''' Get numpy array after the rescale
'''
_pq.Quantity.ur = lambda self, unit: self.rescale(unit).dimensionality.string
_pq.Quantity.us = lambda self: self.simplified.dimensionality.string
_pq.Quantity.us = property(_pq.Quantity.us)
''' Return the string expression of the simpilied unit.
'''
[docs]class k2:
''' MKSA scalar constant.
.. deprecated:: 4.2.5
Use :mod:`irfpy.util.constant` module.
>>> print('%.3e' % k2.amu)
1.661e-27
'''
amu = k.amu.nr(kg)
''' AMU'''
e = k.e.nr(C)
''' Elementary charge'''
m_p = k.m_p.nr(kg)
''' Proton mass'''
CUnit = CompoundUnit
''' Alias to ``CompoundUnit``
'''
cu_pflux = CUnit('1 /cm**2 /s')
''' Unit of particle flux
'''
cu_flux = CUnit('1 /cm**2 /s')
''' Unit of particle flux
'''
cu_eflux = CUnit('eV/cm**2 /s')
''' Unit of energy flux
'''
cu_diffflux = CUnit('1/cm**2 /sr /eV / s')
''' Unit of differential particle flux
'''
cu_ediffflux = CUnit('eV/cm**2 /sr /eV/s')
''' Unit of differential energy flux
'''
cu_gfactor = CUnit('cm**2 * sr * eV/eV')
''' Unit of g-factor
'''
eVT = eV_asT = UnitQuantity('electron_volt_temperature',
K * 11604.5052897, symbol='eVT')
''' Electron volt used as temperature
'''
eV_asE = eV
''' Electron volt used as energy (quantities defined object)
'''
ster = sr
''' Steredian, an alias
'''
qe = q = e
''' Unit charge, an alias
'''
rydberg_energy = UnitQuantity('rydberg_energy', _pq.hartree_energy / 2., 'E_Ry')
""" Rydberg energy. It is about 13.6 eV, and half of the Hartree.
"""
si_prefix = {
'Y': 1e24, 'yotta': 1e24,
'Z': 1e21, 'zetta': 1e21,
'E': 1e18, 'exa': 1e18,
'P': 1e15, 'peta': 1e15,
'T': 1e12, 'tera': 1e12,
'G': 1e9, 'giga': 1e9,
'M': 1e6, 'mega': 1e6,
'k': 1e3, 'kilo': 1e3,
'h': 1e2, 'hecto': 1e2,
'da': 1e1, 'deca': 1e1, 'D':1e1,
'd': 0.1, 'deci': 1e-1,
'c': 1e-2, 'centi': 1e-2,
'm': 1e-3, 'milli': 1e-3,
'u': 1e-6, 'micro': 1e-6,
'n': 1e-9, 'nano': 1e-9,
'p': 1e-12, 'pico': 1e-12,
'f': 1e-15, 'femto': 1e-15,
'a': 1e-18, 'atto': 1e-18,
'z': 1e-21, 'zepto': 1e-21,
'y': 1e-24, 'yocto': 1e-24,
}
''' SI prefix'''
known_unit = {}
'''A dictionary of unit.
It maps the name of the unit and the UnitQuantity object.
'''
for var in list(globals().keys()):
if isinstance(globals()[var], _pq.UnitQuantity):
known_unit[var] = globals()[var]
def __token2unit(token, verbose=False):
r''' Helper function for nit
Convert token to unit.
If token is a number (float), it is returned.
Otherwise, token is thought to be the unit.
First, token's ^ is omitted, and then,
Token should be regexp as:
match = re.match('^([A-z]+)?(-?[0-9]+(\.\d+)?)?$', token2)
If the unit part cannot be found, a SI prefixed unit
is searched.
'''
logger = _logging.getLogger(__name__ + '.token2unit')
if verbose:
logger.setLevel(_logging.DEBUG)
else:
logger.setLevel(_logging.WARN)
try:
### If the token is just a number, simple.
return float(token)
except ValueError as e:
pass
token2 = token.replace('^', '') # Remove '^', then alphabet only, or alphabet followed by number
match = _re.match(r'^([A-z]+)?(-?[0-9]+(\.\d+)?)?$', token2)
if match is None:
raise ValueError('Unknown format, %s' % token)
un, mul, mulin = match.groups()
logger.debug('Match: The alphabet part: %s' % str(un))
logger.debug('Match: The number part: %s' % str(mul))
logger.debug('Match: The fraction part: %s' % str(mulin))
if mul is None:
mul = 1 # No number means the ^1
else:
if mulin is None: # Fraction part is missing then int
mul = int(mul)
else:
mul = float(mul)
try:
### If the system knows the name (alphabet part),
### the unit multiplied by the number is returned.
unitobject = known_unit[un] ** mul
except (AttributeError, TypeError, KeyError) as e:
### Try finding the prefix.
logger.debug('No conversion found. Try to find prefix')
for pref in si_prefix:
# logger.debug('Trying to find prefix "%s"' % pref)
if un.startswith(pref): # Ok, this is not very efficient way, but
logger.debug(' -- Prefix "%s" found' % pref)
unpref = un[:len(pref)] # Unit prefix
unun = un[len(pref):] # Main unit
logger.debug(' -- Prefix "%s" Unit "%s"' % (unpref, unun))
try:
unitobject = known_unit[unun] ** mul # Try to find the main unit
unitobject = unitobject * si_prefix[unpref] # Multiply by the prefix
try:
newunit = _pq.UnitQuantity(un, known_unit[unun] * si_prefix[unpref], symbol=un)
globals()[un] = newunit
known_unit[un] = newunit
except KeyError:
pass
return unitobject
except (AttributeError, TypeError, KeyError, IndexError) as e:
logger.debug(e)
pass
raise ValueError('Unknown unit, %s' % token)
return unitobject
[docs]def nit(string_units, verbose=False):
''' Return the unit object corresponding to the given expression
>>> import irfpy.util.unitq as u
>>> print(3 * u.nit('kg'))
3.0 kg
As you see the above example, the name comes for above convenstion :).
Each unit should be one word: ``'kg'`` or ``'s'``.
For the multiple of units, separate by space or aster or dot.:
e.g. ``'J s'`` or ``'J*s'`` or ``'J.s'``. ``'Js'`` not allowed.
Devision is ``'/'``. ``'1/s'``, ``m/s``.
For the power, you may use ``**`` or ``^``. No space before and after the operator.
e.g. ``'m2'``, ``'m^2'`` or ``'m**2``'
Anyway, it is highly recommended to make ``*`` or space between units
and numbers to multiply. Also, '^' is recommended to be added, but not spaces around.
(``'m2'`` and ``'m 2'`` provide different result. ``'m^2'`` and ``'m*2'``)
Parenthesis ``()`` not supported.
>>> print(u.nit('m'))
1 m (meter)
>>> print(u.nit(' 3 m'))
3.0 m
>>> print(u.nit('m/s'))
1.0 m/s
>>> print(u.nit('J s'))
1.0 s*J
>>> print(u.nit('J2 s'))
1.0 s*J**2
>>> print(u.nit('J2 s-2'))
1.0 J**2/s**2
>>> print(u.nit('J3 s-1 .8'))
0.8 J**3/s
>>> print(u.nit('3 / 2')) # Even usual number calculation.
1.5
>>> # print u.nit('3**2') ## But this is not possible
>>> print(u.nit('s0.5')) # Fractional unit ok.
1.0 s**0.5
>>> #print u.nit('3.s') # Will be error
>>> #print u.nit('.3s') # Will be error
>>> print(u.nit('s.K 3.'))
3.0 s*K
>>> #print u.nit('s.K3.') # Will be error
>>> print(u.nit('m2'))
1.0 m**2
>>> print(u.nit('m 2'))
2.0 m
>>> print(u.nit('m^2'))
1.0 m**2
>>> print(u.nit('m*2'))
2.0 m
It will guess the SI prefix.
>>> print(u.nit('pC')) # Pico coulon
1e-12 C
>>> pprint(u.nit('TeV').rescale('eV'), fmt='%.2e') # Tera eV
1.00e+12 eV
Differential flux, for example
>>> print(u.nit('3e8 /cm2 /sr /s /eV'))
300000000.0 1/(cm**2*s*sr*eV)
>>> pprint(u.nit('3e8 /cm2 /sr /s /eV').simplified, fmt='%.3e')
1.872e+31 s/(kg*m**4)
>>> print(nit('3 /cm2 nm sr/s J/A cd')) # I don't know this unit...
3.0 nm*cd*sr*J/(cm**2*s*A)
>>> pprint(nit('3 /cm2 nm sr/s J/A cd').simplified, fmt='%.2e') # I don't know this unit...
3.00e-05 kg*m*cd/(s**3*A)
'''
try:
# If it is known, just return the unit
return known_unit[string_units]
except KeyError:
pass
logger = _logging.getLogger(__name__ + '.nit')
if verbose:
logger.setLevel(_logging.DEBUG)
else:
logger.setLevel(_logging.WARN)
string2 = string_units
logger.debug('Input= ``%s``' % string2)
if string2.find('** ') >= 0 or string2.find(' **') >= 0:
raise ValueError('No space before/after **')
if string2.find('^ ') >= 0 or string2.find(' ^') >= 0:
raise ValueError('No space before/after ^')
if _re.search(r'[A-z]\.\d', string2) is not None:
raise ValueError('Format ambiguous. Dot is multiple or delimiter?')
if _re.search(r'\d\.[A-z]', string2) is not None:
raise ValueError('Format ambiguous. Dot is multiple or delimiter?')
string2 = string2.replace('**', '^') # ** to be ^
string2 = string2.replace('*', ' ') # * to be space.
string2 = _re.sub(r'([A-z]\s*)\.([A-z]\s*)', r'\1 \2', string2) # . to be space.
string2 = string2.replace('/', ' / ') # / can be split by split()
logger.debug('Input modified= ``%s``' % string2)
string2 = string2.split()
logger.debug('Input split= ``%s``' % str(string2))
numer = []
denom = []
isnumer = True
for string in string2:
if string == '/':
isnumer = False # Next is demoninator
else: # For usual word
if isnumer:
numer.append(string)
else:
denom.append(string)
isnumer = True # Next is numer
logger.debug('Numerator: %s' % str(numer))
logger.debug('Denominat: %s' % str(denom))
numer_u = 1
denom_u = 1
for token in numer:
numer_u *= __token2unit(token, verbose=verbose)
for token in denom:
denom_u *= __token2unit(token, verbose=verbose)
return numer_u / denom_u
nit('nT') # Define nano-tesla here.
nit('nPa') # Define nano-Pascal here.
unit = nit
""" An alias
"""
[docs]def make_unitful_function(func, unitlist, runit, offset=0):
''' Make unitful function.
:param func: Physical function (formula)
:param unitful: List of the unit in the order of the argument of ``func``.
:param runit: Unit of the returned value of ``func``.
:keyword offset: If set, the unitlist is started to be applied to
the argument after ``offset``. This is used for
when this function is applied for an object method, whose
0-th argument is the object.
When you implement a physical function (or formula),
it is usually a good idea to have internal expression
of ``float`` or ``numpy.array``.
This is usually very fast, and easy to extend.
On the other hand, you may also want to make "unitful" version,
particularly for step-by-step calculation in the interactive shell.
This function produces a function of "unitful" version
from the ``float`` or ``numpy.array`` version of function.
Example follows.
>>> import irfpy.util.unitq as u
>>> def PE(m_si, h_si):
... return 9.8 * m_si * h_si
The above function calculate the gravity potential.
It is implincitly assumes float (or numpy array) argument in the MKSA system.
>>> print(PE(1.5, 20.0))
294.0
Something with 1.5 kg mass is placed at 20.0 m, the gravity potential
is 294.0 J (at the Earth).
If you want to calculate with unit, PE cannot handle it properly.
>>> print(PE(1.5 * u.kg, 20.0 * u.m))
294.0 kg*m
The value is ok, but the unit is not in Joule, so that I cannot
convert it.
>>> pprint(PE(1500 * u.g, 2000 * u.cm), fmt='%.4e')
2.9400e+07 g*cm
The example gives more complicated situation. Again, the unit is not
consistent.
To solve this problem, I suggest to make the unitful function.
The unitful version of the function can be made via
>>> PE_u = make_unitful_function(PE, [u.kg, u.m], u.J)
The first argument is the unitless function, the second argument is
a list of unit assumed in the PE function in the order,
and the third argument is the unit returned in the KE function.
>>> print(PE_u(1500 * u.g, 2000 * u.cm).rescale(J))
294.0 J
The :make_unitful_function: returns a wrapped function (PE_u)
of the original function (PE).
When called the wrapped function, the argument is first
converted to float/np.array by stripping the given units
used in the original function (mass is 'kg' and height is 'm' here).
Then, the original function is executed with float/np.array.
Finally, the given unit (the result is 'J') is added
to the output float/np.array, to be returned by the new function.
*Use of ``offset``*
If you have formula as an method of some class.
Let's try to find out the Moon gravity calculator class.
>>> class MoonGravity:
... def __init__(self):
... self.g = 1.625 # m/s2, gravitation of Moon
... def PE(self, mass, distance):
... return self.g * mass * distance
>>> mg = MoonGravity()
>>> print(mg.PE(1.5, 20))
48.75
A new class, namely the unitful class may be generated as follows.
>>> class MoonGravity_u(MoonGravity):
... def __init__(self):
... MoonGravity.__init__(self)
... PE = make_unitful_function(MoonGravity.PE, [kg, m], J, offset=1)
>>> mg_u = MoonGravity_u()
>>> print(mg_u.PE(1500 * g, 2000 * cm))
48.75 J
Here to define (override) new PE, offset=1 is given.
This is because MoonGravity.PE method will take the object as the 0-th
argument, which should not be unit-converted.
'''
def unitful_function(*args):
if len(args) != len(unitlist) + offset:
raise IndexError('make_unitful_function: Input parameter length differ. Contact to the developer.')
newargs = list(args[:offset]) + [np.array(a.rescale(u), copy=True) for (a, u) in zip(args[offset:], unitlist)]
return runit * func(*newargs)
unitful_function.__doc__ = func.__doc__
return unitful_function
[docs]def make_unitful_function2(func, argunitarr, retunitarr):
''' Improved factory of :func:`make_unitful_function`, allowing tuple of return.
:param func: Physical function (formula)
:param argunitarr: List of unit in the order of the argument of ``func``. Units or ``None``.
:param retunitarr: List of unit of the returned value of ``func``. Units or ``None``.
When you implement a physical function (or formula),
it is usually a good idea to have internal expression
of ``float`` or ``numpy.array``.
This is usually very fast, and easy to extend.
On the other hand, you may also want to make "unitful" version,
particularly for step-by-step calculation in the interactive shell.
This function produces a function of "unitful" version
from the ``float`` or ``numpy.array`` version of function.
Example follows.
Suppose you make a function to calculate the kinetic energy and the momentum.
The argument is the mass and the velocity. Also, you may want to give an
argument for which unit conversion is not needed, and to return
a status flag, which is also not a physical unit.
The function should look like ``kene, mom, stat = enemom(mass, vel, name)``.
The list of the units of the argument is ``[u.kg, u.m / u.s, None]`` and
that for return values is ``[u.J, u.kg * u.m / u.s, None]``.
>>> def enemom(mass, vel, name):
... ke = mass * vel * vel / 2.
... mom = mass * vel
... stat = len(name) # Immitate the status.
... return ke, mom, stat
>>> print(enemom(10., 30., 'sample'))
(4500.0, 300.0, 6)
If you want to make a unitful version,
>>> enemom_u = make_unitful_function2(enemom, [kg, m / s, None], [J, kg * m / s, None])
>>> print(enemom(10. * kg, 30. * m / s, 'unitful sample')) # doctest: +NORMALIZE_WHITESPACE
(array(4500.) * kg*m**2/s**2, array(300.) * kg*m/s, 14)
.. todo::
**Known issue**. If the retunitarr is scalar, the argument
is transfered to :func:`make_unitful_function`. However,
the ``None`` argument is not supported for
:func:`make_unitful_function`.
'''
try:
len(retunitarr)
except:
return make_unitful_function(func, argunitarr, retunitarr)
def unitful_function(*args):
if len(args) != len(argunitarr):
raise IndexError('Argument length differs!')
newargs = []
for (a, u) in zip(args, argunitarr):
if u is not None:
newargs.append(np.array(a.rescale(u), copy=True))
else:
newargs.append(a)
retorgs = func(*newargs)
retvals = []
for idx in range(len(retunitarr)):
if retunitarr[idx] is not None:
retvals.append(retorgs[idx] * retunitarr[idx])
else:
retvals.append(retorgs[idx])
return retvals
unitful_function.__doc__ = func.__doc__
return unitful_function
[docs]def use_plasma():
''' Use space plasma familier system as default.
Length in cm, mass in AMU, temperature in eV, and time in s.
>>> import irfpy.util.unitq as u
>>> print((1 * u.km ** 3).simplified)
1000000000.0 m**3
>>> u.use_plasma()
>>> pprint((1 * u.km ** 3).simplified, fmt='%.1e')
1.0e+15 cm**3
>>> pprint((30000 * u.K).simplified)
2.585 eVT
>>> u.use_mksa()
>>> print((30000 * u.K).simplified)
30000.0 K
'''
_pq.set_default_units(length=cm, mass=amu, temperature=eVT, time=s)
[docs]def use_mksa():
_pq.set_default_units(system='SI')
use = set_default_units
''' Alias to set_default_units provided by qunatities
'''
[docs]def pprint(value, fmt='%.3f', *args, **kwds):
''' Print the value in easily-readable way.
:param value: A quantity
:keyword fmt: Format. Only effective if the value has shape ``()``, namely
float.
:keyword kwds: Keywords given to
:func:`irfpy.util.with_context.printoptions`
Only effective if ``value`` is array.
``precision`` and ``suppress`` are frequently used.
>>> pprint(np.sqrt(np.arange(3000)) * kg, precision=3, suppress=True) # doctest: +NORMALIZE_WHITESPACE
[ 0. 1. 1.414 ... 54.745 54.754 54.763] kg
>>> pprint(1.537732343e-10 * g, fmt='%.5e')
1.53773e-10 g
'''
number = value.n
unit = value.u
if number.shape == ():
print(fmt % number, unit)
else:
from irfpy.util.with_context import printoptions
with printoptions(*args, **kwds):
print(number, end=' ')
print(unit)
pp = pprint
""" Alias for :meth:`pprint`"""