How to work with physical units in Python?¶
Autor: Eduardo Martín Calleja
Imports and references¶
from __future__ import division
import quantities as pq
import numpy as np
import inspect
# This IPython magic generates a table with version information
#https://github.com/jrjohansson/version_information
%load_ext version_information
%version_information numpy, quantities
There are several Python packages for working with units and physical constants. I tried some of them and my personal preference is the package quantities I find its use particularly simple and natural, and the level of functionality is enough for me, so it will be that I will use in the following examples.
Why use physical units in Python?¶
The first question that may arise is: why I need my preferred programming language to support physical units?. Well, suppose that at a certain time I assign to a variable the value of a length equal to 3 measured in meters. Another time I create another variable with a length of 5 feet extent. And later I add both quantities. In the absence of a control of the physical units involved, I'd get a result of 8, meaningless and wrong in any system of units. Imagine instead that I'm working with quantities, I would write:
x = 3 * pq.m
y = 5 * pq.foot
print x + y
As you can see, we get a correct result in meters.
But ... How could it happen to add meters and feet?. Well, in 1998 the spacecraft Mars Climate Orbiter (MCO) of NASA disintegrated in the Martian atmosphere because the software in the Earth based tracking station made use of the Anglo-Saxon system of units (inches, feet and pounds) while the software embarked on the ship used, according to specifications, the metric system (millimeters, meters, kilometers and kilos). No one made the necessary conversions.
Of course there are many other possibilities for error when working with numbers expressing quantities with physical units. In particular the loss of dimensional consistency. Suppose for example that we intend to add two quantities inadvertently, one representing a magnitude of force (in Newtons) and other of energy (in Joules):
x = 3.5 * pq.newton
y = 0.7 * pq.joule
try:
x + y # we get an error if the operation is invalid with this units
except Exception as e:
print u'*** Attention! ***: ',e
In this case quantities has generated an exception that warns us of the dimensional inconsistency of the operation we want to perform.
The usefulness of quantities goes much further, by keeping track of the composite units used to express the result of a calculation and automatically performing the appropriate conversions so that the result is correct. This means we can enter the values of the variables in the units we are most comfortable with at each time:
x = 0.3 * pq.km # length in Km
t = 250 * pq.ms # time in milliseconds
accel = x/t**2
mass = 7 * pq.g # mass in grams
force = mass * accel
print force
We can ask that the result is expressed in the base units (SI by default):
print force.simplified
Or to be expressed in any unit derived, while respecting dimensional consistency:
print force.rescale(pq.newton)
print force.rescale(pq.dyne)
What is most important is that the correct units will go along all intermediate calculations, appropriate conversions being realized automatically, and letting us know if we ever make a dimensional mistake. Furthermore, the final outcome can be expressed in the units we want, provided they are dimensionally compatible with the units computed.
The most used units in astronomy¶
We have seen before some simple examples of using physical units with the package quantities. Its use in astronomy can be particularly useful due to the diversity of units. Below we will indicate some of the most common. In addition We will show the expression "simplified" in order to check the numerical value in SI units:
# Time and distance units
print pq.year # time in years
print pq.au, pq.au.simplified # distance in astronomical units
print pq.light_year, pq.light_year.simplified # distance in light-years
print pq.pc, pq.pc.simplified # distance in parsecs
# Temperature units
print pq.K
print pq.celsius
But be very careful because quantities does not convert between scales of temperatures and absolute temperatures in degrees, as this means considering a reference point (absolute zero) for which the package is not ready. Temperatures are always considered as temperature differences.
# In this example we see that the conversion is not performed
temp = 345 * pq.celsius
print temp.rescale(pq.K)
# Neither is correct the conversion between Farhrenheit and Celsius scales
# However, if the value is considered as an increase in temperature,
# then, the result is the correct answer
temp = 170 * pq.fahrenheit
print temp.rescale(pq.celsius)
# Measurement of angles
print pq.deg
print pq.rad
print pq.arcmin
print pq.arcsec
Examples of using angular units:
print (np.pi/2 * pq.rad).rescale(pq.deg)
print (1 * pq.deg).rescale(pq.arcmin)
You can define new units that do not exist in the package:
# Definition of new units
Mly = pq.UnitQuantity('megalight_year', 1e6*pq.ly, symbol = 'Mly')
Mpc = pq.UnitQuantity('megaparsec', 1e6*pq.pc, symbol = 'Mpc')
For example, how many light years are there in a megaparsec?
print (1 * Mpc).rescale(pq.ly)
The previous example shows that the new unit is not part of the quantities namespace, ie we write Mpc, not pq.Mpc
Constants¶
The package quantities also presents a wide range of physical constants. As an example we quote the following, of interest in astronomy:
print pq.c, pq.c.simplified # Speed of light
print pq.constants.G, pq.constants.G.simplified # Newtonian constant of gravitation
And you can define new constants. I personally use the following ones:
# Definition of new constants:
earth_radius = 6378160.0 * pq.m
moon_radius = 1740000.0 * pq.m
sun_radius = 695000000.0 * pq.m
earth_mass = 5.97219e24 * pq.kg
sun_mass = 332946 * earth_mass
Hubble_constant = 67.80 * (pq.km/pq.s)/(1e6*pq.pc)
NumPy Support¶
A major strength of quantities is that it is fully compatible with NumPy arrays. You can define the type of units in an array:
xa = np.arange(10) * pq.N
Each element of the above array will have units in Newtons
print xa[5]
And array operations are supported: In the following example, by dividing the former array by a magnitude of mass, all the elements of the array will have acceleration units:
print (xa / (3.0 * pq.kg)).simplified
Sometimes it is necessary to obtain the numerical value of the variable without including the unit type. This can be done in two ways, depending on that we want to obtain a numeric value or an array of NumPy:
# A single numerical value
f = 30 * pq.N
print f.item(), type(f.item())
# An array of values
print xa.magnitude, type(f.magnitude)
Exploring the contents of quantities¶
The quantities.units package includes a large number of units arranged in a series of sub-modules, one for each type of unit. To list the different sub-modules we can write as follows:
for nombre, datos in inspect.getmembers(pq.units, inspect.ismodule):
print nombre, ',',
Then you can get a list of the names of the units in each package. For example, to see what are the units of type 'mass':
for u in dir(pq.mass):
print u,',',
The other important package is quantities.constants:
for nombre, datos in inspect.getmembers(pq.constants, inspect.ismodule):
print nombre, ',',
And to see for example the content of the module of astronomical constants:
dir(pq.constants.astronomy)
No hay comentarios:
Publicar un comentario