Computational Science and Data Science

Hans Fangohr

Physical quantities (numerical value with units) in Python

The idea to combine units (such as metres, seconds and kilograms) with numbers to represent physical quantities is not new: I saw this first in the late 1980s in the MS DOS based MathCad 2.5. A similar feature has been provided by Konrad Hinsen's Scientific.Physics.PhysicalQuantities tools. Our own open-source software Nmag has provided related functionality as a very minor part of the overall package although we should have used Konrad's code instead of re-inventing the wheel. The idea of carrying information regarding the "unit of measure" has also been picked up my Microsoft's F# and probably many other tools.

Let us agree to refer to something as a physical quantity $q$ if it is made of up of a numerical value (such as 10) and (some product of powers of) units (such as metre to the power one): for example $q = 10\mathrm{m}$.

It took me a while to realise that the excellent IPython provides computation with physical quantities (based on Konrad Hinsen's code and kindly contributed by Georg Brandl). This physics module can be used as an ipython notebook extension (and then provides extra syntactic sugar), but can also be used outside the ipython notebook: either in a Python program or in a Python interpreter session. So it is a generic standalone module (at the moment just a file physics.py) which allows to work conveniently with physical quantities in Python.

I'll try to summarise basic installation and usage instructions here to provide an idea of what the tool can be used for (corrections and improvements welcome).

Using the physics module

At the time of writing (Feb 2013), there is only one file that we need to download, and that is https://bitbucket.org/birkenfeld/ipython-physics/raw/default/physics.py

In a terminal on Linux or OS X, we could just type

wget https://bitbucket.org/birkenfeld/ipython-physics/raw/default/physics.py

and the file should be downloaded to our machine into the current directory and be saved under the name physics.py. Or we need to save it from the browser. Let us assume we have done this, and the file is available in the current working directory.

Importing the physics module

The following commands are carried out within the ipython notebook but work exactly the same within a normal Python prompt (which shows >>> while waiting for user input), or in a Python program that is saved in a .py file.

We start by importing the physics module - mostly to check that the physics.py file is in the right place:

In [1]:
import physics

The main idea of the physics module can be demonstrated by using the physics.PhysicalQuantity object. The package provides the name Q as a shortcut for PhysicalQuantity, and we import that object directly:

In [2]:
from physics import Q

First steps

We can create PhysicalQuantity objects by initialising them with a number, and some units. Imagine we have the distance 10 metres to work with, and want to refer to this as s:

In [3]:
s = Q(10, 'm')
In [4]:
type(s)
Out[4]:
physics.PhysicalQuantity
In [5]:
print(s)

10 m

We now have a PhysicalQuantity object s and can carry out the usual manipulations we do with the number 10, and the units will be exposed to the same operations, and updated accordingly. For example, if we ask what the area of a square with side length s would be, we can evaluate s**2:

In [6]:
s**2
Out[6]:
100 m^2

Of course we know that 10**2=$10^2=100$, but if we square the physical quantity s, we square the number and the unit m correctly, and the output above shows 100 m^2. When dealing with more complex units and more complicated expressions, the extra information on the units can be extremely useful to detect errors.

Different PhysicalQuantity objects can be combined, for example to compute the velocity v for something that moves the distance s in the time t:

In [7]:
t = Q(2, 's')  # create physical quantity to represent the duration of 2 seconds
print(t)

2 s
In [8]:
v = s / t      # and compute the velocity when travelling distance s in time t
print(v)

5 m/s

If we do want to access the numerical value of the quantity, we can use the value attribute of the the PhysicalQuantity object:

In [9]:
v.value
Out[9]:
5

Similarly, we can ask for only the units of the quantity using .unit:

In [10]:
v.unit
Out[10]:
<PhysicalUnit m/s>

SI or not SI?

The nearly globally used system of measurements is the International System of units (so-called SI units), which goes back a long time and is based on the metre, second, kilogram, ampere, kelvin, candela and mole as the base units.

As long as we express all quantities in these base units, and compute other entities (correctly) purely based on their numerical value, we can rely on the resulting quantity having the right number for these base units.

Example (in base units)

Here is an example: we measure the distance in metres, say 10,000 metres (the base unit for distance is metre), and we establish that it took as 1800 seconds (the base unit for time is seconds) to travel that distance by bicycle. We can compute the average velocity by dividing the distance by time, and get the average velocity in the base units:

In [11]:
s1 = 10000     # we know this is expressed in metres, even though s1 does not know.
t1 = 1800      # we knows this is expressed in seconds. Of course t1 does not know.
v1 = s1 / t1   # we know that the result will be expressed in m/s. And v1 does not know.
print(v1)

5

So our average velocity, while travelling 10,000 metres in 1,800 seconds, is 5 metres per second. The calculation only returns the value 5 but we know the units are metres per second because this is the way to express velocity in base units).

