Commit 65b9d0b6 authored by Peter Eastman's avatar Peter Eastman
Browse files

Python API wrappers are now part of OpenMM

parent abb19052
"""
Physical quantities with units for dimensional analysis and automatic unit conversion.
"""
__docformat__ = "epytext en"
__author__ = "Christopher M. Bruns"
__copyright__ = "Copyright 2010, Stanford University and Christopher M. Bruns"
__credits__ = []
__license__ = "MIT"
__maintainer__ = "Christopher M. Bruns"
__email__ = "cmbruns@stanford.edu"
from unit import Unit, is_unit
from quantity import Quantity, is_quantity
from unit_math import *
from unit_definitions import *
from constants import *
#!/bin/env python
"""
Module simtk.unit.basedimension
BaseDimension class for use by units and quantities.
BaseDimensions are things like "length" and "mass".
"""
__author__ = "Christopher M. Bruns"
__version__ = "0.6"
class BaseDimension(object):
'''
A physical dimension such as length, mass, or temperature.
It is unlikely the user will need to create new ones.
'''
# Keep deterministic order of dimensions
_index_by_name = {
'mass': 1,
'length': 2,
'time': 3,
'temperature': 4,
'amount': 5,
'charge': 6,
'luminous intensity': 7,
'angle': 8,
}
_next_unused_index = 9
def __init__(self, name):
"""Create a new BaseDimension.
Each new BaseDimension is assumed to be independent of all other BaseDimensions.
Use the existing BaseDimensions in simtk.dimension instead of creating
new ones.
"""
self.name = name
if not self.name in BaseDimension._index_by_name.keys():
BaseDimension._index_by_name[name] = BaseDimension._next_unused_index
BaseDimension._next_unused_index += 1
def __cmp__(self, other):
"""
The implicit order of BaseDimensions is the order in which they were created.
This method is used for using BaseDimensions as hash keys, and also affects
the order in which units appear in multi-dimensional Quantities.
Returns 0 if self == other, -1 if self < other, and 1 if self > other.
"""
return cmp(BaseDimension._index_by_name[self.name], BaseDimension._index_by_name[other.name])
def __hash__(self):
"""
Needed for using BaseDimensions as hash keys.
"""
return hash(BaseDimension._index_by_name[self.name])
def __repr__(self):
return 'BaseDimension("%s")' % self.name
# run module directly for testing
if __name__=='__main__':
# Test the examples in the docstrings
import doctest, sys
doctest.testmod(sys.modules[__name__])
#!/bin/env python
"""
Module simtk.unit.baseunit
Contains BaseUnit class, which is a component of the
Unit class.
"""
__author__ = "Christopher M. Bruns"
__version__ = "0.6"
class BaseUnit(object):
'''
Physical unit expressed in exactly one BaseDimension.
For example, meter_base_unit could be a BaseUnit for the length dimension.
The BaseUnit class is used internally in the more general Unit class.
'''
def __init__(self, base_dim, name, symbol):
"""Creates a new BaseUnit.
Parameters
- self: The newly created BaseUnit.
- base_dim: (BaseDimension) The dimension of the new unit, e.g. 'mass'
- name: (string) Name of the unit, e.g. "kilogram"
- symbol: (string) Symbol for the unit, e.g. 'kg'. This symobol will
Quantity string descriptions.
"""
self.dimension = base_dim
self.name = name
self.symbol = symbol
self._conversion_factor_to = {}
self._conversion_factor_to[self] = 1.0
self._conversion_factor_to_by_name = {}
self._conversion_factor_to_by_name[self.name] = 1.0
def __cmp__(self, other):
"""
Comparison function that sorts BaseUnits by BaseDimension
"""
# First sort on dimension
c = cmp(self.dimension, other.dimension)
if c != 0: return c
# Second on conversion factor
return cmp(self.conversion_factor_to(other), 1.0)
def iter_base_dimensions(self):
"""
Returns a dictionary of BaseDimension:exponent pairs, describing the dimension of this unit.
"""
yield (self.dimension, 1)
def iter_base_units(self):
yield (self, 1)
def get_dimension_tuple(self):
"""
Returns a sorted tuple of (BaseDimension, exponent) pairs, that can be used as a dictionary key.
"""
l = list(self.iter_base_dimensions())
l.sort()
return tuple(l)
def __str__(self):
"""Returns a string with the name of this BaseUnit
"""
return self.name
def __repr__(self):
return 'BaseUnit(base_dim=%s, name="%s", symbol="%s")' % (self.dimension, self.name, self.symbol)
def define_conversion_factor_to(self, other, factor):
"""
Defines a conversion factor between two BaseUnits.
self * factor = other
Parameters:
- self: (BaseUnit) 'From' unit in conversion.
- other: (BaseUnit) 'To' unit in conversion.
- factor: (float) Conversion factor.
After calling this method, both self and other will have stored
conversion factors for one another, plus all other BaseUnits which
self and other have previously defined.
Both self and other must have the same dimension, otherwise a TypeError
will be raised.
Returns None.
"""
if self.dimension != other.dimension:
raise TypeError('Cannot define conversion for BaseUnits with different dimensions.')
assert(factor != 0)
assert(not self is other)
# import all transitive conversions
self._conversion_factor_to[other] = factor
self._conversion_factor_to_by_name[other.name] = factor
for (unit, cfac) in other._conversion_factor_to.items():
if unit is self: continue
if self._conversion_factor_to.has_key(unit): continue
self._conversion_factor_to[unit] = factor * cfac
unit._conversion_factor_to[self] = pow(factor * cfac, -1)
self._conversion_factor_to_by_name[unit.name] = factor * cfac
unit._conversion_factor_to_by_name[self.name] = pow(factor * cfac, -1)
# and for the other guy
invFac = pow(factor, -1.0)
other._conversion_factor_to[self] = invFac
other._conversion_factor_to_by_name[self.name] = invFac
for (unit, cfac) in self._conversion_factor_to.items():
if unit is other: continue
if other._conversion_factor_to.has_key(unit): continue
other._conversion_factor_to[unit] = invFac * cfac
unit._conversion_factor_to[other] = pow(invFac * cfac, -1)
other._conversion_factor_to_by_name[unit.name] = invFac * cfac
unit._conversion_factor_to_by_name[other.name] = pow(invFac * cfac, -1)
def conversion_factor_to(self, other):
"""Returns a conversion factor from this BaseUnit to another BaseUnit.
It does not matter which existing BaseUnit you define the conversion factor to.
Conversions for all other known BaseUnits will be computed at the same time.
Raises TypeError if dimension does not match.
Raises LookupError if no conversion has been defined. (see define_conversion_factor_to).
"""
if self is other: return 1.0
if self.dimension != other.dimension:
raise TypeError('Cannot get conversion for BaseUnits with different dimensions.')
if not other.name in self._conversion_factor_to_by_name:
raise LookupError('No conversion defined from BaseUnit "%s" to "%s".' % (self, other))
return self._conversion_factor_to_by_name[other.name]
# run module directly for testing
if __name__=='__main__':
# Test the examples in the docstrings
import doctest, sys
doctest.testmod(sys.modules[__name__])
#!/bin/env python
"""
Module simtk.unit.constants
"""
__author__ = "Christopher M. Bruns"
__version__ = "0.5"
from unit_definitions import *
#################
### CONSTANTS ###
#################
# codata 2006
AVOGADRO_CONSTANT_NA = 6.02214179e23 / mole
BOLTZMANN_CONSTANT_kB = 1.3806504e-23 * joule / kelvin
MOLAR_GAS_CONSTANT_R = AVOGADRO_CONSTANT_NA * BOLTZMANN_CONSTANT_kB
# From simtkcommon
SPEED_OF_LIGHT_C = 2.99792458e8 * meter / second
GRAVITATIONAL_CONSTANT_G = 6.6742e-11 * newton * meter**2 / kilogram**2
# run module directly for testing
if __name__=='__main__':
# Test the examples in the docstrings
import doctest, sys
doctest.testmod(sys.modules[__name__])
This diff is collapsed.
"""
Pure python inversion of small matrices, to avoid requiring numpy or similar in SimTK.
"""
import sys
def eye(size):
"""
Returns identity matrix.
>>> print eye(3)
[[1, 0, 0]
[0, 1, 0]
[0, 0, 1]]
"""
result = []
for row in range(0, size):
r = []
for col in range(0, size):
if row == col:
r.append(1)
else:
r.append(0)
result.append(r)
return MyMatrix(result)
def zeros(m, n=None):
"""
Returns matrix of zeroes
>>> print zeros(3)
[[0, 0, 0]
[0, 0, 0]
[0, 0, 0]]
"""
if n == None:
n = m
result = []
for row in range(0, m):
r = []
for col in range(0, n):
r.append(0)
result.append(r)
return MyMatrix(result)
class MyVector(object):
"""
Parent class of MyMatrix and type of Matrix Row.
"""
def __init__(self, collection):
if isinstance(collection, MyVector):
self.data = collection.data
else:
self.data = collection
def __str__(self):
return str(self.data)
def __repr__(self):
return self.__class__.__name__ + "(" + repr(self.data) + ")"
def __getitem__(self, key):
return self.data[key]
def __contains__(self, item):
return item in self.data
def __delitem__(self, key):
del self.data[key]
def __iter__(self):
for item in self.data:
yield item
def __len__(self):
return len(self.data)
def __setitem__(self, key, value):
self.data[key] = value
def __rmul__(self, lhs):
try:
len(lhs)
# left side is not scalar, delegate mul to that class
return NotImplemented
except TypeError:
new_vec = []
for element in self:
new_vec.append(lhs * element)
return self.__class__(new_vec)
class MyMatrix(MyVector):
"""
Pure python linear algebra matrix for internal matrix inversion in UnitSystem.
>>> m = MyMatrix([[1,0,],[0,1,]])
>>> print m
[[1, 0]
[0, 1]]
>>> print ~m
[[1.0, 0.0]
[0.0, 1.0]]
>>> print eye(5)
[[1, 0, 0, 0, 0]
[0, 1, 0, 0, 0]
[0, 0, 1, 0, 0]
[0, 0, 0, 1, 0]
[0, 0, 0, 0, 1]]
>>> m = eye(5)
>>> m[1][1]
1
>>> m[1:4]
MyMatrixTranspose([[0, 0, 0],[1, 0, 0],[0, 1, 0],[0, 0, 1],[0, 0, 0]])
>>> print m[1:4]
[[0, 0, 0]
[1, 0, 0]
[0, 1, 0]
[0, 0, 1]
[0, 0, 0]]
>>> print m[1:4][0:2]
[[0, 1]
[0, 0]
[0, 0]]
>>> m[1:4][0:2] = [[9,8],[7,6],[5,4]]
>>> print m
[[1, 0, 0, 0, 0]
[9, 8, 0, 0, 0]
[7, 6, 1, 0, 0]
[5, 4, 0, 1, 0]
[0, 0, 0, 0, 1]]
"""
def numRows(self):
return len(self.data)
def numCols(self):
if len(self.data) == 0:
return 0
else:
return len(self.data[0])
def __len__(self):
return self.numRows()
def __str__(self):
result = ""
start_char = "["
for m in range(0, self.numRows()):
result += start_char
result += str(self[m])
if m < self.numRows() - 1:
result += "\n"
start_char = " "
result += "]"
return result
def __repr__(self):
return 'MyMatrix(' + MyVector.__repr__(self) + ')'
def is_square(self):
return self.numRows() == self.numCols()
def __iter__(self):
for item in self.data:
yield MyVector(item)
def __getitem__(self, m):
if isinstance(m, slice):
return MyMatrixTranspose(self.data[m])
else:
return MyVector(self.data[m])
def __setitem__(self, key, rhs):
if isinstance(key, slice):
self.data[key] = rhs
else:
assert len(rhs) == self.numCols()
self.data[key] = MyVector(rhs)
def __mul__(self, rhs):
"""
Matrix multiplication.
>>> a = MyMatrix([[1,2],[3,4]])
>>> b = MyMatrix([[5,6],[7,8]])
>>> print a
[[1, 2]
[3, 4]]
>>> print b
[[5, 6]
[7, 8]]
>>> print a*b
[[19, 22]
[43, 50]]
"""
m = self.numRows()
n = len(rhs[0])
r = len(rhs)
if self.numCols() != r:
raise ArithmeticError("Matrix multplication size mismatch (%d vs %d)" % (self.numCols(), r))
result = zeros(m, n)
for i in range(0, m):
for j in range(0, n):
for k in range(0, r):
result[i][j] += self[i][k]*rhs[k][j]
return result
def __add__(self, rhs):
"""
Matrix addition.
>>> print MyMatrix([[1, 2],[3, 4]]) + MyMatrix([[5, 6],[7, 8]])
[[6, 8]
[10, 12]]
"""
m = self.numRows()
n = self.numCols()
assert len(rhs) == m
assert len(rhs[0]) == n
result = zeros(m,n)
for i in range(0,m):
for j in range(0,n):
result[i][j] = self[i][j] + rhs[i][j]
return result
def __sub__(self, rhs):
"""
Matrix subtraction.
>>> print MyMatrix([[1, 2],[3, 4]]) - MyMatrix([[5, 6],[7, 8]])
[[-4, -4]
[-4, -4]]
"""
m = self.numRows()
n = self.numCols()
assert len(rhs) == m
assert len(rhs[0]) == n
result = zeros(m,n)
for i in range(0,m):
for j in range(0,n):
result[i][j] = self[i][j] - rhs[i][j]
return result
def __pos__(self):
return self
def __neg__(self):
m = self.numRows()
n = self.numCols()
result = zeros(m, n)
for i in range(0,m):
for j in range(0,n):
result[i][j] = -self[i][j]
return result
def __invert__(self):
"""
>>> m = MyMatrix([[1,1],[0,1]])
>>> print m
[[1, 1]
[0, 1]]
>>> print ~m
[[1.0, -1.0]
[0.0, 1.0]]
>>> print m*~m
[[1.0, 0.0]
[0.0, 1.0]]
>>> print ~m*m
[[1.0, 0.0]
[0.0, 1.0]]
>>> m = MyMatrix([[1,0,0],[0,0,1],[0,-1,0]])
>>> print m
[[1, 0, 0]
[0, 0, 1]
[0, -1, 0]]
>>> print ~m
[[1.0, 0.0, 0.0]
[0.0, 0.0, -1.0]
[0.0, 1.0, 0.0]]
>>> print m*~m
[[1.0, 0.0, 0.0]
[0.0, 1.0, 0.0]
[0.0, 0.0, 1.0]]
>>> print ~m*m
[[1.0, 0.0, 0.0]
[0.0, 1.0, 0.0]
[0.0, 0.0, 1.0]]
"""
assert self.is_square()
if self.numRows() == 0:
return self
elif self.numRows() == 1:
val = self[0][0]
val = 1.0/val
return MyMatrix([[val]])
elif self.numRows() == 2: # 2x2 is analytic
# http://en.wikipedia.org/wiki/Invertible_matrix#Inversion_of_2.C3.972_matrices
a = self[0][0]
b = self[0][1]
c = self[1][0]
d = self[1][1]
determinant = a*d - b*c
if determinant == 0:
raise ArithmeticError("Cannot invert 2x2 matrix with zero determinant")
else:
return 1.0/(a*d - b*c) * MyMatrix([[d, -b],[-c, a]])
else:
# Gauss Jordan elimination from numerical recipes
n = self.numRows()
m1 = self.numCols()
assert n == m1
# Copy initial matrix into result matrix
a = zeros(n, n)
for i in range (0,n):
for j in range (0,n):
a[i][j] = self[i][j]
# These arrays are used for bookkeeping on the pivoting
indxc = [0] * n
indxr = [0] * n
ipiv = [0] * n
for i in range (0,n):
big = 0.0
for j in range (0,n):
if ipiv[j] != 1:
for k in range (0,n):
if ipiv[k] == 0:
if abs(a[j][k]) >= big:
big = abs(a[j][k])
irow = j
icol = k
ipiv[icol] += 1
# We now have the pivot element, so we interchange rows...
if irow != icol:
for l in range(0,n):
temp = a[irow][l]
a[irow][l] = a[icol][l]
a[icol][l] = temp
indxr[i] = irow
indxc[i] = icol
if a[icol][icol] == 0:
raise ArithmeticError("Cannot invert singular matrix")
pivinv = 1.0/a[icol][icol]
a[icol][icol] = 1.0
for l in range(0,n):
a[icol][l] *= pivinv
for ll in range(0,n): # next we reduce the rows
if ll == icol:
continue # except the pivot one, of course
dum = a[ll][icol]
a[ll][icol] = 0.0
for l in range(0,n):
a[ll][l] -= a[icol][l]*dum
# Unscramble the permuted columns
for l in range(n-1, -1, -1):
if indxr[l] == indxc[l]:
continue
for k in range(0,n):
temp = a[k][indxr[l]]
a[k][indxr[l]] = a[k][indxc[l]]
a[k][indxc[l]] = temp
return a
def transpose(self):
return MyMatrixTranspose(self.data)
class MyMatrixTranspose(MyMatrix):
def transpose(self):
return MyMatrix(self.data)
def numRows(self):
if len(self.data) == 0:
return 0
else:
return len(self.data[0])
def numCols(self):
return len(self.data)
def __getitem__(self, key):
result = []
for row in self.data:
result.append(row[key])
if isinstance(key, slice):
return MyMatrix(result)
else:
return MyVector(result)
def __setitem__(self, key, rhs):
for n in range(0, len(self.data)):
self.data[n][key] = rhs[n]
def __str__(self):
if len(self.data) == 0:
return "[[]]"
start_char = "["
result = ""
for m in range(0, len(self.data[0])):
result += start_char
result += "["
sep_char = ""
for n in range(0, len(self.data)):
result += sep_char
result += str(self.data[n][m])
sep_char = ", "
result += "]"
if m < len(self.data[0]) - 1:
result += "\n"
start_char = " "
result += "]"
return result
def __repr__(self):
if len(self.data) == 0:
return "MyMatrixTranspose([[]])"
start_char = "["
result = 'MyMatrixTranspose('
for m in range(0, len(self.data[0])):
result += start_char
result += "["
sep_char = ""
for n in range(0, len(self.data)):
result += sep_char
result += repr(self.data[n][m])
sep_char = ", "
result += "]"
if m < len(self.data[0]) - 1:
result += ","
start_char = ""
result += '])'
return result
# run module directly for testing
if __name__=='__main__':
# Test the examples in the docstrings
import doctest, sys
doctest.testmod(sys.modules[__name__])
#!/bin/env python
"""
Module simtk.unit.prefix
"""
__author__ = "Christopher M. Bruns"
__version__ = "0.6"
from baseunit import BaseUnit
from unit import Unit, ScaledUnit
import sys
###################
### SI PREFIXES ###
###################
class SiPrefix(object):
"""
Unit prefix that can be multiplied by a unit to yield a new unit.
e.g. millimeter = milli*meter
"""
def __init__(self, prefix, factor, symbol):
self.prefix = prefix
self.factor = factor
self.symbol = symbol
def __mul__(self, unit):
"""
SiPrefix * BaseUnit yields new BaseUnit
SiPrefix * ScaledUnit yields new ScaledUnit
SiPrefix * Unit with exactly one BaseUnit or ScaledUnit yields new Unit
"""
if isinstance(unit, BaseUnit):
# BaseUnit version
symbol = self.symbol + unit.symbol
name = self.prefix + unit.name
factor = self.factor
# TODO - check for existing BaseUnit with same name, symbol, and factor
new_base_unit = BaseUnit(unit.dimension, name, symbol)
new_base_unit.define_conversion_factor_to(unit, factor)
return new_base_unit
elif isinstance(unit, ScaledUnit):
# ScaledUnit version
symbol = self.symbol + unit.symbol
name = self.prefix + unit.name
factor = self.factor * unit.factor
# TODO - check for existing BaseUnit with same name, symbol, and factor
return ScaledUnit(factor, unit.master, name, symbol)
elif isinstance(unit, Unit):
base_units = list(unit.iter_base_or_scaled_units())
if 1 != len(base_units):
raise TypeError('Unit prefix "%s" can only be with simple Units containing one component.' % self.prefix)
if 1 != base_units[0][1]:
raise TypeError('Unit prefix "%s" can only be with simple Units with an exponent of 1.' % self.prefix)
base_unit = base_units[0][0]
# Delegate to Base or Scaled Unit multiply
new_base_unit = self * base_unit
new_unit = Unit({new_base_unit: 1.0})
return new_unit
else:
raise TypeError('Unit prefix "%s" can only be applied to a Unit, BaseUnit, or ScaledUnit.' % self.prefix)
yotto = SiPrefix("yotto", 1e-24, 'y')
zepto = SiPrefix("zepto", 1e-21, 'z')
atto = SiPrefix("atto", 1e-18, 'a')
femto = SiPrefix("femto", 1e-15, 'f')
pico = SiPrefix("pico", 1e-12, 'p')
nano = SiPrefix("nano", 1e-9, 'n')
micro = SiPrefix("micro", 1e-6, 'u')
milli = SiPrefix("milli", 1e-3, 'm')
centi = SiPrefix("centi", 1e-2, 'c')
deci = SiPrefix("deci", 1e-1, 'd')
# two spellings for deka prefix
deka = SiPrefix("deka", 1e1, 'da')
deca = SiPrefix("deca", 1e1, 'da')
hecto = SiPrefix("hecto", 1e2, 'h')
kilo = SiPrefix("kilo", 1e3, 'k')
mega = SiPrefix("mega", 1e6, 'M')
giga = SiPrefix("giga", 1e9, 'G')
tera = SiPrefix("tera", 1e12, 'T')
peta = SiPrefix("peta", 1e15, 'P')
exa = SiPrefix("exa", 1e18, 'E')
zetta = SiPrefix("zetta", 1e21, 'Z')
yotta = SiPrefix("yotta", 1e24, 'Y')
si_prefixes = ( yotto
, zepto
, atto
, femto
, pico
, nano
, micro
, milli
, centi
, deci
, deka
, deca
, hecto
, kilo
, mega
, giga
, tera
, peta
, exa
, zetta
, yotta)
def define_prefixed_units(base_unit, module = sys.modules[__name__]):
"""
Create attributes for prefixed units derived from a particular BaseUnit, e.g. "kilometer" from "meter_base_unit"
Parameters
- base_unit (BaseUnit) existing base unit to use as a basis for prefixed units
- module (Module) module which will contain the new attributes. Defaults to simtk.unit module.
"""
for prefix in si_prefixes:
new_base_unit = prefix * base_unit
name = new_base_unit.name
new_unit = Unit({new_base_unit: 1.0})
# Create base_unit attribute, needed for creating UnitSystems
module.__dict__[name + '_base_unit'] = new_base_unit # e.g. "kilometer_base_unit"
# Create attribue in this module
module.__dict__[name] = new_unit # e.g. "kilometer"
# And plural version
module.__dict__[name + 's'] = new_unit # e.g. "kilometers"
# Binary prefixes
kibi = SiPrefix("kibi", 2.0**10, 'Ki')
mebi = SiPrefix("mebi", 2.0**20, 'Mi')
gibi = SiPrefix("gibi", 2.0**30, 'Gi')
tebi = SiPrefix("tebi", 2.0**40, 'Ti')
pebi = SiPrefix("pebi", 2.0**50, 'Pi')
exbi = SiPrefix("exbi", 2.0**60, 'Ei')
zebi = SiPrefix("zebi", 2.0**70, 'Zi')
yobi = SiPrefix("yobi", 2.0**80, 'Yi')
binary_prefixes = ( kibi
, mebi
, gibi
, tebi
, pebi
, exbi
, zebi
, yobi)
# run module directly for testing
if __name__=='__main__':
# Test the examples in the docstrings
import doctest, sys
doctest.testmod(sys.modules[__name__])
This diff is collapsed.
#!/bin/env python
"""
Module simtk.unit.standard_dimensions
Definition of principal dimensions: mass, length, time, etc.
"""
__author__ = "Christopher M. Bruns"
__version__ = "0.6"
from basedimension import BaseDimension
##################
### DIMENSIONS ###
##################
mass_dimension = BaseDimension('mass')
length_dimension = BaseDimension('length')
time_dimension = BaseDimension('time')
temperature_dimension = BaseDimension('temperature')
amount_dimension = BaseDimension('amount')
charge_dimension = BaseDimension('charge')
luminous_intensity_dimension = BaseDimension('luminous intensity')
angle_dimension = BaseDimension('angle')
information_dimension = BaseDimension('information')
# run module directly for testing
if __name__=='__main__':
# Test the examples in the docstrings
import doctest, sys
doctest.testmod(sys.modules[__name__])
This diff is collapsed.
#!/bin/env python
"""
Module simtk.unit.unit_definitions
"""
__author__ = "Christopher M. Bruns"
__version__ = "0.6"
from baseunit import BaseUnit
from standard_dimensions import *
from unit import Unit, ScaledUnit, UnitSystem, dimensionless
from unit_operators import * ; # needed for manipulation of units
from prefix import *
import math
import sys
#####################
### DIMENSIONLESS ###
#####################
# dimensionless = Unit({}); # defined in unit.py
##############
### LENGTH ###
##############
meter_base_unit = BaseUnit(length_dimension, "meter", "m")
meters = meter = Unit({meter_base_unit: 1.0})
define_prefixed_units(meter_base_unit, module = sys.modules[__name__])
angstrom_base_unit = BaseUnit(length_dimension, "angstrom", "A")
angstrom_base_unit.define_conversion_factor_to(meter_base_unit, 1e-10)
angstroms = angstrom = Unit({angstrom_base_unit: 1.0})
planck_length_base_unit = BaseUnit(length_dimension, "Planck length", "l_P")
planck_length_base_unit.define_conversion_factor_to(meter_base_unit, 1.61625281e-35)
inch_base_unit = BaseUnit(length_dimension, "inch", "in")
inch_base_unit.define_conversion_factor_to(centimeter_base_unit, 2.5400)
inch = inches = Unit({inch_base_unit: 1.0})
foot_base_unit = BaseUnit(length_dimension, "foot", "ft")
foot_base_unit.define_conversion_factor_to(inch_base_unit, 12.0)
foot = feet = Unit({foot_base_unit: 1.0})
yard_base_unit = BaseUnit(length_dimension, "yard", "yd")
yard_base_unit.define_conversion_factor_to(foot_base_unit, 3.0)
yard = yards = Unit({yard_base_unit: 1.0})
furlongs = furlong = yard.create_unit(scale=220.0, name="furlong", symbol="furlong")
miles = mile = furlong.create_unit(scale=8.0, name="mile", symbol="mi")
############
### MASS ###
############
gram_base_unit = BaseUnit(mass_dimension, "gram", "g")
grams = gram = Unit({gram_base_unit: 1.0})
define_prefixed_units(gram_base_unit, module = sys.modules[__name__])
planck_mass_base_unit = BaseUnit(mass_dimension, "Planck mass", "m_P")
planck_mass_base_unit.define_conversion_factor_to(kilogram_base_unit, 2.1764411e-8)
# pound can be mass, force, or currency
pound_mass_base_unit = BaseUnit(mass_dimension, "pound", "lb")
pound_mass_base_unit.define_conversion_factor_to(kilogram_base_unit, 0.3732)
pound_mass = pounds_mass = Unit({pound_mass_base_unit: 1.0})
stone_base_unit = BaseUnit(mass_dimension, "stone", "stone")
stone_base_unit.define_conversion_factor_to(pound_mass_base_unit, 14.0)
stone = stones = Unit({stone_base_unit: 1.0})
############
### TIME ###
############
second_base_unit = BaseUnit(time_dimension, "second", "s")
seconds = second = Unit({second_base_unit: 1.0})
define_prefixed_units(second_base_unit, module = sys.modules[__name__])
planck_time_base_unit = BaseUnit(time_dimension, "Planck time", "t_P")
planck_time_base_unit.define_conversion_factor_to(second_base_unit, 5.3912427e-44)
minutes = minute = second.create_unit(scale=60.0, name="minute", symbol="min")
hours = hour = minute.create_unit(scale=60.0, name="hour", symbol="hr")
days = day = hour.create_unit(scale=24.0, name="day", symbol="day")
weeks = week = day.create_unit(scale=7.0, name="week", symbol="week")
years = year = day.create_unit(scale=365.25, name="julian year", symbol="a")
centuries = centurys = century = year.create_unit(scale=100.0, name="century", symbol="century")
millenia = milleniums = millenium = century.create_unit(scale=10.0, name="millenium", symbol="ka")
fortnights = fortnight = day.create_unit(scale=14.0, name="fortnight", symbol="fortnight")
###################
### TEMPERATURE ###
###################
kelvin_base_unit = BaseUnit(temperature_dimension, "kelvin", "K")
kelvins = kelvin = Unit({kelvin_base_unit: 1.0})
planck_temperature_base_unit = BaseUnit(temperature_dimension, "Planck temperature", "T_p")
planck_temperature_base_unit.define_conversion_factor_to(kelvin_base_unit, 1.41678571e32)
##############
### CHARGE ###
##############
elementary_charge_base_unit = BaseUnit(charge_dimension, "elementary charge", "q")
elementary_charges = elementary_charge = Unit({elementary_charge_base_unit: 1.0})
coulomb_base_unit = BaseUnit(charge_dimension, "elementary charge", "q")
# Exact conversion factor
coulomb_base_unit.define_conversion_factor_to(elementary_charge_base_unit, 6.24150962915265e18)
coulombs = coulomb = Unit({coulomb_base_unit: 1.0})
planck_charge_base_unit = BaseUnit(charge_dimension, "Planck charge", "q_P")
planck_charge_base_unit.define_conversion_factor_to(elementary_charge_base_unit, 11.706237639840)
##############
### AMOUNT ###
##############
mole_base_unit = BaseUnit(amount_dimension, "mole", "mol")
moles = mole = Unit({mole_base_unit: 1.0})
single_item_amount_base_unit = BaseUnit(amount_dimension, "item", "")
mole_base_unit.define_conversion_factor_to(single_item_amount_base_unit, 6.0221417930e23)
items = item = Unit({single_item_amount_base_unit: 1.0})
##########################
### Luminous Intensity ###
##########################
candela_base_unit = BaseUnit(luminous_intensity_dimension, "candela", "cd")
candelas = candela = Unit({candela_base_unit: 1.0})
#############
### ANGLE ###
#############
radian_base_unit = BaseUnit(angle_dimension, "radian", "rad")
radians = radian = Unit({radian_base_unit: 1.0})
degree_base_unit = BaseUnit(angle_dimension, "degree", "deg")
degree_base_unit.define_conversion_factor_to(radian_base_unit, math.pi/180.0)
degrees = degree = Unit({degree_base_unit: 1.0})
arcminutes = arcminute = degree.create_unit(scale=1.0/60.0, name="arcminute", symbol="'")
arcseconds = arcsecond = arcminute.create_unit(scale=1.0/60.0, name="arcsecond", symbol='"')
###################
### INFORMATION ###
###################
bit_base_unit = BaseUnit(information_dimension, "bit", "bit")
bits = bit = Unit({bit_base_unit: 1.0})
byte_base_unit = BaseUnit(information_dimension, "byte", "B")
byte_base_unit.define_conversion_factor_to(bit_base_unit, 8.0)
bytes = byte = Unit({byte_base_unit: 1.0})
nat_base_unit = BaseUnit(information_dimension, "nat", "nat")
nat_base_unit.define_conversion_factor_to(bit_base_unit, 1.0/math.log(2.0))
nats = nat = nits = nit = nepits = nepit = Unit({nat_base_unit: 1.0})
ban_base_unit = BaseUnit(information_dimension, "ban", "ban")
ban_base_unit.define_conversion_factor_to(bit_base_unit, math.log(10.0, 2.0))
bans = ban = hartleys = hartley = dits = dit = Unit({ban_base_unit: 1.0})
###############
### DERIVED ###
###############
# Molar mass
# daltons = dalton = grams / mole
daltons = dalton = Unit({ScaledUnit(1.0, gram/mole, "dalton", "Da"): 1.0})
amus = amu = dalton
atom_mass_units = atomic_mass_unit = dalton
# Volume
liter_base_unit = ScaledUnit(1.0, decimeter**3, "liter", "l")
liter = liters = litre = litres = Unit({liter_base_unit: 1.0})
define_prefixed_units(liter_base_unit, module = sys.modules[__name__])
# Concentration
molar_base_unit = ScaledUnit(1.0, mole/liter, "molar", "M")
molar = molal = Unit({molar_base_unit: 1.0})
define_prefixed_units(molar_base_unit, module = sys.modules[__name__])
# Force
newton_base_unit = ScaledUnit(1.0, kilogram * meter / second / second, "newton", "N")
newtons = newton = Unit({newton_base_unit: 1.0})
define_prefixed_units(newton_base_unit, module = sys.modules[__name__])
# pound can be mass, force, or currency
pound_force_base_unit = ScaledUnit(4.448, newton, "pound", "lb")
pound_force = pounds_force = Unit({pound_force_base_unit: 1.0})
dyne_base_unit = ScaledUnit(1.0, gram * centimeter / second**2, "dyne", "dyn")
dyne = dynes = Unit({dyne_base_unit: 1.0})
# Energy
joule_base_unit = ScaledUnit(1.0, newton * meter, "joule", "J")
joules = joule = Unit({joule_base_unit: 1.0})
define_prefixed_units(joule_base_unit, module = sys.modules[__name__])
erg_base_unit = ScaledUnit(1.0, dyne * centimeter, "erg", "erg")
erg = ergs = Unit({erg_base_unit: 1.0})
# In molecular simulations, "kilojoules" are in microscopic units
# And you really only want to use kilojoules/mole.
md_kilojoule_raw = gram * nanometer**2 / picosecond**2
md_kilojoules = md_kilojoule = Unit({ScaledUnit(1.0, md_kilojoule_raw, "kilojoule", "kJ"): 1.0})
kilojoules_per_mole = kilojoule_per_mole = md_kilojoule / mole
calorie_base_unit = ScaledUnit(4.184, joule, "calorie", "cal")
calories = calorie = Unit({calorie_base_unit: 1.0})
define_prefixed_units(calorie_base_unit, module = sys.modules[__name__])
md_kilocalories = md_kilocalorie = Unit({ScaledUnit(4.184, md_kilojoule, "kilocalorie", "kcal"): 1.0})
kilocalories_per_mole = kilocalorie_per_mole = md_kilocalorie / mole
# Power
watt_base_unit = ScaledUnit(1.0, joule / second, "watt", "W")
watt = watts = Unit({watt_base_unit: 1.0})
# Current
ampere_base_unit = ScaledUnit(1.0, coulomb / second, "ampere", "A")
ampere = amperes = amp = amps = Unit({ampere_base_unit: 1.0})
# Electrical potential
volt_base_unit = ScaledUnit(1.0, watt / ampere, "volt", "V")
volt = volts = Unit({volt_base_unit: 1.0})
# Magnetic field
tesla_base_unit = ScaledUnit(1.0, newton / (ampere * meter), "tesla", "T")
tesla = teslas = Unit({tesla_base_unit: 1.0})
gauss_base_unit = ScaledUnit(10.0**-4, tesla, "gauss", "G")
gauss = Unit({gauss_base_unit: 1.0})
# Electrical resistance
ohm_base_unit = ScaledUnit(1.0, volt / ampere, "ohm", "O")
ohm = ohms = Unit({ohm_base_unit: 1.0})
# Capacitance
farad_base_unit = ScaledUnit(1.0, coulomb / volt, "farad", "F")
farad = farads = Unit({farad_base_unit: 1.0})
# Inductance
henry_base_unit = ScaledUnit(1.0, volt * second / ampere, "henry", "H")
henry = henries = henrys = Unit({henry_base_unit: 1.0})
# Pressure
pascal_base_unit = ScaledUnit(1.0, newton / (meter**2), "pascal", "Pa")
pascals = pascal = Unit({pascal_base_unit: 1.0})
define_prefixed_units(pascal_base_unit, module = sys.modules[__name__])
psi_base_unit = ScaledUnit(1.0, pound_force / (inch**2), "psi", "psi")
psi = Unit({psi_base_unit: 1.0})
bar_base_unit = ScaledUnit(10.0**5, pascal, "bar", "bar")
bar = bars = Unit({bar_base_unit: 1.0})
atmosphere_base_unit = ScaledUnit(101325.0, pascal, "atmosphere", "atm")
atmosphere = atmospheres = Unit({atmosphere_base_unit: 1.0})
torr_base_unit = ScaledUnit(1.0/760.0, atmosphere, "torr", "Torr")
torr = Unit({torr_base_unit: 1.0})
mmHg_base_unit = ScaledUnit(133.322, pascal, "mmHg", "mmHg")
mmHg = Unit({mmHg_base_unit: 1.0})
####################
### Unit Systems ###
####################
ampere_base_unit = ScaledUnit(1.0, coulomb/second, "ampere", "A")
si_unit_system = UnitSystem([\
meter_base_unit,\
kilogram_base_unit,\
second_base_unit,\
ampere_base_unit,
kelvin_base_unit,
mole_base_unit,
candela_base_unit,
radian_base_unit])
cgs_unit_system = UnitSystem([\
centimeter_base_unit,\
gram_base_unit,\
second_base_unit,\
ampere_base_unit,
kelvin_base_unit,
mole_base_unit,
radian_base_unit])
dalton_base_unit = ScaledUnit(1.0, gram/mole, "dalton", "Da")
md_unit_system = UnitSystem([\
nanometer_base_unit,\
dalton_base_unit,
picosecond_base_unit,\
elementary_charge_base_unit,
kelvin_base_unit,
mole_base_unit,
radian_base_unit])
planck_unit_system = UnitSystem([\
planck_length_base_unit,
planck_mass_base_unit,
planck_time_base_unit,
planck_charge_base_unit,
planck_temperature_base_unit,
single_item_amount_base_unit,
radian_base_unit])
########################
### TESTING/EXAMPLES ###
########################
# run module directly for testing
if __name__=='__main__':
# Test the examples in the docstrings
import doctest, sys
doctest.testmod(sys.modules[__name__])
#!/bin/env python
"""
Module simtk.unit.math
Arithmetic methods on Quantities and Units
"""
__author__ = "Christopher M. Bruns"
__version__ = "0.5"
import math
from quantity import is_quantity
from unit_definitions import *
####################
### TRIGONOMETRY ###
####################
def sin(angle):
"""
Examples
>>> sin(90*degrees)
1.0
"""
if is_quantity(angle):
return math.sin(angle/radians)
else:
return math.sin(angle)
def sinh(angle):
if is_quantity(angle):
return math.sinh(angle/radians)
else:
return math.sinh(angle)
def cos(angle):
"""
Examples
>>> cos(180*degrees)
-1.0
"""
if is_quantity(angle):
return math.cos(angle/radians)
else:
return math.cos(angle)
def cosh(angle):
if is_quantity(angle):
return math.cosh(angle/radians)
else:
return math.cosh(angle)
def tan(angle):
if is_quantity(angle):
return math.tan(angle/radians)
else:
return math.tan(angle)
def tanh(angle):
if is_quantity(angle):
return math.tanh(angle/radians)
else:
return math.tanh(angle)
def acos(x):
"""
>>> acos(1.0)
Quantity(value=0.0, unit=radian)
>>> print acos(1.0)
0.0 rad
"""
return math.acos(x) * radians
def acosh(x):
return math.acosh(x) * radians
def asin(x):
return math.asin(x) * radians
def asinh(x):
return math.asinh(x) * radians
def atan(x):
return math.atan(x) * radians
def atanh(x):
return math.atanh(x) * radians
def atan2(x, y):
return math.atan2(x, y) * radians
###################
### SQUARE ROOT ###
###################
def sqrt(val):
"""
>>> sqrt(9.0)
3.0
>>> print sqrt(meter*meter)
meter
>>> sqrt(9.0*meter*meter)
Quantity(value=3.0, unit=meter)
>>> sqrt(9.0*meter*meter*meter)
Traceback (most recent call last):
...
ArithmeticError: Exponents in Unit.sqrt() must be even.
"""
try:
return val.sqrt()
except AttributeError:
return math.sqrt(val)
# run module directly for testing
if __name__=='__main__':
# Test the examples in the docstrings
import doctest, sys
doctest.testmod(sys.modules[__name__])
#!/bin/env python
"""
Module simtk.unit.unit_operators
Physical quantities with units, intended to produce similar functionality
to Boost.Units package in C++ (but with a runtime cost).
Uses similar API as Scientific.Physics.PhysicalQuantities
but different internals to satisfy our local requirements.
In particular, there is no underlying set of 'canonical' base
units, whereas in Scientific.Physics.PhysicalQuantities all
units are secretly in terms of SI units. Also, it is easier
to add new fundamental dimensions to simtk.dimensions. You
might want to make new dimensions for, say, "currency" or
"information".
Two possible enhancements that have not been implemented are
1) Include uncertainties with propagation of errors
2) Incorporate offsets for celsius <-> kelvin conversion
"""
__author__ = "Christopher M. Bruns"
__version__ = "0.5"
from unit import Unit, is_unit
from quantity import Quantity, is_quantity
# Attach methods of Unit class that return a Quantity to Unit class.
# I put them here to avoid circular dependence in imports.
# i.e. Quantity depends on Unit, but not vice versa
def _unit_class_rdiv(self, other):
"""
Divide another object type by a Unit.
Returns a new Quantity with a value of other and units
of the inverse of self.
"""
if is_unit(other):
raise NotImplementedError('programmer is surprised __rdiv__ was called instead of __div__')
else:
# print "R scalar / unit"
unit = pow(self, -1.0)
value = other
return Quantity(value, unit).reduce_unit(self)
Unit.__rdiv__ = _unit_class_rdiv
def _unit_class_mul(self, other):
"""Multiply a Unit by an object.
If other is another Unit, returns a new composite Unit.
Exponents of similar dimensions are added. If self and
other share similar BaseDimension, but
with different BaseUnits, the resulting BaseUnit for that
BaseDimension will be that used in self.
If other is a not another Unit, this method returns a
new Quantity... UNLESS other is a Quantity and the resulting
unit is dimensionless, in which case the underlying value type
of the Quantity is returned.
"""
if is_unit(other):
# print "unit * unit"
result1 = {} # dictionary of dimensionTuple: (BaseOrScaledUnit, exponent)
for unit, exponent in self.iter_base_or_scaled_units():
d = unit.get_dimension_tuple()
if d not in result1:
result1[d] = {}
assert unit not in result1[d]
result1[d][unit] = exponent
for unit, exponent in other.iter_base_or_scaled_units():
d = unit.get_dimension_tuple()
if d not in result1:
result1[d] = {}
if unit not in result1[d]:
result1[d][unit] = 0
result1[d][unit] += exponent
result2 = {} # stripped of zero exponents
for d in result1:
for unit in result1[d]:
exponent = result1[d][unit]
if exponent != 0:
assert unit not in result2
result2[unit] = exponent
return Unit(result2)
elif is_quantity(other):
# print "unit * quantity"
value = other._value
unit = self * other.unit
return Quantity(value, unit).reduce_unit(self)
else:
# print "scalar * unit"
value = other
unit = self
return Quantity(other, self).reduce_unit(self)
Unit.__mul__ = _unit_class_mul
Unit.__rmul__ = Unit.__mul__
# run module directly for testing
if __name__=='__main__':
# Test the examples in the docstrings
import doctest, sys
doctest.testmod(sys.modules[__name__])
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment