Trabajando con unidades y constantes físicas en Python

Autor: Eduardo Martín Calleja

Importaciones y referencias¶

In [1]:

from __future__ import division

import quantities as pq
import numpy as np
import inspect

# Generar un cuadro con versiones de las librerías utilizadas en este notebook
#https://github.com/jrjohansson/version_information
%load_ext version_information
%version_information numpy, quantities

Out[1]:

Software	Version
Python	2.7.9 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]
IPython	2.3.1
OS	Linux 3.13.0 45 generic x86_64 with debian jessie sid
numpy	1.9.1
quantities	0.10.1
Thu Feb 05 12:26:14 2015 CET

¿Por qué utilizar unidades físicas en Python?¶

Existen varios paquetes de Python que permiten trabajar con unidades y constantes físicas. He probado algunos de ellos y mi preferencia personal es el paquete quantities. Su utilización me resulta particularmente sencilla y natural, y para mí posee el nivel de funcionalidad suficiente, por lo que será el que utilizaré en los siguientes ejemplos.

La primera cuestión que se puede plantear es: ¿por que necesito que mi lenguaje de programación preferido soporte unidades físicas?. Bien, supongamos que en un momento determinado asigno a una variable el valor de una longitud igual a 3 medida en metros. En otro momento creo otra variable con una longitud de 5 medida en pies. Y más tarde sumo ambas cantidades. Sin un control de unidades físicas obtendría un resultado de 8, carente de sentido y erroneo en cualquier sistema de unidades. Imaginemos por el contrario que estoy trabajando con quantities:

In [2]:

x = 3 * pq.m
y = 5 * pq.foot
print x + y

4.524 m

Como se ve, obtenemos un resultado correcto expresado en metros.

Pero... ¿a quién va a ocurrirsele sumar metros y pies?. Bueno, en el año 1998 la sonda espacial Mars Climate Orbiter (MCO) de la NASA se desintegró en la atmósfera marciana porque el software de la estación de seguimiento en Tierra hacía uso del sistema de unidades anglosajón (pulgadas, pies y libras) mientras que el software embarcado en la nave utilizaba, de acuerdo con las especificaciones, el sistema métrico decimal (milímetros, metros, kilómetros y kilos), sin que se llevaran a cabo las conversiones requeridas.

Por supuesto existen muchas otras posibilidades de error cuando se trabaja con cantidades que expresan magnitudes con unidades físicas. En particular puede no respetarse la coherencia dimensional. Supongamos por ejemplo que inadvertidamente pretendemos sumar dos cantidades, una que representa una magnitud de fuerza (en Newtons) y otra de energía (en Julios):

In [3]:

x = 3.5 * pq.newton
y = 0.7 * pq.joule
try:
    x + y # da un error si se intenta realizar una operación no válida con esas unidades
except Exception as e:
    print u'*** ¡Atención! ***: ',e

*** ¡Atención! ***:  Unable to convert between units of "J" and "N"

En este caso, quantities ha generado una excepción que nos avisa de la incoherencia dimensional de la operación que queremos llevar a cabo.

La utilidad de quantities va mucho más allá, al permitir que sea el intérprete el que lleva el control de las unidades compuestas en que se expresa el resultado de un cálculo y realiza automáticamente las conversiones adecuadas para que el resultado sea correcto. Esto significa que podemos introducir los valores de las magnitudes en las unidades que en cada momento nos resulten más cómodas:

In [4]:

x = 0.3 * pq.km    # longitud en Km
t = 250 * pq.ms     # tiempo en milisegundos
acel = x/t**2
masa = 7 * pq.g  # masa en gramos
fuerza = masa * acel
print fuerza

3.36e-05 g*km/ms**2

Podemos pedir que el resultado se exprese en el sistema de unidades de base (SI por defecto):

In [5]:

print fuerza.simplified

33.6 kg*m/s**2

O bien que se exprese en alguna unidad derivada, siempre que se respete la coherencia dimensional:

In [6]:

print fuerza.rescale(pq.newton)
print fuerza.rescale(pq.dyne)

33.6 N
3360000.0 dyn

Lo importante es que las unidades correctas se irán computando en todos los cálculos intermedios, realizandose las conversiones adecuadas de forma automática, y avisandonos si en algún momento cometemos un error dimensional. Además el resultado final se podrá expresar en las unidades que deseemos, siempre que sean dimensionalmente compatibles con las unidades computadas.

Unidades más utilizadas en astronomía¶

Hemos visto anteriormente algunos ejemplos elementales de utilización de unidades físicas con el paquete quantities. En astronomía su utilización puede ser particularmente útil debido a la diversidad de las unidades utilizadas. A continuación indicaremos algunas de las más habituales. Mostraremos además la expresión "simplificada" para poder comprobar el valor numérico en unidades del SI:

In [7]:

# Unidades de tiempo y distancias
print pq.year # tiempo en años
print pq.au, pq.au.simplified # distancia en unidades astronómicas
print pq.light_year, pq.light_year.simplified # distancia en años luz
print pq.pc, pq.pc.simplified # distancia en parsecs

1 yr (year)
1 au (astronomical_unit) 1.49597870691e+11 m
1 ly (light_year) 9.46073047258e+15 m
1 pc (parsec) 3.08568025e+16 m

In [8]:

# Unidades de temperatura
print pq.K
print pq.celsius

1 K (Kelvin)
1 degC (Celsius)

Pero, mucho cuidado, ya que quantities no realiza conversiones entre escalas de temperaturas absolutas y temperaturas en grados, ya que eso implica tener en cuenta un punto de referencia (el cero absoluto) para lo que no está preparado el paquete. Las temperaturas se consideran siempre diferencias de temperaturas.