Example (not in base units)

We tend to speak about velocities of vehicles in kilometre per hour, km/h, (or miles per hour) which are not SI base units. So we need to convert the result of 5 m/s into, say, km/h. We can do this manually by starting from 5 m/s, and then we substitute m with km/1000 and substitute s with h/3600:

$$ 5 \frac{m}{s} = 5 \frac{m}{1} \frac{1}{s} = 5 \frac{km}{1000}\frac{3600}{h} = 5 \frac{km}{h}3.6 = 5\cdot 3.6 \frac{km}{h} = 18 \frac{km}{h}$$

We can repeat the same exercise using the PhysicalQuantity objects, which will make the conversion to km/h much easier:

In [12]:
s2 = Q(10000, 'm')
t2 = Q(1800, 's')
v2 = s2 / t2
print(v2)

5 m/s

The PhysicalQuantity provides a convenient convert method, that can change the way the quantity is expressed into any (compatible) units. So to learn what v2 is in km/h, we can use

In [13]:
v2.convert('km/h')
print(v2)

18 km/h

We can change the way that v2 represents its quantity as often as we like. In particular, we may go back to m/s:

In [14]:
v2.convert('m/s')
print(v2)

5 m/s

Or to nanometre per nanosecond (useful in a nanotechnology context, for example the velocity of a read/write head moving relative to the rotating disk in a hard drive is of the order of nanometres per nanosecond):

In [15]:
v2.convert('nm/ns')
print(v2)

5 nm/ns

Or to miles per hour - admittedly we need to know that miles are abbreviated mi within the physics.py file. The best way to find out is to look into the source code if necessary (see below). So here we go in (British) miles per hour:

In [16]:
v2.convert('mi/h')
print(v2)

11.184681 mi/h

The various unit identifiers such as mi for miles are created in the physics.py file, often with a useful comment. So if we have that file handy, we could also just grep for mile and see what names come up:

In [17]:
!grep -i mile physics.py

_addUnit('mi', '5280.*ft', '(British) mile')
_addUnit('nmi', '1852.*m', 'Nautical mile')

Convert quantity to SI units

Finally, the PhysicalQuantity object has a convenient property base which returns the quantity expressed in base units which are SI units for the physics package:

In [18]:
v2  # we had previously expressed this velocity in units of miles per hour
Out[18]:
11.184681 mi/h
In [19]:
v2.base
Out[19]:
5 m/s

And for those who are missing the CGS system, there is a also a property cgs that expresses the physical quantity in base units of Centimetres, Grams and Seconds (CGS):

In [20]:
v2.cgs
Out[20]:
500 cm/s

PhysicalQuantities stick to their units

We note that PhysicalQuantity objects remember their representation: if we create a distance as an Astronomical Unit au (roughly the distance Earth to Sun), then we do so because this is a sensible length scale for the given problem. It thus makes sense (as a design decision for the behaviour of PhysicalQuantities) to stick to the representation of the quantity in the same units when we use it or a derived quantity. For example

In [21]:
longdistance = Q(2, 'au')
print("Twice the distance Earth-Sun is %s." % longdistance)
print("which is %s in SI base units." % longdistance.base)
print("But an astronomer may argue that %s looks more friendly." % longdistance)

Twice the distance Earth-Sun is 2 au.
which is 2.9919574e+11 m in SI base units.
But an astronomer may argue that 2 au looks more friendly.

Predefined physics constants

Now that we are looking at astronomical distances, we may as well ask: how long does it take light to travel 2 astronomical units? The physics package has a number of important physical constants available which may help us here and which we can access:

In [22]:
physics._constants
Out[22]:
[('pi', 3.141592653589793),
 ('e', 2.718281828459045),
 ('c0', 2.9979246e+08 m/s),
 ('mu0', 1.2566371e-06 kg*m/A^2/s^2),
 ('eps0', 8.8541878e-12 A^2*s^4/kg/m^3),
 ('Grav', 6.67384e-11 m^3/s^2/kg),
 ('hpl', 6.6260696e-34 s*J),
 ('hbar', 1.0545717e-34 s*J),
 ('e0', 1.6021766e-19 C),
 ('me', 9.1093829e-31 kg),
 ('mp', 1.6726218e-27 kg),
 ('mn', 1.6749274e-27 kg),
 ('NA', 6.0221413e+23 1/mol),
 ('kb', 1.3806488e-23 J/K),
 ('g0', 9.80665 m/s^2),
 ('R', 8.3144621 J/K/mol),
 ('alpha', 0.0072973525698),
 ('Ry', 10973732 1/m),
 ('mu_n', -9.6623647e-27 J/T),
 ('gamma', 183.24718 MHz/T),
 ('h0', 0.704),
 ('sigmaT', 6.652453e-29 m^2)]

