Source code for irfpy.util.unitq

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

``irfpy.util.unitq`` handles the unit (dimension) based on the quantities package.
See http://pythonhosted.org/quantities/ for detailed documentation on the quantities package.

``irfpy.util.unitq`` adds some functionalies:

- Shorter names (aliases) of the dimension.

 - :meth:`r` (alias of :meth:`rescale`)

 - :meth:`n` (get the value - remove units)

 - :meth:`s` (get the simplified form based on the current system -- default MKSA system)

 - :meth:`u` (get the unit in the string format)

- Simple conversion of the string formatted unit to the unit object.
  (Functions :func:`unit`, or :func:`nit`)

- Some new units commonly used in plasma physics. (For example, :attr:`cu_flux`.)

- Dynamic definition of the new unit with 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`.)

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

- Converting the "list of Quantities" to the single "Quantities" instance.
  (:func:`asarray`)


*Usage*

It is recommended to start with

>>> import irfpy.util.unitq as u

to use the unitq module.

Then, you can make a quantity of "8.57 m" simply as

>>> print(8.57 * u.m)
8.57 m

Numerical calculation is simply done.

>>> l = 8.57 * u.m   # 8.57 meter
>>> t = 2.48 * u.s   # 2.48 seconds
>>> v = l / t    # Get the speed
>>> print('{:.3f}'.format(v))
3.456 m/s

Or you can make a unit object using :func:`unit` function (or :func:`nit`, with the module name ``u``, it can be called ``u.nit``).

>>> print(u.nit('300 MPa'))    # 300 Mega pascal
300.0 MPa

The attribute ``s`` (or ``simplified`` for long) convert the unit.

>>> print(u.nit('3.6 km/h').s)
1.0 m/s

You can remove the unit to get the normal numpy ndarray with the attribute of ``n``.

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

The new array is a new instance, so changes are not propageted.

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

The solid angle is supported. But when you simplify to MKSA, it will go to dimensionless.

>>> 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 is used both energy and temperature. This could be ambiguous.
The default ``eV`` is for energy. A new unit of ``eV_asT`` (or ``eVT`` in short) is defined.
The function :func:`r` rescales the unit.

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

In summary, the following functions / attributes can be used.

* ``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.

In the following, a series of sample is seen.

Constant with units are defined under the namespace of "k".
For example the Planck constant can be obtained as.

>>> h = u.k.Planck_constant   # Planck constant
>>> print(h)
1 h (Planck_constant)

>>> print(h.n)   # Get the value in the unit of Planck_constant
1.0
>>> print(h.u)   # The unit is ``h`` (Planck_constant).
h

Let's see this value in the MKSA system with simplified attribute.

>>> print(h.simplified)
6.62606896e-34 kg*m**2/s

>>> print(h.ns)   # This expression simplifies the quantity, then return as a value.
6.62606896e-34
>>> print(h.us)    # This returns the unit of simplified quantity.
kg*m**2/s

Rescaling to another unit is also done using the function ``r``, ``nr`` and ``ur``.

>>> print('{:.3e}'.format(h.r('g*cm^2/s')))   # Rescale to the g cm^2 / s
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

*Dynamic definition of the unit*

Dynamic definition means that based on the defined unit, SI prefix is added and dynamically defined.
For example, Giga year is not defined by default.

>>> print(u.Gyr)    #  doctest: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
...
AttributeError: module 'irfpy.util.unitq' has no attribute 'Gyr'

However, once you use :func:`unit` function, "Gyr" is defined.
>>> three_giga_year = u.nit('3 Gyr')
>>> print(three_giga_year)
3000000000.0 yr
>>> print(u.Gyr)    # Gyr was defined automatically
1 Gyr

*Convert list of unit array to the single unit array*

Sometimes, one need to iterate some calculation then you may get a list of unit array.
As a simple example,

>>> l = [1, 2, 3] * u.m
>>> t = 1 * u.s
>>> print(l / t)
[1. 2. 3.] m/s

For this simple example of calculating the velocity, the above division is sufficient.
But if a calculation becomes more complicated, iteration may not be avoided.

Such a case, one can rewrite the velocity calculation as

>>> v = [l_ / t for l_ in l]    # This is not equal to "l / t".
>>> print(v)
[array(1.) * m/s, array(2.) * m/s, array(3.) * m/s]

The result is not very easy to handle as unit array.

>>> print(v.rescale('cm/s'))   # Rescale does not work.   doctest: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
...
AttributeError: 'list' object has no attribute 'rescale'

The function :func:`asarray` can be used to convert a list of unit arrays to a single unit array.

>>> v2 = asarray(v)
>>> print(v2)
[1. 2. 3.] m/s
>>> print(v2.rescale('cm/s'))
[100. 200. 300.] cm/s
'''
import logging as _logging
import re as _re

#import numpy as _np

import quantities as _pq
import quantities.constants as k
from quantities import *

_logger = _logging.getLogger(__name__)

_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]def asarray(a, unit=None): """ From the input array-like to make as unit array. :param a: Array like input. All elements should be values with units, equivalent for all. :type a: Array like. :keyword unit: The unit to be used. It should not be a text representation. Use to wrap ``u.nit('m')`` instead of 'm'. :return: An array fully with a single unit. A simple one >>> from irfpy.util import unitq as u >>> a_in = [3, 4, 5] * u.m >>> a_out = u.asarray(a_in, unit=u.cm) # This should work as if the rescale function >>> print(a_out) [300. 400. 500.] cm >>> a_out = u.asarray([3, 4, 5], unit=u.km) # This also should work as if the ``[3, 4, 5] * u.km`` >>> print(a_out) [3. 4. 5.] km >>> a_in = [3 * u.cm, 4 * u.m, 5 * u.km] # A list. Input unit shall be connected to each component. >>> a_out = asarray(a_in) # In this case, the input array shall be converted to the first element unit, and converted to the unit array >>> print(a_out) [3.e+00 4.e+02 5.e+05] cm >>> a_in = [[3 * u.cm, 4 * u.m, 5 * u.mm], [15, 20, 30] * u.mm] # (3, 2) shape >>> a_out = asarray(a_in, unit=u.nit('mm')) >>> print(a_out) [[ 30. 4000. 5.] [ 15. 20. 30.]] mm """ # First, create a copied array of input. This eliminates the units. a_arr = np.array(a) _logger.debug(a_arr) if a_arr.ndim == 0: # Scalar case return _asarray_0d(a, unit=unit) elif a_arr.ndim == 1: # 1 dim return _asarray_1d(a, unit=unit) else: return _asarray_nd(a, unit=unit)
def _asarray_1d(a, unit=None): """ 1d input array to be converted to the unit array. :param a: :param unit: :return: >>> from irfpy.util import unitq as u Given normal list + no unit given -> Dimensionless unit >>> print(asarray([1, 3, 5])) [1. 3. 5.] dimensionless Given normal list + unit specified -> Unit array >>> print(asarray([1, 3, 5], unit=u.km / u.s)) [1. 3. 5.] km/s Given array + unit specified >>> print(asarray(np.array([1, 3, 5]), unit=u.cm**2)) [1. 3. 5.] cm**2 Given unitarray + different unit >>> print(asarray([1, 3, 5] * u.m, unit=u.cm)) [100. 300. 500.] cm Given unitarray + No dimension > As is >>> print(asarray([1, 3, 5] * u.cm / u.h)) [1. 3. 5.] cm/h Given list of unit scalar, no dimension -> Rescaled to the first element >>> print(asarray([1 * u.m, 3 * u.km, 5 * u.cm])) [1.e+00 3.e+03 5.e-02] m Given list of unit scalar, dimension specified >>> print(asarray([1 * u.m, 3 * u.km, 5 * u.cm], unit=u.mm)) [1.e+03 3.e+06 5.e+01] mm """ # Does a have unit? try: in_unit = a[0].units except: in_unit = None if in_unit is None: # a has no specific dimension if unit is None: # In and out unit is None. -> Dimensionless value to be returned. unit = dimensionless return np.asarray(a) * unit # In unit is None, out unit is specified. else: # a has dimension if unit is None: # In specified (a is unit array), but out not -> Keep in. unit = in_unit b = [_.nr(unit) for _ in a] * unit return b def _asarray_nd(a, unit=None): """ >>> from irfpy.util import unitq as u Given normal list + no unit given -> Dimensionless unit >>> print(asarray([[1, 3, 5], [7, 6, 4]])) [[1. 3. 5.] [7. 6. 4.]] dimensionless Given normal list + unit specified -> Unit array >>> print(asarray([[1, 3, 5], np.array([7, 6, 4])], unit=u.km / u.s)) [[1. 3. 5.] [7. 6. 4.]] km/s Given array + unit specified >>> print(asarray(np.array([[1, 3, 5], [7, 6, 4]]), unit=u.cm**2)) [[1. 3. 5.] [7. 6. 4.]] cm**2 Given unitarray + different unit >>> print(asarray([[1, 3, 5], [7, 6, 4]] * u.m, unit=u.cm)) [[100. 300. 500.] [700. 600. 400.]] cm Given unitarray + No dimension >>> print(asarray([[1, 3, 5], [7, 6, 4]] * u.cm / u.h)) [[1. 3. 5.] [7. 6. 4.]] cm/h Given list of unit scalar, no dimension -> Rescaled to the first element >>> print(asarray([[1 * u.m, 3 * u.km, 5 * u.cm], [7, 6, 4] * u.cm])) [[1.e+00 3.e+03 5.e-02] [7.e-02 6.e-02 4.e-02]] m Given list of unit scalar, dimension specified >>> print(asarray([[1 * u.m, 3 * u.km, 5 * u.cm], [7 * u.mm, 6 * u.km, 4 * u.cm]], unit=u.mm)) [[1.e+03 3.e+06 5.e+01] [7.e+00 6.e+06 4.e+01]] mm """ ### a cannot be flatteded, as it is not necessarily ndarray. It is not impossible to convert to ndarray keeping the dimension. #a_ = np.asanyarray(a, dtype=object, copy=True) # Copy the input array to object array, keeping the information of a_ = np.array(a) # Copy the input array, but this loses unit information. a_shape = a_.shape # The shape. a_ndim = a_.ndim # The dimension it_a = np.nditer(a_, flags=['multi_index']) b = np.zeros_like(a_) # Output array. Initialized by zero ### Extremely slow possibly, but it should work while not it_a.finished: idx = it_a.multi_index val = a for _ in idx: # Multi index recursively get until each element val = val[_] _logger.debug(str(idx) + str(val)) # This should be the element if unit is None: # Out unit should only initialized once try: # Refer to the input unit for the first element in_unit = val.units except: in_unit = None if in_unit is None: unit = dimensionless else: unit = in_unit try: b[idx] = val.nr(unit) # If the input has unit, it works except AttributeError: b[idx] = val # Otherwise, it is just a value. it_a.iternext() return b * unit def _asarray_0d(a, unit=None): """ 0d input array/scalar to be converted to the unit-array. # Given a scalar. >>> from irfpy.util import unitq as u >>> print(asarray(5)) 5.0 dimensionless >>> print(asarray(u.nit('5 km'))) 5.0 km >>> print(asarray(5, unit=u.m)) 5.0 m >>> print(asarray(5 * m, unit=u.cm)) 500.0 cm """ # Does a have unit? try: in_unit = a.units except: in_unit = None if in_unit is None: if unit is None: # In and out unit is None. -> Dimensionless value to be returned. unit = dimensionless return np.asarray(a) * unit # In unit is None, out unit is specified. else: if unit is None: # In specified (a is unit array), but out not -> Keep in. return a # Return as it is. else: return a.r(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`"""