In [9]:

# En este ejemplo se vé que la conversión no se lleva a cabo
temp = 345 * pq.celsius
print temp.rescale(pq.K)

345.0 K

In [10]:

# Tampoco es correcta la conversión entre las escalas Farhrenheit y Celsius
# Sin embargo, si se considera el valor como un incremento de temperaturas,
# el resultado es el correcto
temp = 170 * pq.fahrenheit
print temp.rescale(pq.celsius)

94.4444444444 degC

In [11]:

# Medida de ángulos
print pq.deg
print pq.rad
print pq.arcmin
print pq.arcsec

1 deg (arcdegree)
1 rad (radian)
1 arcmin (arcminute)
1 arcsec (arcsecond)

Ejemplos de utilización de unidades angulares:

In [12]:

print (np.pi/2 * pq.rad).rescale(pq.deg)

90.0 deg

In [13]:

print (1 * pq.deg).rescale(pq.arcmin)

60.0 arcmin

Es posible definir nuevas unidades que no existen en el paquete:

In [14]:

# Definición de nuevas unidades
Mly = pq.UnitQuantity('megalight_year', 1e6*pq.ly, symbol = 'Mly')
Mpc = pq.UnitQuantity('megaparsec', 1e6*pq.pc, symbol = 'Mpc')

Por ejemplo, ¿cuantos años luz hay en un megaparsec?

In [15]:

print (1 * Mpc).rescale(pq.ly)

3261566.59778 ly

En el ejemplo anterior puede verse que la nueva unidad definida no forma parte del espacio de nombres de quantities, es decir, escribimos Mpc, no pq.Mpc

Constantes¶

El paquete quantities también posee un amplio repertorio de constantes físicas. A título de ejmplo citaremos las siguientes, de interés en astronomía:

In [16]:

print pq.c, pq.c.simplified # Velocidad de la luz
print pq.constants.G, pq.constants.G.simplified # Constante newtoniana de gravitación

1 c (speed_of_light) 299792458.0 m/s
1 G (Newtonian_constant_of_gravitation) 6.67428e-11 m**3/(kg*s**2)

Y se pueden definir nuevas constantes. Personalmente utilizo las siguientes:

In [17]:

# Definición de nuevas constantes:
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)

Soporte de NumPy¶

Una gran cualidad de quantities es que es totalmente compatible con los arrays de NumPy. Puede definirse el tipo de unidades de un array:

In [18]:

xa = np.arange(10) * pq.N

Cada elemento del array anterior tendrá unidades en Newton

In [19]:

print xa[5]

5.0 N

Y se soportan operaciones con arrays: en el ejemplo siguiente, al dividir el array anterior por una megnitud de masa, todos los elementos del array tendrán unidades de aceleración:

In [20]:

print (xa / (3.0 * pq.kg)).simplified

[ 0.          0.33333333  0.66666667  1.          1.33333333  1.66666667
  2.          2.33333333  2.66666667  3.        ] m/s**2

En ocasiones es necesario obtener el valor numérico de la variable sin incluir el tipo de unidad. Esto se puede hacer de dos maneras, según queramos obtener un valor numérico o un array de NumPy:

In [21]:

# Un solo valor numérico
f = 30 * pq.N
print f.item(), type(f.item())

# Un array de valores
print xa.magnitude, type(f.magnitude)

30.0 <type 'float'>
[ 0.  1.  2.  3.  4.  5.  6.  7.  8.  9.] <type 'numpy.ndarray'>

Explorando el contenido del paquete quantities¶

El paquete quantities.units incluye un gran número de unidades organizadas en una serie de submódulos, uno para cada tipo de unidades. Para listar los diferentes submódulos se puede escribir lo siguiente:

In [22]:

for nombre, datos in inspect.getmembers(pq.units, inspect.ismodule):
    print nombre, ',',

acceleration , angle , area , compound , electromagnetism , energy , force , frequency , heat , information , length , mass , power , prefixes , pressure , radiation , substance , temperature , time , velocity , viscosity , volume ,

A continuación, se puede obtener una lista de los nombres de las unidades en cada paquete. Por ejemplo, para ver cuales son las unidades de tipo 'mass':

In [23]:

for u in dir(pq.mass):
    print u,',',

Da , UnitMass , __builtins__ , __doc__ , __file__ , __name__ , __package__ , absolute_import , amu , apdram , apothecary_ounce , apothecary_pound , apounce , appound , atomic_mass_unit , avoirdupois_ounce , avoirdupois_pound , bag , carat , dalton , denier , dr , drachm , dram , dtex , dwt , g , gr , grain , gram , kg , kilogram , lb , long_hundredweight , long_ton , metric_ton , mg , milligram , ounce , oz , pennyweight , pound , scruple , short_hundredweight , short_ton , slug , slugs , st , stone , t , tex , ton , tonne , toz , troy_ounce , troy_pound , u ,

El otro package importante es el de constantes:

In [24]:

for nombre, datos in inspect.getmembers(pq.constants, inspect.ismodule):
    print nombre, ',',

_codata , _utils , astronomy , atomicunits , deuteron , electromagnetism , electron , helion , mathematical , muon , naturalunits , neutron , proton , quantum , relationships , statisticalmechanics , tau , triton , weak , xray ,

Y para ver por ejemplo el contenido del módulo de constantes astronómicas:

In [25]:

dir(pq.constants.astronomy)

Out[25]:

['G',
 'Newtonian_constant_of_gravitation',
 '__builtins__',
 '__doc__',
 '__file__',
 '__name__',
 '__package__',
 'absolute_import',
 'astronomical_unit',
 'au']