The Kitchin Research Group

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

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

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

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)

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.

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)

t = 4.61 seconds
<type 'numpy.float64'>

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

Copyright (C) 2013 by John Kitchin. See the License for information about copying.

org-mode source