One of those is the velocity of light 'c0', so we can divide the distance by the velocity to get the travel time:

In [23]:
longdistance / Q(1, 'c0')
Out[23]:
2 au/c0

The resulting PhysicalQuantity sticks to the use of c0 and au as the units, because we started from them. Of course we can convert this into any (compatible) units we like. The easiest way is to use the base property:

In [24]:
traveltime = longdistance / Q(1, 'c0')
traveltime.base
Out[24]:
998.00957 s

So the light needs about 1000 seconds to travel 2 astronomical units. How much is that in minutes?

In [25]:
traveltime.convert('min')
print traveltime

16.633493 min

We could of course have defined the velocity of light ourselves, for example

In [26]:
c = Q(2.99792458e8, 'm/s')
print (longdistance / c).base

998.00957 s

Equality of PhysicalQuantities

For clarity, let's create two objects x and y that describe the same velocity but expressed as 18 km/h and 5 m/s, respectively:

In [27]:
x = Q(5, 'm/s')
y = Q(18, 'km/h')
print("x = %s" % x)
print("y = %s" % y)

x = 5 m/s
y = 18 km/h

The velocities appear different when printed, yet they describe the same velocity, i.e. they represent the same quantity. The PhysicalQuantity comparison confirms this:

In [28]:
x == y
Out[28]:
True

Another way to compare quantities is to refer back to their unique representation using base units:

In [29]:
print(x.base)
print(y.base)

5 m/s
5 m/s

Using the physics extension in the IPython Notebook

Installation IPython-Physics

We need to install the physics packages as an IPython extension using this command:

In [30]:
%install_ext https://bitbucket.org/birkenfeld/ipython-physics/raw/default/physics.py

Installed physics.py. To use it, type:
  %load_ext physics

