** Using units in python
:PROPERTIES:
:categories: python, units
:date: 2013/01/19 09:00:00
:updated: 2013/03/23 09:45:20
:END:
[[http://matlab.cheme.cmu.edu/2011/08/05/using-cmu-units-in-matlab-for-basic-calculations/][Units in Matlab]]
I think an essential feature in an engineering computational environment is properly handling units and unit conversions. Mathcad supports that pretty well. I wrote a [[https://github.com/jkitchin/matlab-cmu][package]] for doing it in Matlab. Today I am going to explore units in python. Here are some of the packages that I have found which support units to some extent
1. http://pypi.python.org/pypi/units/
2. http://packages.python.org/quantities/user/tutorial.html
3. http://dirac.cnrs-orleans.fr/ScientificPython/ScientificPythonManual/Scientific.Physics.PhysicalQuantities-module.html
4. http://home.scarlet.be/be052320/Unum.html
5. https://simtk.org/home/python_units
6. http://docs.enthought.com/scimath/units/intro.html
The last one looks most promising.
#+BEGIN_SRC python
import numpy as np
from scimath.units.volume import liter
from scimath.units.substance import mol
q = np.array([1, 2, 3]) * mol
print q
P = q / liter
print P
#+END_SRC
#+RESULTS:
: [1.0*mol 2.0*mol 3.0*mol]
: [1000.0*m**-3*mol 2000.0*m**-3*mol 3000.0*m**-3*mol]
That doesn't look too bad. It is a little clunky to have to import every unit, and it is clear the package is saving everything in SI units by default. Let us try to solve an equation.
Find the time that solves this equation.
$0.01 = C_{A0} e^{-kt}$
First we solve without units. That way we know the answer.
#+BEGIN_SRC python
import numpy as np
from scipy.optimize import fsolve
CA0 = 1.0 # mol/L
CA = 0.01 # mol/L
k = 1.0 # 1/s
def func(t):
z = CA - CA0 * np.exp(-k*t)
return z
t0 = 2.3
t, = fsolve(func, t0)
print 't = {0:1.2f} seconds'.format(t)
#+END_SRC
#+RESULTS:
: t = 4.61 seconds
Now, with units. I note here that I tried the obvious thing of just importing the units, and adding them on, but the package is unable to work with floats that have units. For some functions, there must be an ndarray with units which is practically what the UnitScalar code below does.
#+BEGIN_SRC python
import numpy as np
from scipy.optimize import fsolve
from scimath.units.volume import liter
from scimath.units.substance import mol
from scimath.units.time import second
from scimath.units.api import has_units, UnitScalar
CA0 = UnitScalar(1.0, units = mol / liter)
CA = UnitScalar(0.01, units = mol / liter)
k = UnitScalar(1.0, units = 1 / second)
@has_units(inputs="t::units=s",
outputs="result::units=mol/liter")
def func(t):
z = CA - CA0 * float(np.exp(-k*t))
return z
t0 = UnitScalar(2.3, units = second)
t, = fsolve(func, t0)
print 't = {0:1.2f} seconds'.format(t)
print type(t)
#+END_SRC
#+RESULTS:
: t = 4.61 seconds
:
This is some heavy syntax that in the end does not preserve the units. In my Matlab package, we had to "wrap" many functions like fsolve so they would preserve units. Clearly this package will need that as well. Overall, in its current implementation this package does not do what I would expect all the time.[fn:units]
[fn:units] Then again no package does yet!