(If this fails, or your IPython version is older, one can also do a manual install using the instructions from https://bitbucket.org/birkenfeld/ipython-physics).

I had to start a new IPython Notebook to be able to load the physics extension after running the install_ext command.

Activation of physics extension in IPython Notebook

Let us activate the installed physics extension to be available in our IPython Notebook session:

In [31]:
%load_ext physics

Unit calculation and physics extensions activated.

Syntactic niceties of Physics extension in IPython Notebook

The functionality of PhysicalQuantities in the IPython Notebook is exactly the same as without it, and all the examples shown above will work in the Ipython Notebook.

The extra syntactic sugar that IPython offers, is that instead of x = Q(5, 'm') we can simply write

In [32]:
x = 5 m
In [33]:
print(x)

5 m

Another additional syntax feature is the // operator which is used to express the quantity on the left of the operator // in units given as the operand to the right of //. Here are some examples:

In [34]:
x // cm
Out[34]:
500 cm
In [35]:
x // km
Out[35]:
0.005 km

which show how the quantity x can be expressed in units of cm or km.

Note that for this convenient notation, the line has to complete immediately after the units: any spaces or comments being appended within the same line result in a Syntax error.

We further note that we may need to write statements that make use of this syntactic sugar in separate notebook cells (to avoid a Syntax error with the current version of the physics module.)

Application example in IPython Notebook: kinetic and potential energy (and Tom Daley)

We review the features introduced above by computing the kinetic energy of a car that has a mass of 1000 kg and moves at a velocity of 50 km/h:

In [36]:
m = 1000 kg
In [37]:
v = 50 km/h

The PhysicalQuantity object v remembers that it has been created with the value 50 and the unit km/h:

In [38]:
v
Out[38]:
50 km/h

We can convert it into SI units (of metre and seconds here), using the .base attribute which returns the same physical quantity in SI base units:

In [39]:
v.base
Out[39]:
13.888889 m/s

Or we can use the convert method that will change the way the object v is presented from km/h to m/s:

In [40]:
v.convert('m/s')
v
Out[40]:
13.888889 m/s

The kinetic energy $E_\mathrm{k}$ is $E_\mathrm{k} = \frac{1}{2} m v^2$, and we can conveniently compute this as

In [41]:
Ek = 0.5 * m * v **2
Ek
Out[41]:
96450.617 m^2*kg/s^2

The physics extension has used the information it had about the mass [kg] and the velocity [m/s] to work out the correct units for the kinetic energy that contains the mass [kg] and the velocity squared [m^2 / s^2].

Or to get this in kilo Joule (kJ):

In [42]:
Ek // kJ
Out[42]:
96.450617 kJ

We may now ask from what height the same car would have to fall from to reach the same kinetic energy (on the surface of planet Earth). We do this by comparing the potential energy $E_\mathrm{p} = m g h$ with the kinetic energy $E_\mathrm{k}$ that we have at a velocity of 50 km/h, where $g = 9.81\mathrm{kg/s^2}$ is the Earth's gravitational acceleration at the surface of the planet. And thus $h = \frac{E_{\mathrm{k}}}{m g}$:

In [43]:
g = 9.81 m/s^2
In [44]:
h = Ek / (m * g)
In [45]:
h
Out[45]:
9.8318672 m

What does this quantity $h$ mean? It means that the kinetic energy of a car travelling with 50 km/h is the same as the energy of that car falling (in vacuum) from a height of $\approx 9.83\\mathrm{m}$. This may appear to be a fairly large height: jumping from a ten metre tower is generally seen as dangerous (unless one is jumping into water, and even then this is not risk free). However, what it actually shows is that the travelling at (only) 50 km/h is also a somewhat risky procedure.

The fact that our calculation not only returned a numerical value (i.e. 9.83) but also a unit (i.e. metres) gives some extra confidence in the result: at least we have not made a mistake which would lead to the result having the wrong units (and almost certainly the wrong value as well).

We should say for completeness, that there are better ways of obtaining this result: for example we can see that the mass of the object (and thus the total kinetic energy) does not actually matter for the question posed, and that generally the fall height $h$ relates to a travelling velocity $v$ through

$$ h = \frac{v^2}{2g}$$

so that for the particular numbers we have

In [46]:
v**2 / (2 * g)
Out[46]:
9.8318672 m

Another way of cross checking this result is that for objects under the constant acceleration $g$, the velocity changes as $v(t) = g t$, and the travelled distance is given by $s(t) = \frac{1}{2} g t^2$. Solving the latter equation for $t$ provides $t(s) = \sqrt{2\frac{s}{g}}$, i.e. for our example where $s=h$:

In [47]:
t = sqrt(2*h/g)
In [48]:
t
Out[48]:
1.4157889 s

This means it takes 1.416 seconds to fall a height of 9.83 metres. The velocity is then obtain as $v(t) = g t$

In [49]:
v2 = g * t
In [50]:
v2
Out[50]:
13.888889 m/s
In [51]:
v2.convert('km/h')
In [52]:
v2
Out[52]:
50 km/h

The velocity v2 that we have computed independently based on the fall height $h$ matches the velocity v we started from to determine the right fall height; giving us more confidence in the correctness of the procedure.

A quick reminder regarding the material presented in this section: the IPython Notebook notation (i.e. x = 5 m instead of x = Q(5, 'm')) can only be used within the Ipython Notebook. To save a Python program that carries out calculations of quantities, we need to use the PhysicalQuantity constructor Q as shown in the first part of this post.

Summary

In the same way that computer-assisted symbolic calculation is often undervalued, I believe that we as computational scientists and engineers should make more use of calculations with physical quantities that carry information about the units of a numerical value. I have certainly spend a lot of time getting units right (often getting them wrong first), and see others making similar mistakes in long handwritten calculations. Mixup of units can lead to the failure of major projects such as the so-called metric mixup of the Mars Climate Orbiter leading to its crash in 1999.

In practical simulation and modelling we are often given physical quantities from experimentalists or sensors or theoretical scientists that need to be converted, for example from non-SI to SI units or into simulation units. Even independent of computer simulations, we regularly need to convert physical quantities between different unit systems for calculations and further processing of data.

Where appropriate and possible in our software, we should explicitly carry around the units in which a quantity is expressed together with the numerical value (for example using the physics module) - instead of assuming we know in which units a quantity is expressed. This avoids errors introduced by incorrect/missing documentation or misunderstandings regarding the units for the numerical values. Furthermore, instead of working out on paper the conversion factor for a numerical value based on one set of units to another set of units, we should delegate this algorithmic task to the computer. This makes the process more reproducible, more flexible, and documents the conversion.

Software versions used: the ipython-physics repository commit used for the testing above was 66ece9cb6b0916502ccf8f9984c7a1a28e8b1a17 from 22 Oct 2012. The IPython version was 0.13.1.

You can download this notebook.

Comments