Merge branch 'feature/_offset_units' of git://github.com/dalito/pint into dalito-feature/_offset_units

12 files changed, 1388 insertions, 174 deletions

diff --git a/AUTHORS b/AUTHORS
index 4300c6c..f10ab38 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -5,9 +5,9 @@ Other contributors, listed alphabetically, are:
 * Alexander Böhn <fish2000@gmail.com>
 * Brend Wanders <b.wanders@utwente.nl>
 * choloepus
-* coutinho <coutinho@esrf.fr>
 * Daniel Sokolowski <daniel.sokolowski@danols.com>
 * Dave Brooks <dave@bcs.co.nz>
+* David Linke
 * Eduard Bopp <eduard.bopp@aepsil0n.de>
 * Felix Hummel <felix@felixhummel.de>
 * Giel van Schijndel <me@mortis.eu>
diff --git a/docs/nonmult.rst b/docs/nonmult.rst
index 6af34af..3d20c79 100644
--- a/docs/nonmult.rst
+++ b/docs/nonmult.rst
@@ -4,17 +4,17 @@
 Temperature conversion
 ======================
 
-Unlike meters and seconds, fahrenheits, celsius and kelvin are not
-multiplicative units. Temperature is expressed in a system with a
-reference point, and relations between temperature units include
-not only an scaling factor but also an offset. Pint supports these
-type of units and conversions between them. The default definition
-file includes fahrenheits, celsius, kelvin and rankine abbreviated
-as degF, degC, degK, and degR.
+Unlike meters and seconds, the temperature units fahrenheits and
+celsius are non-multiplicative units. These temperature units are 
+expressed in a system with a reference point, and relations between
+temperature units include not only a scaling factor but also an offset.
+Pint supports these type of units and conversions between them. 
+The default definition file includes fahrenheits, celsius,
+kelvin and rankine abbreviated as degF, degC, degK, and degR.
 
 For example, to convert from celsius to fahrenheit:
 
-.. testsetup:: *
+.. testsetup::
 
    from pint import UnitRegistry
    ureg = UnitRegistry()
@@ -24,21 +24,23 @@ For example, to convert from celsius to fahrenheit:
 
    >>> from pint import UnitRegistry
    >>> ureg = UnitRegistry()
-   >>> home = 25.4 * ureg.degC
+   >>> Q_ = ureg.Quantity
+   >>> home = Q_(25.4, ureg.degC)
    >>> print(home.to('degF'))
-   77.72000039999993 degF
+   77.7200004 degF
 
 or to other kelvin or rankine:
 
 .. doctest::
 
-    >>> print(home.to('degK'))
-    298.54999999999995 degK
+    >>> print(home.to('kelvin'))
+    298.55 kelvin
     >>> print(home.to('degR'))
     537.39 degR
 
-Additionally, for every temperature unit in the registry,
-there is also a *delta* counterpart to specify differences.
+Additionally, for every non-multiplicative temperature unit
+in the registry, there is also a *delta* counterpart to specify 
+differences. Absolute units have no *delta* counterpart.
 For example, the change in celsius is equal to the change
 in kelvin, but not in fahrenheit (as the scaling factor
 is different).
@@ -46,20 +48,61 @@ is different).
 .. doctest::
 
    >>> increase = 12.3 * ureg.delta_degC
-   >>> print(increase.to(ureg.delta_degK))
-   12.3 delta_degK
+   >>> print(increase.to(ureg.kelvin))
+   12.3 kelvin
    >>> print(increase.to(ureg.delta_degF))
-   6.833333333333334 delta_degF
+   22.14 delta_degF
 
 ..
-    Subtraction of two temperatures also yields a *delta* unit.
+Subtraction of two temperatures given in offset units yields a *delta* unit:
 
-    .. doctest::
+.. doctest::
+
+    >>> Q_(25.4, ureg.degC) - Q_(10., ureg.degC)
+    <Quantity(15.4, 'delta_degC')>
+
+You can add or subtract a quantity with *delta* unit and a quantity with 
+offset unit:
+
+.. doctest::
+
+    >>> Q_(25.4, ureg.degC) + Q_(10., ureg.delta_degC)
+    <Quantity(35.4, 'degC')>
+    >>> Q_(25.4, ureg.degC) - Q_(10., ureg.delta_degC)
+    <Quantity(15.4, 'degC')>
+
+If you want to add a quantity with absolute unit to one with offset unit, like here
+
+.. doctest::
+
+    >>> heating_rate = 0.5 * ureg.kelvin/ureg.min
+    >>> Q_(10., ureg.degC) + heating_rate * Q_(30, ureg.min)
+    Traceback (most recent call last):
+            ...
+    pint.unit.OffsetUnitCalculusError: Ambiguous operation with offset unit (degC, kelvin).
+
+you have to avoid the ambiguity by either converting the offset unit to the 
+absolute unit before addition
+
+.. doctest::
+
+    >>> Q_(10., ureg.degC).to(ureg.kelvin) + heating_rate * Q_(30, ureg.min)
+    <Quantity(298.15, 'kelvin')>
+
+or convert the absolute unit to a *delta* unit:
+
+.. doctest::
 
-        >>> 25.4 * ureg.degC - 10. * ureg.degC
-        15.4 delta_degC
+    >>> Q_(10., ureg.degC) + heating_rate.to('delta_degC/min') * Q_(30, ureg.min)
+    <Quantity(25.0, 'degC')>
 
-Differences in temperature are multiplicative:
+In contrast to subtraction, the addition of quantities with offset units 
+is ambiguous, e.g. for *10 degC + 100 degC* two different result are reasonable
+depending on the context, *110 degC* or *383.15 °C (= 283.15 K + 373.15 K)*.
+Because of this ambiguity pint raises an error for the addition of two 
+quantities with offset units (since pint-0.6).
+
+Quantities with *delta* units are multiplicative:
 
 .. doctest::
 
@@ -67,7 +110,55 @@ Differences in temperature are multiplicative:
     >>> print(speed.to('delta_degC/second'))
     1.0 delta_degC / second
 
-The parser knows about *delta* units and use them when a temperature unit
+However, multiplication, division and exponentiation of quantities with 
+offset units is problematic just like addition. Pint (since version 0.6)
+will by default raise an error when a quantity with offset unit is used in
+these operations. Due to this quantities with offset units cannot be created
+like other quantities by multiplication of magnitude and unit but have
+to be explicitly created:
+
+.. doctest::
+
+    >>> home = 25.4 * ureg.degC
+    Traceback (most recent call last):
+        ...
+    pint.unit.OffsetUnitCalculusError: Ambiguous operation with offset unit (degC).
+    >>> Q_(25.4, ureg.degC)
+    <Quantity(25.4, 'degC')>
+
+As an alternative to raising an error, pint can be configured to work more
+relaxed via setting the UnitRegistry parameter *autoconvert_offset_to_baseunit*
+to true. In this mode, pint behaves differently:
+
+* Multiplication of a quantity with a single offset unit with order +1 by
+  a number or ndarray yields the quantity in the given unit.
+
+.. doctest::
+
+    >>> ureg = UnitRegistry(autoconvert_offset_to_baseunit = True)
+    >>> T = 25.4 * ureg.degC
+    >>> T
+    <Quantity(25.4, 'degC')>
+    
+* Before all other multiplications, all divisions and in case of 
+  exponentiation [#f1]_ involving quantities with offset-units, pint 
+  will convert the quantities with offset units automatically to the 
+  corresponding base unit before performing the operation. 
+
+    >>> 1/T
+    <Quantity(0.00334952269302, '1 / kelvin')>
+    >>> T * 10 * ureg.meter
+    <Quantity(527.15, 'kelvin * meter')>
+
+You can change the behaviour at any time:
+
+    >>> ureg.autoconvert_offset_to_baseunit = False
+    >>> 1/T
+    Traceback (most recent call last):
+        ...
+    pint.unit.OffsetUnitCalculusError: Ambiguous operation with offset unit (degC).
+
+The parser knows about *delta* units and uses them when a temperature unit
 is found in a multiplicative context. For example, here:
 
 .. doctest::
@@ -92,9 +183,10 @@ You can override this behaviour:
 Note that the magnitude is left unchanged:
 
 .. doctest::
+
     >>> Q_(10, 'degC/meter')
     <Quantity(10, 'delta_degC / meter')>
-    
+
 To define a new temperature, you need to specify the offset. For example,
 this is the definition of the celsius and fahrenheit::
 
@@ -103,3 +195,5 @@ this is the definition of the celsius and fahrenheit::
 
 You do not need to define *delta* units, as they are defined automatically.
 
+.. [#f1] If the exponent is +1, the quantity will not be converted to base
+         unit but remains unchanged.
\ No newline at end of file
diff --git a/docs/serialization.rst b/docs/serialization.rst
index 6e7c133..c1f2a89 100644
--- a/docs/serialization.rst
+++ b/docs/serialization.rst
@@ -19,7 +19,7 @@ The easiest way to do this is by converting the quantity to a string:
    >>> import pint
    >>> ureg = pint.UnitRegistry()
    >>> duration = 24.2 * ureg.years
-   >>> print(duration)
+   >>> duration
    <Quantity(24.2, 'year')>
    >>> serialized = str(duration)
    >>> print(serialized)
@@ -37,7 +37,7 @@ to recover it in another process/machine, you just:
    >>> ureg = pint.UnitRegistry()
    >>> duration = ureg('24.2 year')
    >>> print(duration)
-   <Quantity(24.2, 'year')>
+   24.2 year
 
 Notice that the serialized quantity is likely to be parsed in **another** registry
 as shown in this example. Pint Quantities do not exist on their own but they are
@@ -72,6 +72,7 @@ To unpickle, just
 
     >>> magnitude, units = pickle.loads(serialized)
     >>> ureg.Quantity(magnitude, units)
+    <Quantity(24.2, 'year')>
 
 You can use the same mechanism with any serialization protocol, not only with binary ones.
 (In fact, version 0 of the Pickle protocol is ascii). Other common serialization protocols/packages
diff --git a/docs/tutorial.rst b/docs/tutorial.rst
index 8ad01cd..f89d2cf 100644
--- a/docs/tutorial.rst
+++ b/docs/tutorial.rst
@@ -96,7 +96,6 @@ There are also methods 'to_base_units' and 'ito_base_units' which automatically
    >>> print(height)
    5.75 foot
    >>> height.ito_base_units()
-   <Quantity(1.7526, 'meter')>
    >>> print(height)
    1.7526 meter
 
diff --git a/pint/__init__.py b/pint/__init__.py
index 70adb49..e5452ce 100644
--- a/pint/__init__.py
+++ b/pint/__init__.py
@@ -16,7 +16,8 @@ import os
 import subprocess
 import pkg_resources
 from .formatting import formatter
-from .unit import UnitRegistry, DimensionalityError, UndefinedUnitError, LazyRegistry
+from .unit import (UnitRegistry, DimensionalityError, OffsetUnitCalculusError,
+                   UndefinedUnitError, LazyRegistry)
 from .util import pi_theorem, logger
 
 from .context import Context
diff --git a/pint/quantity.py b/pint/quantity.py
index 305d240..509c0a3 100644
--- a/pint/quantity.py
+++ b/pint/quantity.py
@@ -15,7 +15,8 @@ import operator
 import functools
 
 from .formatting import remove_custom_flags
-from .unit import DimensionalityError, UnitsContainer, UnitDefinition, UndefinedUnitError
+from .unit import (DimensionalityError, OffsetUnitCalculusError,
+                   UnitsContainer, UnitDefinition, UndefinedUnitError)
 from .compat import string_types, ndarray, np, _to_magnitude
 from .util import logger
 
@@ -55,21 +56,8 @@ def _check(q1, other):
     return False
 
 
-def _only_multiplicative_units(q):
-    """Check if the quantity has non-multiplicative units.
-    """
-
-    # Compound units are never multiplicative
-    if len(q.units) != 1:
-        return True
-
-    unit = list(q.units.keys())[0]
-
-    return q._REGISTRY._units[unit].is_multiplicative
-
-
 class _Quantity(object):
-    """Implements a class to describe a physical quantities:
+    """Implements a class to describe a physical quantity:
     the product of a numerical value and a unit of measurement.
 
     :param value: value of the physical quantity to be created.
@@ -294,15 +282,8 @@ class _Quantity(object):
         :param op: operator function (e.g. operator.add, operator.isub)
         :type op: function
         """
-        if _check(self, other):
-            if not self.dimensionality == other.dimensionality:
-                raise DimensionalityError(self.units, other.units,
-                                          self.dimensionality, other.dimensionality)
-            if self._units == other._units:
-                self._magnitude = op(self._magnitude, other._magnitude)
-            else:
-                self._magnitude = op(self._magnitude, other.to(self)._magnitude)
-        else:
+        if not _check(self, other):
+            # other not from same Registry or not a Quantity
             try:
                 other_magnitude = _to_magnitude(other, self.force_ndarray)
             except TypeError:
@@ -318,6 +299,75 @@ class _Quantity(object):
                 self._magnitude = op(self._magnitude, other_magnitude)
             else:
                 raise DimensionalityError(self.units, 'dimensionless')
+            return self
+
+        if not self.dimensionality == other.dimensionality:
+            raise DimensionalityError(self.units, other.units,
+                                      self.dimensionality,
+                                      other.dimensionality)
+
+        # Next we define some variables to make if-clauses more readable.
+        self_non_mul_units = self._get_non_multiplicative_units()
+        is_self_multiplicative = len(self_non_mul_units) == 0
+        if len(self_non_mul_units) == 1:
+            self_non_mul_unit = self_non_mul_units[0]
+        other_non_mul_units = other._get_non_multiplicative_units()
+        is_other_multiplicative = len(other_non_mul_units) == 0
+        if len(other_non_mul_units) == 1:
+            other_non_mul_unit = other_non_mul_units[0]
+
+        # Presence of non-multiplicative units gives rise to several cases.
+        if is_self_multiplicative and is_other_multiplicative:
+            if self._units == other._units:
+                self._magnitude = op(self._magnitude, other._magnitude)
+            # If only self has a delta unit, other determines unit of result.
+            elif self._get_delta_units() and not other._get_delta_units():
+                self._magnitude = op(self._convert_magnitude(other.units),
+                                     other._magnitude)
+                self._units = copy.copy(other.units)
+            else:
+                self._magnitude = op(self._magnitude,
+                                     other.to(self.units)._magnitude)
+
+        elif (op == operator.isub and len(self_non_mul_units) == 1
+                and self.units[self_non_mul_unit] == 1
+                and not other._has_compatible_delta(self_non_mul_unit)):
+            if self.units == other.units:
+                self._magnitude = op(self._magnitude, other._magnitude)
+            else:
+                self._magnitude = op(self._magnitude,
+                                     other.to(self.units)._magnitude)
+            self.units['delta_' + self_non_mul_unit
+                       ] = self.units.pop(self_non_mul_unit)
+
+        elif (op == operator.isub and len(other_non_mul_units) == 1
+                and other.units[other_non_mul_unit] == 1
+                and not self._has_compatible_delta(other_non_mul_unit)):
+            # we convert to self directly since it is multiplicative
+            self._magnitude = op(self._magnitude,
+                                other.to(self.units)._magnitude)
+
+        elif (len(self_non_mul_units) == 1
+                # order of the dimension of offset unit == 1 ?
+                and self._units[self_non_mul_unit] == 1
+                and other._has_compatible_delta(self_non_mul_unit)):
+            tu = copy.copy(self.units)
+            # Replace offset unit in self by the corresponding delta unit.
+            # This is done to prevent a shift by offset in the to()-call.
+            tu['delta_' + self_non_mul_unit] = tu.pop(self_non_mul_unit)
+            self._magnitude = op(self._magnitude, other.to(tu)._magnitude)
+        elif (len(other_non_mul_units) == 1
+                # order of the dimension of offset unit == 1 ?
+                and other._units[other_non_mul_unit] == 1
+                and self._has_compatible_delta(other_non_mul_unit)):
+            tu = copy.copy(other.units)
+            # Replace offset unit in other by the corresponding delta unit.
+            # This is done to prevent a shift by offset in the to()-call.
+            tu['delta_' + other_non_mul_unit] = tu.pop(other_non_mul_unit)
+            self._magnitude = op(self._convert_magnitude(tu), other._magnitude)
+            self._units = copy.copy(other.units)
+        else:
+            raise OffsetUnitCalculusError(self.units, other.units)
 
         return self
 
@@ -329,17 +379,8 @@ class _Quantity(object):
         :param op: operator function (e.g. operator.add, operator.isub)
         :type op: function
         """
-        if _check(self, other):
-            if not self.dimensionality == other.dimensionality:
-                raise DimensionalityError(self.units, other.units,
-                                          self.dimensionality, other.dimensionality)
-            if self._units == other._units:
-                magnitude = op(self._magnitude, other._magnitude)
-            else:
-                magnitude = op(self._magnitude, other.to(self)._magnitude)
-
-            units = copy.copy(self.units)
-        else:
+        if not _check(self, other):
+            # other not from same Registry or not a Quantity
             if _eq(other, 0, True):
                 # If the other value is 0 (but not Quantity 0)
                 # do the operation without checking units.
@@ -354,6 +395,79 @@ class _Quantity(object):
                                _to_magnitude(other, self.force_ndarray))
             else:
                 raise DimensionalityError(self.units, 'dimensionless')
+            return self.__class__(magnitude, units)
+
+        if not self.dimensionality == other.dimensionality:
+            raise DimensionalityError(self.units, other.units,
+                                      self.dimensionality,
+                                      other.dimensionality)
+
+        # Next we define some variables to make if-clauses more readable.
+        self_non_mul_units = self._get_non_multiplicative_units()
+        is_self_multiplicative = len(self_non_mul_units) == 0
+        if len(self_non_mul_units) == 1:
+            self_non_mul_unit = self_non_mul_units[0]
+        other_non_mul_units = other._get_non_multiplicative_units()
+        is_other_multiplicative = len(other_non_mul_units) == 0
+        if len(other_non_mul_units) == 1:
+            other_non_mul_unit = other_non_mul_units[0]
+
+        # Presence of non-multiplicative units gives rise to several cases.
+        if is_self_multiplicative and is_other_multiplicative:
+            if self._units == other._units:
+                magnitude = op(self._magnitude, other._magnitude)
+                units = copy.copy(self.units)
+            # If only self has a delta unit, other determines unit of result.
+            elif self._get_delta_units() and not other._get_delta_units():
+                magnitude = op(self._convert_magnitude(other.units),
+                               other._magnitude)
+                units = copy.copy(other.units)
+            else:
+                units = copy.copy(self.units)
+                magnitude = op(self._magnitude,
+                               other.to(self.units).magnitude)
+
+        elif (op == operator.sub and len(self_non_mul_units) == 1
+                and self.units[self_non_mul_unit] == 1
+                and not other._has_compatible_delta(self_non_mul_unit)):
+            if self.units == other.units:
+                magnitude = op(self._magnitude, other._magnitude)
+            else:
+                magnitude = op(self._magnitude,
+                               other.to(self.units)._magnitude)
+            units = copy.copy(self.units)
+            units['delta_' + self_non_mul_unit] = units.pop(self_non_mul_unit)
+
+        elif (op == operator.sub and len(other_non_mul_units) == 1
+                and other.units[other_non_mul_unit] == 1
+                and not self._has_compatible_delta(other_non_mul_unit)):
+            # we convert to self directly since it is multiplicative
+            magnitude = op(self._magnitude,
+                           other.to(self.units)._magnitude)
+            units = copy.copy(self.units)
+
+        elif (len(self_non_mul_units) == 1
+                # order of the dimension of offset unit == 1 ?
+                and self._units[self_non_mul_unit] == 1
+                and other._has_compatible_delta(self_non_mul_unit)):
+            tu = copy.copy(self.units)
+            # Replace offset unit in self by the corresponding delta unit.
+            # This is done to prevent a shift by offset in the to()-call.
+            tu['delta_' + self_non_mul_unit] = tu.pop(self_non_mul_unit)
+            magnitude = op(self._magnitude, other.to(tu).magnitude)
+            units = copy.copy(self.units)
+        elif (len(other_non_mul_units) == 1
+                # order of the dimension of offset unit == 1 ?
+                and other._units[other_non_mul_unit] == 1
+                and self._has_compatible_delta(other_non_mul_unit)):
+            tu = copy.copy(other.units)
+            # Replace offset unit in other by the corresponding delta unit.
+            # This is done to prevent a shift by offset in the to()-call.
+            tu['delta_' + other_non_mul_unit] = tu.pop(other_non_mul_unit)
+            magnitude = op(self._convert_magnitude(tu), other._magnitude)
+            units = copy.copy(other.units)
+        else:
+            raise OffsetUnitCalculusError(self.units, other.units)
 
         return self.__class__(magnitude, units)
 
@@ -393,24 +507,40 @@ class _Quantity(object):
         if units_op is None:
             units_op = magnitude_op
 
-        if self.__used:
-            if not _only_multiplicative_units(self):
-                self.ito_base_units()
-        else:
-            self.__used = True
-
-        if _check(self, other):
-            if not _only_multiplicative_units(other):
-                other = other.to_base_units()
-            self._magnitude = magnitude_op(self._magnitude, other._magnitude)
-            self._units = units_op(self._units, other._units)
-        else:
+        offset_units_self = self._get_non_multiplicative_units()
+        no_offset_units_self = len(offset_units_self)
+
+        if not _check(self, other):
+            if not self._ok_for_muldiv(no_offset_units_self):
+                raise OffsetUnitCalculusError(self.units,
+                                              getattr(other, 'units', ''))
+            if len(offset_units_self) == 1:
+                if (self.units[offset_units_self[0]] != 1
+                        or magnitude_op not in [operator.mul, operator.imul]):
+                    raise OffsetUnitCalculusError(self.units,
+                                                  getattr(other, 'units', ''))
             try:
                 other_magnitude = _to_magnitude(other, self.force_ndarray)
             except TypeError:
                 return NotImplemented
             self._magnitude = magnitude_op(self._magnitude, other_magnitude)
             self._units = units_op(self._units, UnitsContainer())
+            return self
+
+        if not self._ok_for_muldiv(no_offset_units_self):
+            raise OffsetUnitCalculusError(self.units, other.units)
+        elif no_offset_units_self == 1 and len(self.units) == 1:
+                self.ito_base_units()
+
+        no_offset_units_other = len(other._get_non_multiplicative_units())
+
+        if not other._ok_for_muldiv(no_offset_units_other):
+            raise OffsetUnitCalculusError(self.units, other.units)
+        elif no_offset_units_other == 1 and len(other.units) == 1:
+            other.ito_base_units()
+
+        self._magnitude = magnitude_op(self._magnitude, other._magnitude)
+        self._units = units_op(self._units, other._units)
 
         return self
 
@@ -427,27 +557,46 @@ class _Quantity(object):
         if units_op is None:
             units_op = magnitude_op
 
-        new_self = self
-        if self.__used:
-            if not _only_multiplicative_units(self):
-                new_self = self.to_base_units()
-
-        if _check(self, other):
-            if not _only_multiplicative_units(other):
-                other = other.to_base_units()
-            magnitude = magnitude_op(new_self._magnitude, other._magnitude)
-            units = units_op(new_self._units, other._units)
-        else:
+        offset_units_self = self._get_non_multiplicative_units()
+        no_offset_units_self = len(offset_units_self)
+
+        if not _check(self, other):
+            if not self._ok_for_muldiv(no_offset_units_self):
+                raise OffsetUnitCalculusError(self.units,
+                                              getattr(other, 'units', ''))
+            if len(offset_units_self) == 1:
+                if (self.units[offset_units_self[0]] != 1
+                        or magnitude_op not in [operator.mul, operator.imul]):
+                    raise OffsetUnitCalculusError(self.units,
+                                                  getattr(other, 'units', ''))
             try:
                 other_magnitude = _to_magnitude(other, self.force_ndarray)
             except TypeError:
                 return NotImplemented
-            magnitude = magnitude_op(new_self._magnitude, other_magnitude)
-            units = units_op(new_self._units, UnitsContainer())
 
-        ret = self.__class__(magnitude, units)
-        ret.__used = True
-        return ret
+            magnitude = magnitude_op(self._magnitude, other_magnitude)
+            units = units_op(self._units, UnitsContainer())
+
+            return self.__class__(magnitude, units)
+
+        new_self = self
+
+        if not self._ok_for_muldiv(no_offset_units_self):
+            raise OffsetUnitCalculusError(self.units, other.units)
+        elif no_offset_units_self == 1 and len(self.units) == 1:
+            new_self = self.to_base_units()
+
+        no_offset_units_other = len(other._get_non_multiplicative_units())
+
+        if not other._ok_for_muldiv(no_offset_units_other):
+            raise OffsetUnitCalculusError(self.units, other.units)
+        elif no_offset_units_other == 1 and len(other.units) == 1:
+            other = other.to_base_units()
+
+        magnitude = magnitude_op(new_self._magnitude, other._magnitude)
+        units = units_op(new_self._units, other._units)
+
+        return self.__class__(magnitude, units)
 
     def __imul__(self, other):
         if not isinstance(self._magnitude, ndarray):
@@ -483,6 +632,13 @@ class _Quantity(object):
             other_magnitude = _to_magnitude(other, self.force_ndarray)
         except TypeError:
             return NotImplemented
+
+        no_offset_units_self = len(self._get_non_multiplicative_units())
+        if not self._ok_for_muldiv(no_offset_units_self):
+            raise OffsetUnitCalculusError(self.units, '')
+        elif no_offset_units_self == 1 and len(self.units) == 1:
+            self = self.to_base_units()
+
         return self.__class__(other_magnitude / self._magnitude, 1 / self._units)
 
     def __rfloordiv__(self, other):
@@ -490,6 +646,13 @@ class _Quantity(object):
             other_magnitude = _to_magnitude(other, self.force_ndarray)
         except TypeError:
             return NotImplemented
+
+        no_offset_units_self = len(self._get_non_multiplicative_units())
+        if not self._ok_for_muldiv(no_offset_units_self):
+            raise OffsetUnitCalculusError(self.units, '')
+        elif no_offset_units_self == 1 and len(self.units) == 1:
+            self = self.to_base_units()
+
         return self.__class__(other_magnitude // self._magnitude, 1 / self._units)
 
     __div__ = __truediv__
@@ -505,10 +668,36 @@ class _Quantity(object):
         except TypeError:
             return NotImplemented
         else:
-            if not _only_multiplicative_units(self):
-                self.ito_base_units()
+            if not self._ok_for_muldiv:
+                raise OffsetUnitCalculusError(self.units)
+
+            if isinstance(getattr(other, '_magnitude', other), ndarray):
+                # arrays are refused as exponent, because they would create
+                #  len(array) quanitites of len(set(array)) different units
+                if np.size(other) > 1:
+                    raise DimensionalityError(self.units, 'dimensionless')
+
+            new_self = self
+            if other == 1:
+                return self
+            elif other == 0:
+                self._units = UnitsContainer()
+            else:
+                if not self._is_multiplicative:
+                    if self._REGISTRY.autoconvert_offset_to_baseunit:
+                        self.ito_base_units()
+                    else:
+                        raise OffsetUnitCalculusError(self.units)
+
+                if getattr(other, 'dimensionless', False):
+                    other = other.to_base_units()
+                    self._units **= other.magnitude
+                elif not getattr(other, 'dimensionless', True):
+                    raise DimensionalityError(self.units, 'dimensionless')
+                else:
+                    self._units **= other
+
             self._magnitude **= _to_magnitude(other, self.force_ndarray)
-            self._units **= other
             return self
 
     def __pow__(self, other):
@@ -517,14 +706,51 @@ class _Quantity(object):
         except TypeError:
             return NotImplemented
         else:
+            if not self._ok_for_muldiv:
+                raise OffsetUnitCalculusError(self.units)
+
+            if isinstance(getattr(other, '_magnitude', other), ndarray):
+                # arrays are refused as exponent, because they would create
+                #  len(array) quantities of len(set(array)) different units
+                if np.size(other) > 1:
+                    raise DimensionalityError(self.units, 'dimensionless')
+
             new_self = self
-            if not _only_multiplicative_units(self):
-                new_self = self.to_base_units()
+            if other == 1:
+                return self
+            elif other == 0:
+                units = UnitsContainer()
+            else:
+                if not self._is_multiplicative:
+                    if self._REGISTRY.autoconvert_offset_to_baseunit:
+                        new_self = self.to_base_units()
+                    else:
+                        raise OffsetUnitCalculusError(self.units)
+
+                if getattr(other, 'dimensionless', False):
+                    units = new_self._units ** other.to_base_units().magnitude
+                elif not getattr(other, 'dimensionless', True):
+                    raise DimensionalityError(self.units, 'dimensionless')
+                else:
+                    units = new_self._units ** other
 
             magnitude = new_self._magnitude ** _to_magnitude(other, self.force_ndarray)
-            units = new_self._units ** other
             return self.__class__(magnitude, units)
 
+    def __rpow__(self, other):
+        try:
+            other_magnitude = _to_magnitude(other, self.force_ndarray)
+        except TypeError:
+            return NotImplemented
+        else:
+            if not self.dimensionless:
+                raise DimensionalityError(self.units, 'dimensionless')
+            if isinstance(self._magnitude, ndarray):
+                if np.size(self._magnitude) > 1:
+                    raise DimensionalityError(self.units, 'dimensionless')
+            new_self = self.to_base_units()
+            return other**new_self._magnitude
+
     def __abs__(self):
         return self.__class__(abs(self._magnitude), self._units)
 
@@ -939,3 +1165,58 @@ class _Quantity(object):
                 error = error * abs(self.magnitude)
 
         return self._REGISTRY.Measurement(copy.copy(self.magnitude), error, self.units)
+
+    # methods/properties that help for math operations with offset units
+    @property
+    def _is_multiplicative(self):
+        """Check if the Quantity object has only multiplicative units.
+        """
+        # XXX Turn this into a method/property of _Quantity?
+        return not self._get_non_multiplicative_units()
+
+    def _get_non_multiplicative_units(self):
+        """Return a list of the of non-multiplicative units of the Quantity object
+        """
+        offset_units = [unit for unit in self.units.keys()
+                        if not self._REGISTRY._units[unit].is_multiplicative]
+        return offset_units
+
+    def _get_delta_units(self):
+        """Return list of delta units ot the Quantity object
+        """
+        delta_units = [u for u in self.units.keys() if u.startswith("delta_")]
+        return delta_units
+
+    def _has_compatible_delta(self, unit):
+        """"Check if Quantity object has a delta_unit that is compatible with unit
+        """
+        deltas = self._get_delta_units()
+        if 'delta_' + unit in deltas:
+            return True
+        else:  # Look for delta units with same dimension as the offset unit
+            offset_unit_dim = self._REGISTRY._units[unit].reference
+            for d in deltas:
+                if self._REGISTRY._units[d].reference == offset_unit_dim:
+                    return True
+        return False
+
+    def _ok_for_muldiv(self, no_offset_units=None):
+        """Checks if Quantity object can be multiplied or divided
+
+        :q: quantity object that is checked
+        :no_offset_units: number of offset units in q
+        """
+        is_ok = True
+        if no_offset_units is None:
+            no_offset_units = len(self._get_non_multiplicative_units())
+        if no_offset_units > 1:
+            is_ok = False
+        if no_offset_units == 1:
+            if len(self._units) > 1:
+                is_ok = False
+            if (len(self._units) == 1
+                    and not self._REGISTRY.autoconvert_offset_to_baseunit):
+                is_ok = False
+            if next(iter(self._units.values())) != 1:
+                is_ok = False
+        return is_ok
diff --git a/pint/testsuite/parameterized.py b/pint/testsuite/parameterized.py
new file mode 100644
index 0000000..0b452a6
--- /dev/null
+++ b/pint/testsuite/parameterized.py
@@ -0,0 +1,152 @@
+# -*- coding: utf-8 -*-
+#
+# Adds Parameterized tests for Python's unittest module
+#
+# Code from: parameterizedtestcase, version: 0.1.0
+# Homepage: https://github.com/msabramo/python_unittest_parameterized_test_case
+# Author: Marc Abramowitz, email: marc@marc-abramowitz.com
+# License: MIT
+#
+# Fixed for to work in Python 2 & 3 with "add_metaclass" decorator from six
+# https://pypi.python.org/pypi/six
+# Author: Benjamin Peterson
+# License: MIT 
+#
+# Use like this:
+#
+#    from parameterizedtestcase import ParameterizedTestCase
+#
+#    class MyTests(ParameterizedTestCase):
+#        @ParameterizedTestCase.parameterize(
+#            ("input", "expected_output"),
+#            [
+#                ("2+4", 6),
+#                ("3+5", 8),
+#                ("6*9", 54),
+#            ]
+#        )
+#        def test_eval(self, input, expected_output):
+#            self.assertEqual(eval(input), expected_output)
+
+try:
+    import unittest2 as unittest
+except ImportError:  # pragma: no cover
+    import unittest
+
+from functools import wraps
+import collections
+
+
+def add_metaclass(metaclass):
+    """Class decorator for creating a class with a metaclass."""
+    def wrapper(cls):
+        orig_vars = cls.__dict__.copy()
+        orig_vars.pop('__dict__', None)
+        orig_vars.pop('__weakref__', None)
+        slots = orig_vars.get('__slots__')
+        if slots is not None:
+            if isinstance(slots, str):
+                slots = [slots]
+            for slots_var in slots:
+                orig_vars.pop(slots_var)
+        return metaclass(cls.__name__, cls.__bases__, orig_vars)
+    return wrapper
+
+
+def augment_method_docstring(method, new_class_dict, classname,
+                             param_names, param_values, new_method):
+    param_assignments_str = '; '.join(
+        ['%s = %s' % (k, v) for (k, v) in zip(param_names, param_values)])
+    extra_doc = "%s (%s.%s) [with %s] " % (
+        method.__name__, new_class_dict.get('__module__', '<module>'),
+        classname, param_assignments_str)
+
+    try:
+        new_method.__doc__ = extra_doc + new_method.__doc__
+    except TypeError:  # Catches when new_method.__doc__ is None
+        new_method.__doc__ = extra_doc
+
+
+class ParameterizedTestCaseMetaClass(type):
+    method_counter = {}
+
+    def __new__(meta, classname, bases, class_dict):
+        new_class_dict = {}
+
+        for attr_name, attr_value in list(class_dict.items()):
+            if isinstance(attr_value, collections.Callable) and hasattr(attr_value, 'param_names'):
+                # print("Processing attr_name = %r; attr_value = %r" % (
+                #     attr_name, attr_value))
+
+                method = attr_value
+                param_names = attr_value.param_names
+                data = attr_value.data
+                func_name_format = attr_value.func_name_format
+
+                meta.process_method(
+                    classname, method, param_names, data, new_class_dict,
+                    func_name_format)
+            else:
+                new_class_dict[attr_name] = attr_value
+
+        return type.__new__(meta, classname, bases, new_class_dict)
+
+    @classmethod
+    def process_method(
+            cls, classname, method, param_names, data, new_class_dict,
+            func_name_format):
+        method_counter = cls.method_counter
+
+        for param_values in data:
+            new_method = cls.new_method(method, param_values)
+            method_counter[method.__name__] = \
+                method_counter.get(method.__name__, 0) + 1
+            case_data = dict(list(zip(param_names, param_values)))
+            case_data['func_name'] = method.__name__
+            case_data['case_num'] = method_counter[method.__name__]
+
+            new_method.__name__ = func_name_format.format(**case_data)
+
+            augment_method_docstring(
+                method, new_class_dict, classname,
+                param_names, param_values, new_method)
+            new_class_dict[new_method.__name__] = new_method
+
+    @classmethod
+    def new_method(cls, method, param_values):
+        @wraps(method)
+        def new_method(self):
+            return method(self, *param_values)
+
+        return new_method
+
+@add_metaclass(ParameterizedTestCaseMetaClass)
+class ParameterizedTestMixin(object):
+    @classmethod
+    def parameterize(cls, param_names, data,
+                     func_name_format='{func_name}_{case_num:05d}'):
+        """Decorator for parameterizing a test method - example:
+
+        @ParameterizedTestCase.parameterize(
+            ("isbn", "expected_title"), [
+                ("0262033844", "Introduction to Algorithms"),
+                ("0321558146", "Campbell Essential Biology")])
+
+        """
+
+        def decorator(func):
+            @wraps(func)
+            def newfunc(*arg, **kwargs):
+                return func(*arg, **kwargs)
+
+            newfunc.param_names = param_names
+            newfunc.data = data
+            newfunc.func_name_format = func_name_format
+
+            return newfunc
+
+        return decorator
+
+@add_metaclass(ParameterizedTestCaseMetaClass)
+class ParameterizedTestCase(unittest.TestCase, ParameterizedTestMixin):
+    pass
diff --git a/pint/testsuite/test_issues.py b/pint/testsuite/test_issues.py
index 958fc0e..1e23661 100644
--- a/pint/testsuite/test_issues.py
+++ b/pint/testsuite/test_issues.py
@@ -16,6 +16,9 @@ class TestIssues(QuantityTestCase):
 
     FORCE_NDARRAY = False
 
+    def setup(self):
+        self.ureg.autoconvert_offset_to_baseunit = False
+
     @unittest.expectedFailure
     def test_issue25(self):
         x = ParserHelper.from_string('10 %')
@@ -92,9 +95,9 @@ class TestIssues(QuantityTestCase):
     def test_issue66b(self):
         ureg = UnitRegistry()
         self.assertEqual(ureg.get_base_units(ureg.kelvin.units),
-                         (None, UnitsContainer({'kelvin': 1})))
+                         (1.0, UnitsContainer({'kelvin': 1})))
         self.assertEqual(ureg.get_base_units(ureg.degC.units),
-                         (None, UnitsContainer({'kelvin': 1})))
+                         (1.0, UnitsContainer({'kelvin': 1})))
 
     def test_issue69(self):
         ureg = UnitRegistry()
@@ -127,11 +130,12 @@ class TestIssues(QuantityTestCase):
         self.assertQuantityAlmostEqual(va.to_base_units(), vb.to_base_units())
 
     def test_issue86(self):
+        ureg = self.ureg
+        ureg.autoconvert_offset_to_baseunit = True
+
         def parts(q):
             return q.magnitude, q.units
 
-        ureg = UnitRegistry()
-
         q1 = 10. * ureg.degC
         q2 = 10. * ureg.kelvin
 
@@ -158,7 +162,6 @@ class TestIssues(QuantityTestCase):
         self.assertEqual(parts(q1 / q3), (k1m / q3m, k1u / q3u))
         self.assertEqual(parts(q3 * q1), (q3m * k1m, q3u * k1u))
         self.assertEqual(parts(q3 / q1), (q3m / k1m, q3u / k1u))
-        self.assertEqual(parts(q1 **  1), (k1m **  1, k1u **  1))
         self.assertEqual(parts(q1 ** -1), (k1m ** -1, k1u ** -1))
         self.assertEqual(parts(q1 **  2), (k1m **  2, k1u **  2))
         self.assertEqual(parts(q1 ** -2), (k1m ** -2, k1u ** -2))
@@ -177,8 +180,10 @@ class TestIssues(QuantityTestCase):
         self.assertQuantityAlmostEqual(v1.to_base_units(), v2)
         self.assertQuantityAlmostEqual(v1.to_base_units(), v2.to_base_units())
 
+    @unittest.expectedFailure
     def test_issue86c(self):
         ureg = self.ureg
+        ureg.autoconvert_offset_to_baseunit = True
         T = ureg.degC
         T = 100. * T
         self.assertQuantityAlmostEqual(ureg.k*2*T, ureg.k*(2*T))
@@ -263,21 +268,21 @@ class TestIssuesNP(QuantityTestCase):
         q = ureg.meter * x
         self.assertIsInstance(q, ureg.Quantity)
         np.testing.assert_array_equal(q.magnitude, x)
-        self.assertEquals(q.units, ureg.meter.units)
+        self.assertEqual(q.units, ureg.meter.units)
         q = x * ureg.meter
         self.assertIsInstance(q, ureg.Quantity)
         np.testing.assert_array_equal(q.magnitude, x)
-        self.assertEquals(q.units, ureg.meter.units)
+        self.assertEqual(q.units, ureg.meter.units)
 
         m = np.ma.masked_array(2 * np.ones(3,3))
         qq = q * m
         self.assertIsInstance(qq, ureg.Quantity)
         np.testing.assert_array_equal(qq.magnitude, x * m)
-        self.assertEquals(qq.units, ureg.meter.units)
+        self.assertEqual(qq.units, ureg.meter.units)
         qq = m * q
         self.assertIsInstance(qq, ureg.Quantity)
         np.testing.assert_array_equal(qq.magnitude, x * m)
-        self.assertEquals(qq.units, ureg.meter.units)
+        self.assertEqual(qq.units, ureg.meter.units)
 
     @unittest.expectedFailure
     def test_issue39(self):
@@ -286,21 +291,21 @@ class TestIssuesNP(QuantityTestCase):
         q = ureg.meter * x
         self.assertIsInstance(q, ureg.Quantity)
         np.testing.assert_array_equal(q.magnitude, x)
-        self.assertEquals(q.units, ureg.meter.units)
+        self.assertEqual(q.units, ureg.meter.units)
         q = x * ureg.meter
         self.assertIsInstance(q, ureg.Quantity)
         np.testing.assert_array_equal(q.magnitude, x)
-        self.assertEquals(q.units, ureg.meter.units)
+        self.assertEqual(q.units, ureg.meter.units)
 
         m = np.matrix(2 * np.ones(3,3))
         qq = q * m
         self.assertIsInstance(qq, ureg.Quantity)
         np.testing.assert_array_equal(qq.magnitude, x * m)
-        self.assertEquals(qq.units, ureg.meter.units)
+        self.assertEqual(qq.units, ureg.meter.units)
         qq = m * q
         self.assertIsInstance(qq, ureg.Quantity)
         np.testing.assert_array_equal(qq.magnitude, x * m)
-        self.assertEquals(qq.units, ureg.meter.units)
+        self.assertEqual(qq.units, ureg.meter.units)
 
     def test_issue44(self):
         ureg = UnitRegistry()
diff --git a/pint/testsuite/test_numpy.py b/pint/testsuite/test_numpy.py
index 3c4a5e0..5ce0f04 100644
--- a/pint/testsuite/test_numpy.py
+++ b/pint/testsuite/test_numpy.py
@@ -2,6 +2,10 @@
 
 from __future__ import division, unicode_literals, print_function, absolute_import
 
+import copy
+import operator as op
+
+from pint import DimensionalityError
 from pint.compat import np, unittest
 from pint.testsuite import QuantityTestCase, helpers
 from pint.testsuite.test_umath import TestUFuncs
@@ -399,3 +403,39 @@ class TestBitTwiddlingUfuncs(TestUFuncs):
                     (self.qless, 2),
                     (self.q1, self.q2, self.qs, ),
                     'same')
+
+
+class TestNDArrayQunatityMath(QuantityTestCase):
+
+    @helpers.requires_numpy()
+    def test_exponentiation_array_exp(self):
+        arr = np.array(range(3), dtype=np.float)
+        q = self.Q_(arr, None)
+
+        for op_ in [op.pow, op.ipow]:
+            q_cp = copy.copy(q)
+            self.assertRaises(DimensionalityError, op_, 2., q_cp)
+            arr_cp = copy.copy(arr)
+            arr_cp = copy.copy(arr)
+            q_cp = copy.copy(q)
+            self.assertRaises(DimensionalityError, op_, q_cp, arr_cp)
+            q_cp = copy.copy(q)
+            q2_cp = copy.copy(q)
+            self.assertRaises(DimensionalityError, op_, q_cp, q2_cp)
+
+    @unittest.expectedFailure
+    @helpers.requires_numpy()
+    def test_exponentiation_array_exp_2(self):
+        arr = np.array(range(3), dtype=np.float)
+        #q = self.Q_(copy.copy(arr), None)
+        q = self.Q_(copy.copy(arr), 'meter')
+        arr_cp = copy.copy(arr)
+        q_cp = copy.copy(q)
+        # this fails as expected since numpy 1.8.0 but...
+        self.assertRaises(DimensionalityError, op.pow, arr_cp, q_cp)
+        # ..not for op.ipow !
+        # q_cp is treated as if it is an array. The units are ignored.
+        # _Quantity.__ipow__ is never called
+        arr_cp = copy.copy(arr)
+        q_cp = copy.copy(q)
+        self.assertRaises(DimensionalityError, op.ipow, arr_cp, q_cp)
diff --git a/pint/testsuite/test_quantity.py b/pint/testsuite/test_quantity.py
index ed63b87..72122e1 100644
--- a/pint/testsuite/test_quantity.py
+++ b/pint/testsuite/test_quantity.py
@@ -6,10 +6,11 @@ import copy
 import math
 import operator as op
 
-from pint import DimensionalityError, UnitRegistry
+from pint import DimensionalityError, OffsetUnitCalculusError, UnitRegistry
 from pint.unit import UnitsContainer
-from pint.compat import string_types, PYTHON3, np
+from pint.compat import string_types, PYTHON3, np, unittest
 from pint.testsuite import QuantityTestCase, helpers
+from pint.testsuite.parameterized import ParameterizedTestCase
 
 
 class TestQuantity(QuantityTestCase):
@@ -149,7 +150,7 @@ class TestQuantity(QuantityTestCase):
         # Conversions with single units take a different codepath than
         # Conversions with more than one unit.
         src_dst1 = UnitsContainer(meter=1), UnitsContainer(inch=1)
-        src_dst2 = UnitsContainer(meter=1, seconds=-1), UnitsContainer(inch=1, minutes=-1)
+        src_dst2 = UnitsContainer(meter=1, second=-1), UnitsContainer(inch=1, minute=-1)
         for src, dst in (src_dst1, src_dst2):
             a = np.ones((3, 1))
             ac = np.ones((3, 1))
@@ -209,14 +210,12 @@ class TestQuantity(QuantityTestCase):
 
 
     def test_offset_delta(self):
-        self.assertQuantityAlmostEqual(self.Q_(0, 'delta_kelvin').to('delta_kelvin'), self.Q_(0, 'delta_kelvin'))
-        self.assertQuantityAlmostEqual(self.Q_(0, 'delta_degC').to('delta_kelvin'), self.Q_(0, 'delta_kelvin'))
-        self.assertQuantityAlmostEqual(self.Q_(0, 'delta_degF').to('delta_kelvin'), self.Q_(0, 'delta_kelvin'), rtol=0.01)
-
-        self.assertQuantityAlmostEqual(self.Q_(100, 'delta_kelvin').to('delta_kelvin'), self.Q_(100, 'delta_kelvin'))
-        self.assertQuantityAlmostEqual(self.Q_(100, 'delta_kelvin').to('delta_degC'), self.Q_(100, 'delta_degC'))
-        self.assertQuantityAlmostEqual(self.Q_(100, 'delta_kelvin').to('delta_degF'), self.Q_(180, 'delta_degF'), rtol=0.01)
-        self.assertQuantityAlmostEqual(self.Q_(100, 'delta_degF').to('delta_kelvin'), self.Q_(55.55555556, 'delta_kelvin'), rtol=0.01)
+        self.assertQuantityAlmostEqual(self.Q_(0, 'delta_degC').to('kelvin'), self.Q_(0, 'kelvin'))
+        self.assertQuantityAlmostEqual(self.Q_(0, 'delta_degF').to('kelvin'), self.Q_(0, 'kelvin'), rtol=0.01)
+
+        self.assertQuantityAlmostEqual(self.Q_(100, 'kelvin').to('delta_degC'), self.Q_(100, 'delta_degC'))
+        self.assertQuantityAlmostEqual(self.Q_(100, 'kelvin').to('delta_degF'), self.Q_(180, 'delta_degF'), rtol=0.01)
+        self.assertQuantityAlmostEqual(self.Q_(100, 'delta_degF').to('kelvin'), self.Q_(55.55555556, 'kelvin'), rtol=0.01)
         self.assertQuantityAlmostEqual(self.Q_(100, 'delta_degC').to('delta_degF'), self.Q_(180, 'delta_degF'), rtol=0.01)
         self.assertQuantityAlmostEqual(self.Q_(100, 'delta_degF').to('delta_degC'), self.Q_(55.55555556, 'delta_degC'), rtol=0.01)
 
@@ -450,3 +449,547 @@ class TestDimensionsWithDefaultRegistry(TestDimensions):
         from pint import _DEFAULT_REGISTRY
         cls.ureg = _DEFAULT_REGISTRY
         cls.Q_ = cls.ureg.Quantity
+
+
+class TestOffsetUnitMath(QuantityTestCase, ParameterizedTestCase):
+
+    def setup(self):
+        self.ureg.autoconvert_offset_to_baseunit = False
+        self.ureg.default_as_delta = True
+
+    additions = [
+        # --- input tuple -------------------- | -- expected result --
+        (((100, 'kelvin'), (10, 'kelvin')),      (110, 'kelvin')),
+        (((100, 'kelvin'), (10, 'degC')),        'error'),
+        (((100, 'kelvin'), (10, 'degF')),        'error'),
+        (((100, 'kelvin'), (10, 'degR')),        (105.56, 'kelvin')),
+        (((100, 'kelvin'), (10, 'delta_degC')),  (110, 'kelvin')),
+        (((100, 'kelvin'), (10, 'delta_degF')),  (105.56, 'kelvin')),
+
+        (((100, 'degC'), (10, 'kelvin')),      'error'),
+        (((100, 'degC'), (10, 'degC')),        'error'),
+        (((100, 'degC'), (10, 'degF')),        'error'),
+        (((100, 'degC'), (10, 'degR')),        'error'),
+        (((100, 'degC'), (10, 'delta_degC')),  (110, 'degC')),
+        (((100, 'degC'), (10, 'delta_degF')),  (105.56, 'degC')),
+
+        (((100, 'degF'), (10, 'kelvin')),      'error'),
+        (((100, 'degF'), (10, 'degC')),        'error'),
+        (((100, 'degF'), (10, 'degF')),        'error'),
+        (((100, 'degF'), (10, 'degR')),        'error'),
+        (((100, 'degF'), (10, 'delta_degC')),  (118, 'degF')),
+        (((100, 'degF'), (10, 'delta_degF')),  (110, 'degF')),
+
+        (((100, 'degR'), (10, 'kelvin')),      (118, 'degR')),
+        (((100, 'degR'), (10, 'degC')),        'error'),
+        (((100, 'degR'), (10, 'degF')),        'error'),
+        (((100, 'degR'), (10, 'degR')),        (110, 'degR')),
+        (((100, 'degR'), (10, 'delta_degC')),  (118, 'degR')),
+        (((100, 'degR'), (10, 'delta_degF')),  (110, 'degR')),
+
+        (((100, 'delta_degC'), (10, 'kelvin')),     (110, 'kelvin')),
+        (((100, 'delta_degC'), (10, 'degC')),       (110, 'degC')),
+        (((100, 'delta_degC'), (10, 'degF')),       (190, 'degF')),
+        (((100, 'delta_degC'), (10, 'degR')),       (190, 'degR')),
+        (((100, 'delta_degC'), (10, 'delta_degC')), (110, 'delta_degC')),
+        (((100, 'delta_degC'), (10, 'delta_degF')), (105.56, 'delta_degC')),
+
+        (((100, 'delta_degF'), (10, 'kelvin')),     (65.56, 'kelvin')),
+        (((100, 'delta_degF'), (10, 'degC')),       (65.56, 'degC')),
+        (((100, 'delta_degF'), (10, 'degF')),       (110, 'degF')),
+        (((100, 'delta_degF'), (10, 'degR')),       (110, 'degR')),
+        (((100, 'delta_degF'), (10, 'delta_degC')), (118, 'delta_degF')),
+        (((100, 'delta_degF'), (10, 'delta_degF')), (110, 'delta_degF')),
+        ]
+
+    @ParameterizedTestCase.parameterize(("input", "expected_output"),
+                                        additions)
+    def test_addition(self, input_tuple, expected):
+        self.ureg.autoconvert_offset_to_baseunit = False
+        qin1, qin2 = input_tuple
+        q1, q2 = self.Q_(*qin1), self.Q_(*qin2)
+        # update input tuple with new values to have correct values on failure
+        input_tuple = q1, q2
+        if expected == 'error':
+            self.assertRaises(OffsetUnitCalculusError, op.add, q1, q2)
+        else:
+            expected = self.Q_(*expected)
+            self.assertEqual(op.add(q1, q2).units, expected.units)
+            self.assertQuantityAlmostEqual(op.add(q1, q2), expected,
+                                           atol=0.01)
+
+    @helpers.requires_numpy()
+    @ParameterizedTestCase.parameterize(("input", "expected_output"),
+                                        additions)
+    def test_inplace_addition(self, input_tuple, expected):
+        self.ureg.autoconvert_offset_to_baseunit = False
+        (q1v, q1u), (q2v, q2u) = input_tuple
+        # update input tuple with new values to have correct values on failure
+        input_tuple = ((np.array([q1v]*2, dtype=np.float), q1u),
+                       (np.array([q2v]*2, dtype=np.float), q2u))
+        Q_ = self.Q_
+        qin1, qin2 = input_tuple
+        q1, q2 = Q_(*qin1), Q_(*qin2)
+        q1_cp = copy.copy(q1)
+        if expected == 'error':
+            self.assertRaises(OffsetUnitCalculusError, op.iadd, q1_cp, q2)
+        else:
+            expected = np.array([expected[0]]*2, dtype=np.float), expected[1]
+            self.assertEqual(op.iadd(q1_cp, q2).units, Q_(*expected).units)
+            q1_cp = copy.copy(q1)
+            self.assertQuantityAlmostEqual(op.iadd(q1_cp, q2), Q_(*expected),
+                                           atol=0.01)
+
+    subtractions = [
+        (((100, 'kelvin'), (10, 'kelvin')),      (90, 'kelvin')),
+        (((100, 'kelvin'), (10, 'degC')),        (-183.15, 'kelvin')),
+        (((100, 'kelvin'), (10, 'degF')),        (-160.93, 'kelvin')),
+        (((100, 'kelvin'), (10, 'degR')),        (94.44, 'kelvin')),
+        (((100, 'kelvin'), (10, 'delta_degC')),  (90, 'kelvin')),
+        (((100, 'kelvin'), (10, 'delta_degF')),  (94.44, 'kelvin')),
+
+        (((100, 'degC'), (10, 'kelvin')),      (363.15, 'delta_degC')),
+        (((100, 'degC'), (10, 'degC')),        (90, 'delta_degC')),
+        (((100, 'degC'), (10, 'degF')),        (112.22, 'delta_degC')),
+        (((100, 'degC'), (10, 'degR')),        (367.59, 'delta_degC')),
+        (((100, 'degC'), (10, 'delta_degC')),  (90, 'degC')),
+        (((100, 'degC'), (10, 'delta_degF')),  (94.44, 'degC')),
+
+        (((100, 'degF'), (10, 'kelvin')),      (541.67, 'delta_degF')),
+        (((100, 'degF'), (10, 'degC')),        (50, 'delta_degF')),
+        (((100, 'degF'), (10, 'degF')),        (90, 'delta_degF')),
+        (((100, 'degF'), (10, 'degR')),        (549.67, 'delta_degF')),
+        (((100, 'degF'), (10, 'delta_degC')),  (82, 'degF')),
+        (((100, 'degF'), (10, 'delta_degF')),  (90, 'degF')),
+
+        (((100, 'degR'), (10, 'kelvin')),      (82, 'degR')),
+        (((100, 'degR'), (10, 'degC')),        (-409.67, 'degR')),
+        (((100, 'degR'), (10, 'degF')),        (-369.67, 'degR')),
+        (((100, 'degR'), (10, 'degR')),        (90, 'degR')),
+        (((100, 'degR'), (10, 'delta_degC')),  (82, 'degR')),
+        (((100, 'degR'), (10, 'delta_degF')),  (90, 'degR')),
+
+        (((100, 'delta_degC'), (10, 'kelvin')),     (90, 'kelvin')),
+        (((100, 'delta_degC'), (10, 'degC')),       (90, 'degC')),
+        (((100, 'delta_degC'), (10, 'degF')),       (170, 'degF')),
+        (((100, 'delta_degC'), (10, 'degR')),       (170, 'degR')),
+        (((100, 'delta_degC'), (10, 'delta_degC')), (90, 'delta_degC')),
+        (((100, 'delta_degC'), (10, 'delta_degF')), (94.44, 'delta_degC')),
+
+        (((100, 'delta_degF'), (10, 'kelvin')),     (45.56, 'kelvin')),
+        (((100, 'delta_degF'), (10, 'degC')),       (45.56, 'degC')),
+        (((100, 'delta_degF'), (10, 'degF')),       (90, 'degF')),
+        (((100, 'delta_degF'), (10, 'degR')),       (90, 'degR')),
+        (((100, 'delta_degF'), (10, 'delta_degC')), (82, 'delta_degF')),
+        (((100, 'delta_degF'), (10, 'delta_degF')), (90, 'delta_degF')),
+        ]
+
+    @ParameterizedTestCase.parameterize(("input", "expected_output"),
+                                        subtractions)
+    def test_subtraction(self, input_tuple, expected):
+        self.ureg.autoconvert_offset_to_baseunit = False
+        qin1, qin2 = input_tuple
+        q1, q2 = self.Q_(*qin1), self.Q_(*qin2)
+        input_tuple = q1, q2
+        if expected == 'error':
+            self.assertRaises(OffsetUnitCalculusError, op.sub, q1, q2)
+        else:
+            expected = self.Q_(*expected)
+            self.assertEqual(op.sub(q1, q2).units, expected.units)
+            self.assertQuantityAlmostEqual(op.sub(q1, q2), expected,
+                                           atol=0.01)
+
+#    @unittest.expectedFailure
+    @helpers.requires_numpy()
+    @ParameterizedTestCase.parameterize(("input", "expected_output"),
+                                        subtractions)
+    def test_inplace_subtraction(self, input_tuple, expected):
+        self.ureg.autoconvert_offset_to_baseunit = False
+        (q1v, q1u), (q2v, q2u) = input_tuple
+        # update input tuple with new values to have correct values on failure
+        input_tuple = ((np.array([q1v]*2, dtype=np.float), q1u),
+                       (np.array([q2v]*2, dtype=np.float), q2u))
+        Q_ = self.Q_
+        qin1, qin2 = input_tuple
+        q1, q2 = Q_(*qin1), Q_(*qin2)
+        q1_cp = copy.copy(q1)
+        if expected == 'error':
+            self.assertRaises(OffsetUnitCalculusError, op.isub, q1_cp, q2)
+        else:
+            expected = np.array([expected[0]]*2, dtype=np.float), expected[1]
+            self.assertEqual(op.isub(q1_cp, q2).units, Q_(*expected).units)
+            q1_cp = copy.copy(q1)
+            self.assertQuantityAlmostEqual(op.isub(q1_cp, q2), Q_(*expected),
+                                           atol=0.01)
+
+    multiplications = [
+        (((100, 'kelvin'), (10, 'kelvin')),     (1000, 'kelvin**2')),
+        (((100, 'kelvin'), (10, 'degC')),       'error'),
+        (((100, 'kelvin'), (10, 'degF')),       'error'),
+        (((100, 'kelvin'), (10, 'degR')),       (1000, 'kelvin*degR')),
+        (((100, 'kelvin'), (10, 'delta_degC')), (1000, 'kelvin*delta_degC')),
+        (((100, 'kelvin'), (10, 'delta_degF')), (1000, 'kelvin*delta_degF')),
+
+        (((100, 'degC'), (10, 'kelvin')),      'error'),
+        (((100, 'degC'), (10, 'degC')),        'error'),
+        (((100, 'degC'), (10, 'degF')),        'error'),
+        (((100, 'degC'), (10, 'degR')),        'error'),
+        (((100, 'degC'), (10, 'delta_degC')),  'error'),
+        (((100, 'degC'), (10, 'delta_degF')),  'error'),
+
+        (((100, 'degF'), (10, 'kelvin')),      'error'),
+        (((100, 'degF'), (10, 'degC')),        'error'),
+        (((100, 'degF'), (10, 'degF')),        'error'),
+        (((100, 'degF'), (10, 'degR')),        'error'),
+        (((100, 'degF'), (10, 'delta_degC')),  'error'),
+        (((100, 'degF'), (10, 'delta_degF')),  'error'),
+
+        (((100, 'degR'), (10, 'kelvin')),      (1000, 'degR*kelvin')),
+        (((100, 'degR'), (10, 'degC')),        'error'),
+        (((100, 'degR'), (10, 'degF')),        'error'),
+        (((100, 'degR'), (10, 'degR')),        (1000, 'degR**2')),
+        (((100, 'degR'), (10, 'delta_degC')),  (1000, 'degR*delta_degC')),
+        (((100, 'degR'), (10, 'delta_degF')),  (1000, 'degR*delta_degF')),
+
+        (((100, 'delta_degC'), (10, 'kelvin')),     (1000, 'delta_degC*kelvin')),
+        (((100, 'delta_degC'), (10, 'degC')),       'error'),
+        (((100, 'delta_degC'), (10, 'degF')),       'error'),
+        (((100, 'delta_degC'), (10, 'degR')),       (1000, 'delta_degC*degR')),
+        (((100, 'delta_degC'), (10, 'delta_degC')), (1000, 'delta_degC**2')),
+        (((100, 'delta_degC'), (10, 'delta_degF')), (1000, 'delta_degC*delta_degF')),
+
+        (((100, 'delta_degF'), (10, 'kelvin')),     (1000, 'delta_degF*kelvin')),
+        (((100, 'delta_degF'), (10, 'degC')),       'error'),
+        (((100, 'delta_degF'), (10, 'degF')),       'error'),
+        (((100, 'delta_degF'), (10, 'degR')),       (1000, 'delta_degF*degR')),
+        (((100, 'delta_degF'), (10, 'delta_degC')), (1000, 'delta_degF*delta_degC')),
+        (((100, 'delta_degF'), (10, 'delta_degF')), (1000, 'delta_degF**2')),
+        ]
+
+    @ParameterizedTestCase.parameterize(("input", "expected_output"),
+                                        multiplications)
+    def test_multiplication(self, input_tuple, expected):
+        self.ureg.autoconvert_offset_to_baseunit = False
+        qin1, qin2 = input_tuple
+        q1, q2 = self.Q_(*qin1), self.Q_(*qin2)
+        input_tuple = q1, q2
+        if expected == 'error':
+            self.assertRaises(OffsetUnitCalculusError, op.mul, q1, q2)
+        else:
+            expected = self.Q_(*expected)
+            self.assertEqual(op.mul(q1, q2).units, expected.units)
+            self.assertQuantityAlmostEqual(op.mul(q1, q2), expected,
+                                           atol=0.01)
+
+    @helpers.requires_numpy()
+    @ParameterizedTestCase.parameterize(("input", "expected_output"),
+                                        multiplications)
+    def test_inplace_multiplication(self, input_tuple, expected):
+        self.ureg.autoconvert_offset_to_baseunit = False
+        (q1v, q1u), (q2v, q2u) = input_tuple
+        # update input tuple with new values to have correct values on failure
+        input_tuple = ((np.array([q1v]*2, dtype=np.float), q1u),
+                       (np.array([q2v]*2, dtype=np.float), q2u))
+        Q_ = self.Q_
+        qin1, qin2 = input_tuple
+        q1, q2 = Q_(*qin1), Q_(*qin2)
+        q1_cp = copy.copy(q1)
+        if expected == 'error':
+            self.assertRaises(OffsetUnitCalculusError, op.imul, q1_cp, q2)
+        else:
+            expected = np.array([expected[0]]*2, dtype=np.float), expected[1]
+            self.assertEqual(op.imul(q1_cp, q2).units, Q_(*expected).units)
+            q1_cp = copy.copy(q1)
+            self.assertQuantityAlmostEqual(op.imul(q1_cp, q2), Q_(*expected),
+                                           atol=0.01)
+
+    divisions = [
+        (((100, 'kelvin'), (10, 'kelvin')),     (10, '')),
+        (((100, 'kelvin'), (10, 'degC')),       'error'),
+        (((100, 'kelvin'), (10, 'degF')),       'error'),
+        (((100, 'kelvin'), (10, 'degR')),       (10, 'kelvin/degR')),
+        (((100, 'kelvin'), (10, 'delta_degC')), (10, 'kelvin/delta_degC')),
+        (((100, 'kelvin'), (10, 'delta_degF')), (10, 'kelvin/delta_degF')),
+
+        (((100, 'degC'), (10, 'kelvin')),      'error'),
+        (((100, 'degC'), (10, 'degC')),        'error'),
+        (((100, 'degC'), (10, 'degF')),        'error'),
+        (((100, 'degC'), (10, 'degR')),        'error'),
+        (((100, 'degC'), (10, 'delta_degC')),  'error'),
+        (((100, 'degC'), (10, 'delta_degF')),  'error'),
+
+        (((100, 'degF'), (10, 'kelvin')),      'error'),
+        (((100, 'degF'), (10, 'degC')),        'error'),
+        (((100, 'degF'), (10, 'degF')),        'error'),
+        (((100, 'degF'), (10, 'degR')),        'error'),
+        (((100, 'degF'), (10, 'delta_degC')),  'error'),
+        (((100, 'degF'), (10, 'delta_degF')),  'error'),
+
+        (((100, 'degR'), (10, 'kelvin')),      (10, 'degR/kelvin')),
+        (((100, 'degR'), (10, 'degC')),        'error'),
+        (((100, 'degR'), (10, 'degF')),        'error'),
+        (((100, 'degR'), (10, 'degR')),        (10, '')),
+        (((100, 'degR'), (10, 'delta_degC')),  (10, 'degR/delta_degC')),
+        (((100, 'degR'), (10, 'delta_degF')),  (10, 'degR/delta_degF')),
+
+        (((100, 'delta_degC'), (10, 'kelvin')),     (10, 'delta_degC/kelvin')),
+        (((100, 'delta_degC'), (10, 'degC')),       'error'),
+        (((100, 'delta_degC'), (10, 'degF')),       'error'),
+        (((100, 'delta_degC'), (10, 'degR')),       (10, 'delta_degC/degR')),
+        (((100, 'delta_degC'), (10, 'delta_degC')), (10, '')),
+        (((100, 'delta_degC'), (10, 'delta_degF')), (10, 'delta_degC/delta_degF')),
+
+        (((100, 'delta_degF'), (10, 'kelvin')),     (10, 'delta_degF/kelvin')),
+        (((100, 'delta_degF'), (10, 'degC')),       'error'),
+        (((100, 'delta_degF'), (10, 'degF')),       'error'),
+        (((100, 'delta_degF'), (10, 'degR')),       (10, 'delta_degF/degR')),
+        (((100, 'delta_degF'), (10, 'delta_degC')), (10, 'delta_degF/delta_degC')),
+        (((100, 'delta_degF'), (10, 'delta_degF')), (10, '')),
+        ]
+
+    @ParameterizedTestCase.parameterize(("input", "expected_output"),
+                                        divisions)
+    def test_truedivision(self, input_tuple, expected):
+        self.ureg.autoconvert_offset_to_baseunit = False
+        qin1, qin2 = input_tuple
+        q1, q2 = self.Q_(*qin1), self.Q_(*qin2)
+        input_tuple = q1, q2
+        if expected == 'error':
+            self.assertRaises(OffsetUnitCalculusError, op.truediv, q1, q2)
+        else:
+            expected = self.Q_(*expected)
+            self.assertEqual(op.truediv(q1, q2).units, expected.units)
+            self.assertQuantityAlmostEqual(op.truediv(q1, q2), expected,
+                                           atol=0.01)
+
+    @helpers.requires_numpy()
+    @ParameterizedTestCase.parameterize(("input", "expected_output"),
+                                        divisions)
+    def test_inplace_truedivision(self, input_tuple, expected):
+        self.ureg.autoconvert_offset_to_baseunit = False
+        (q1v, q1u), (q2v, q2u) = input_tuple
+        # update input tuple with new values to have correct values on failure
+        input_tuple = ((np.array([q1v]*2, dtype=np.float), q1u),
+                       (np.array([q2v]*2, dtype=np.float), q2u))
+        Q_ = self.Q_
+        qin1, qin2 = input_tuple
+        q1, q2 = Q_(*qin1), Q_(*qin2)
+        q1_cp = copy.copy(q1)
+        if expected == 'error':
+            self.assertRaises(OffsetUnitCalculusError, op.itruediv, q1_cp, q2)
+        else:
+            expected = np.array([expected[0]]*2, dtype=np.float), expected[1]
+            self.assertEqual(op.itruediv(q1_cp, q2).units, Q_(*expected).units)
+            q1_cp = copy.copy(q1)
+            self.assertQuantityAlmostEqual(op.itruediv(q1_cp, q2),
+                                           Q_(*expected), atol=0.01)
+
+    multiplications_with_autoconvert_to_baseunit = [
+        (((100, 'kelvin'), (10, 'degC')),     (28315., 'kelvin**2')),
+        (((100, 'kelvin'), (10, 'degF')),     (26092.78, 'kelvin**2')),
+
+        (((100, 'degC'), (10, 'kelvin')),     (3731.5, 'kelvin**2')),
+        (((100, 'degC'), (10, 'degC')),       (105657.42, 'kelvin**2')),
+        (((100, 'degC'), (10, 'degF')),       (97365.20, 'kelvin**2')),
+        (((100, 'degC'), (10, 'degR')),       (3731.5, 'kelvin*degR')),
+        (((100, 'degC'), (10, 'delta_degC')), (3731.5, 'kelvin*delta_degC')),
+        (((100, 'degC'), (10, 'delta_degF')), (3731.5, 'kelvin*delta_degF')),
+
+        (((100, 'degF'), (10, 'kelvin')),     (3109.28, 'kelvin**2')),
+        (((100, 'degF'), (10, 'degC')),       (88039.20, 'kelvin**2')),
+        (((100, 'degF'), (10, 'degF')),       (81129.69, 'kelvin**2')),
+        (((100, 'degF'), (10, 'degR')),       (3109.28, 'kelvin*degR')),
+        (((100, 'degF'), (10, 'delta_degC')), (3109.28, 'kelvin*delta_degC')),
+        (((100, 'degF'), (10, 'delta_degF')), (3109.28, 'kelvin*delta_degF')),
+
+        (((100, 'degR'), (10, 'degC')),       (28315., 'degR*kelvin')),
+        (((100, 'degR'), (10, 'degF')),       (26092.78, 'degR*kelvin')),
+
+        (((100, 'delta_degC'), (10, 'degC')), (28315., 'delta_degC*kelvin')),
+        (((100, 'delta_degC'), (10, 'degF')), (26092.78, 'delta_degC*kelvin')),
+
+        (((100, 'delta_degF'), (10, 'degC')), (28315., 'delta_degF*kelvin')),
+        (((100, 'delta_degF'), (10, 'degF')), (26092.78, 'delta_degF*kelvin')),
+        ]
+
+    @ParameterizedTestCase.parameterize(
+        ("input", "expected_output"),
+        multiplications_with_autoconvert_to_baseunit)
+    def test_multiplication_with_autoconvert(self, input_tuple, expected):
+        self.ureg.autoconvert_offset_to_baseunit = True
+        qin1, qin2 = input_tuple
+        q1, q2 = self.Q_(*qin1), self.Q_(*qin2)
+        input_tuple = q1, q2
+        if expected == 'error':
+            self.assertRaises(OffsetUnitCalculusError, op.mul, q1, q2)
+        else:
+            expected = self.Q_(*expected)
+            self.assertEqual(op.mul(q1, q2).units, expected.units)
+            self.assertQuantityAlmostEqual(op.mul(q1, q2), expected,
+                                           atol=0.01)
+
+    @helpers.requires_numpy()
+    @ParameterizedTestCase.parameterize(
+        ("input", "expected_output"),
+        multiplications_with_autoconvert_to_baseunit)
+    def test_inplace_multiplication_with_autoconvert(self, input_tuple, expected):
+        self.ureg.autoconvert_offset_to_baseunit = True
+        (q1v, q1u), (q2v, q2u) = input_tuple
+        # update input tuple with new values to have correct values on failure
+        input_tuple = ((np.array([q1v]*2, dtype=np.float), q1u),
+                       (np.array([q2v]*2, dtype=np.float), q2u))
+        Q_ = self.Q_
+        qin1, qin2 = input_tuple
+        q1, q2 = Q_(*qin1), Q_(*qin2)
+        q1_cp = copy.copy(q1)
+        if expected == 'error':
+            self.assertRaises(OffsetUnitCalculusError, op.imul, q1_cp, q2)
+        else:
+            expected = np.array([expected[0]]*2, dtype=np.float), expected[1]
+            self.assertEqual(op.imul(q1_cp, q2).units, Q_(*expected).units)
+            q1_cp = copy.copy(q1)
+            self.assertQuantityAlmostEqual(op.imul(q1_cp, q2), Q_(*expected),
+                                           atol=0.01)
+
+    multiplications_with_scalar = [
+        (((10, 'kelvin'),    2),    (20., 'kelvin')),
+        (((10, 'kelvin**2'), 2),    (20., 'kelvin**2')),
+        (((10, 'degC'),      2),    (20., 'degC')),
+        (((10, '1/degC'),    2),    'error'),
+        (((10, 'degC**0.5'),   2),  'error'),
+        (((10, 'degC**2'),   2),    'error'),
+        (((10, 'degC**-2'),  2),    'error'),
+        ]
+
+    @ParameterizedTestCase.parameterize(
+        ("input", "expected_output"), multiplications_with_scalar)
+    def test_multiplication_with_scalar(self, input_tuple, expected):
+        self.ureg.default_as_delta = False
+        in1, in2 = input_tuple
+        if type(in1) is tuple:
+            in1, in2 = self.Q_(*in1), in2
+        else:
+            in1, in2 = in1, self.Q_(*in2)
+        input_tuple = in1, in2  # update input_tuple for better tracebacks
+        if expected == 'error':
+            self.assertRaises(OffsetUnitCalculusError, op.mul, in1, in2)
+        else:
+            expected = self.Q_(*expected)
+            self.assertEqual(op.mul(in1, in2).units, expected.units)
+            self.assertQuantityAlmostEqual(op.mul(in1, in2), expected,
+                                           atol=0.01)
+
+    divisions_with_scalar = [   # without / with autoconvert to base unit
+        (((10, 'kelvin'), 2),       [(5., 'kelvin'), (5., 'kelvin')]),
+        (((10, 'kelvin**2'), 2),    [(5., 'kelvin**2'), (5., 'kelvin**2')]),
+        (((10, 'degC'), 2),         ['error', 'error']),
+        (((10, 'degC**2'), 2),      ['error', 'error']),
+        (((10, 'degC**-2'), 2),     ['error', 'error']),
+
+        ((2, (10, 'kelvin')),       [(0.2, '1/kelvin'), (0.2, '1/kelvin')]),
+        ((2, (10, 'degC')),         ['error', (2/283.15, '1/kelvin')]),
+        ((2, (10, 'degC**2')),      ['error', 'error']),
+        ((2, (10, 'degC**-2')),     ['error', 'error']),
+        ]
+
+    @ParameterizedTestCase.parameterize(
+        ("input", "expected_output"), divisions_with_scalar)
+    def test_division_with_scalar(self, input_tuple, expected):
+        self.ureg.default_as_delta = False
+        in1, in2 = input_tuple
+        if type(in1) is tuple:
+            in1, in2 = self.Q_(*in1), in2
+        else:
+            in1, in2 = in1, self.Q_(*in2)
+        input_tuple = in1, in2  # update input_tuple for better tracebacks
+        expected_copy = expected[:]
+        for i, mode in enumerate([False, True]):
+            self.ureg.autoconvert_offset_to_baseunit = mode
+            if expected_copy[i] == 'error':
+                self.assertRaises(OffsetUnitCalculusError, op.truediv, in1, in2)
+            else:
+                expected = self.Q_(*expected_copy[i])
+                self.assertEqual(op.truediv(in1, in2).units, expected.units)
+                self.assertQuantityAlmostEqual(op.truediv(in1, in2), expected)
+
+    exponentiation = [                  # resuls without / with autoconvert
+        (((10, 'degC'),    1),          [(10, 'degC'), (10, 'degC')]),
+        (((10, 'degC'),    0.5),        ['error', (283.15**0.5, 'kelvin**0.5')]),
+        (((10, 'degC'),    0),          [(1., ''), (1., '')]),
+        (((10, 'degC'),   -1),          ['error', (1/(10+273.15), 'kelvin**-1')]),
+        (((10, 'degC'),   -2),          ['error', (1/(10+273.15)**2., 'kelvin**-2')]),
+        ((( 0, 'degC'),   -2),          ['error', (1/(273.15)**2, 'kelvin**-2')]),
+        (((10, 'degC'),   (2, '')),     ['error', ((283.15)**2, 'kelvin**2')]),
+        (((10, 'degC'),  (10, 'degK')), ['error', 'error']),
+
+        (((10, 'kelvin'), (2, '')),     [(100., 'kelvin**2'), (100., 'kelvin**2')]),
+
+        ((  2,          (2, 'kelvin')), ['error', 'error']),
+        ((  2,          (500., 'millikelvin/kelvin')), [2**0.5, 2**0.5]),
+        ((  2,          (0.5, 'kelvin/kelvin')),      [2**0.5, 2**0.5]),
+        (((10, 'degC'), (500., 'millikelvin/kelvin')),
+                                        ['error', (283.15**0.5, 'kelvin**0.5')]),
+         ]
+
+    @ParameterizedTestCase.parameterize(
+        ("input", "expected_output"), exponentiation)
+    def test_exponentiation(self, input_tuple, expected):
+        self.ureg.default_as_delta = False
+        in1, in2 = input_tuple
+        if type(in1) is tuple and type(in2) is tuple:
+            in1, in2 = self.Q_(*in1), self.Q_(*in2)
+        elif not type(in1) is tuple and type(in2) is tuple:
+            in2 = self.Q_(*in2)
+        else:
+            in1 = self.Q_(*in1)
+        input_tuple = in1, in2
+        expected_copy = expected[:]
+        for i, mode in enumerate([False, True]):
+            self.ureg.autoconvert_offset_to_baseunit = mode
+            if expected_copy[i] == 'error':
+                self.assertRaises((OffsetUnitCalculusError,
+                                   DimensionalityError), op.pow, in1, in2)
+            else:
+                if type(expected_copy[i]) is tuple:
+                    expected = self.Q_(*expected_copy[i])
+                    self.assertEqual(op.pow(in1, in2).units, expected.units)
+                else:
+                    expected = expected_copy[i]
+                self.assertQuantityAlmostEqual(op.pow(in1, in2), expected)
+
+    @helpers.requires_numpy()
+    @ParameterizedTestCase.parameterize(
+        ("input", "expected_output"), exponentiation)
+    def test_inplace_exponentiation(self, input_tuple, expected):
+        self.ureg.default_as_delta = False
+        in1, in2 = input_tuple
+        if type(in1) is tuple and type(in2) is tuple:
+            (q1v, q1u), (q2v, q2u) = in1, in2
+            in1 = self.Q_(*(np.array([q1v]*2, dtype=np.float), q1u))
+            in2 = self.Q_(q2v, q2u)
+        elif not type(in1) is tuple and type(in2) is tuple:
+            in2 = self.Q_(*in2)
+        else:
+            in1 = self.Q_(*in1)
+
+        input_tuple = in1, in2
+
+        expected_copy = expected[:]
+        for i, mode in enumerate([False, True]):
+            self.ureg.autoconvert_offset_to_baseunit = mode
+            in1_cp = copy.copy(in1)
+            if expected_copy[i] == 'error':
+                self.assertRaises((OffsetUnitCalculusError,
+                                   DimensionalityError), op.ipow, in1_cp, in2)
+            else:
+                if type(expected_copy[i]) is tuple:
+                    expected = self.Q_(np.array([expected_copy[i][0]]*2,
+                                                dtype=np.float),
+                                       expected_copy[i][1])
+                    self.assertEqual(op.ipow(in1_cp, in2).units, expected.units)
+                else:
+                    expected = np.array([expected_copy[i]]*2, dtype=np.float)
+
+
+                in1_cp = copy.copy(in1)
+                self.assertQuantityAlmostEqual(op.ipow(in1_cp, in2), expected)
diff --git a/pint/testsuite/test_unit.py b/pint/testsuite/test_unit.py
index 7131db8..c0f4602 100644
--- a/pint/testsuite/test_unit.py
+++ b/pint/testsuite/test_unit.py
@@ -12,8 +12,9 @@ from pint.unit import (ScaleConverter, OffsetConverter, UnitsContainer,
                        DimensionDefinition, _freeze, Converter, UnitRegistry,
                        LazyRegistry, ParserHelper)
 from pint import DimensionalityError, UndefinedUnitError
-from pint.compat import u, unittest, np
+from pint.compat import u, unittest, np, string_types
 from pint.testsuite import QuantityTestCase, helpers, BaseTestCase
+from pint.testsuite.parameterized import ParameterizedTestCase
 
 
 class TestConverter(BaseTestCase):
@@ -231,6 +232,9 @@ class TestRegistry(QuantityTestCase):
 
     FORCE_NDARRAY = False
 
+    def setup(self):
+        self.ureg.autoconvert_offset_to_baseunit = False
+
     def test_base(self):
         ureg = UnitRegistry(None)
         ureg.define('meter = [length]')
@@ -357,18 +361,6 @@ class TestRegistry(QuantityTestCase):
         self.assertEqual(self.ureg.get_symbol('international_foot'), 'ft')
         self.assertEqual(self.ureg.get_symbol('international_inch'), 'in')
 
-    @unittest.expectedFailure
-    def test_delta_in_diff(self):
-        """This might be supported in future versions
-        """
-        xk = 1 * self.ureg.kelvin
-        yk = 2 * self.ureg.kelvin
-        yf = yk.to('degF')
-        yc = yk.to('degC')
-        self.assertEqual(yk - xk, 1 * self.ureg.kelvin)
-        self.assertEqual(yf - xk, 1 * self.ureg.kelvin)
-        self.assertEqual(yc - xk, 1 * self.ureg.kelvin)
-
     def test_pint(self):
         p = self.ureg.pint
         l = self.ureg.liter
@@ -440,12 +432,13 @@ class TestRegistry(QuantityTestCase):
         self.assertEqual(h2(3, 1), (3 * ureg.meter, 1 * ureg.cm))
 
     def test_to_ref_vs_to(self):
+        self.ureg.autoconvert_offset_to_baseunit = True
         q = 8. * self.ureg.inch
         t = 8. * self.ureg.degF
         dt = 8. * self.ureg.delta_degF
         self.assertEqual(q.to('cm').magnitude, self.ureg._units['inch'].converter.to_reference(8.))
         self.assertEqual(t.to('kelvin').magnitude, self.ureg._units['degF'].converter.to_reference(8.))
-        self.assertEqual(dt.to('delta_kelvin').magnitude, self.ureg._units['delta_degF'].converter.to_reference(8.))
+        self.assertEqual(dt.to('kelvin').magnitude, self.ureg._units['delta_degF'].converter.to_reference(8.))
 
     def test_redefinition(self):
         d = UnitRegistry().define
@@ -473,7 +466,7 @@ class TestRegistry(QuantityTestCase):
         # Conversions with single units take a different codepath than
         # Conversions with more than one unit.
         src_dst1 = UnitsContainer(meter=1), UnitsContainer(inch=1)
-        src_dst2 = UnitsContainer(meter=1, seconds=-1), UnitsContainer(inch=1, minutes=-1)
+        src_dst2 = UnitsContainer(meter=1, second=-1), UnitsContainer(inch=1, minute=-1)
         for src, dst in (src_dst1, src_dst2):
             v = ureg.convert(1, src, dst),
 
@@ -598,3 +591,52 @@ class TestErrors(BaseTestCase):
         msg = "Cannot convert from 'a' (c) to 'b' (d)msg"
         ex = DimensionalityError('a', 'b', 'c', 'd', 'msg')
         self.assertEqual(str(ex), msg)
+
+
+class TestConvertWithOffset(QuantityTestCase, ParameterizedTestCase):
+
+    # The dicts in convert_with_offset are used to create a UnitsContainer.
+    # We create UnitsContainer to avoid any auto-conversion of units.
+    convert_with_offset = [
+        (({'degC': 1}, {'degC': 1}), 10),
+        (({'degC': 1}, {'kelvin': 1}), 283.15),
+        (({'degC': 1}, {'degC': 1, 'millimeter': 1, 'meter': -1}), 'error'),
+        (({'degC': 1}, {'kelvin': 1, 'millimeter': 1, 'meter': -1}), 283150),
+
+        (({'kelvin': 1}, {'degC': 1}), -263.15),
+        (({'kelvin': 1}, {'kelvin': 1}), 10),
+        (({'kelvin': 1}, {'degC': 1, 'millimeter': 1, 'meter': -1}), 'error'),
+        (({'kelvin': 1}, {'kelvin': 1, 'millimeter': 1, 'meter': -1}), 10000),
+
+        (({'degC': 1, 'millimeter': 1, 'meter': -1}, {'degC': 1}), 'error'),
+        (({'degC': 1, 'millimeter': 1, 'meter': -1}, {'kelvin': 1}), 'error'),
+        (({'degC': 1, 'millimeter': 1, 'meter': -1}, {'degC': 1, 'millimeter': 1, 'meter': -1}), 10),
+        (({'degC': 1, 'millimeter': 1, 'meter': -1}, {'kelvin': 1, 'millimeter': 1, 'meter': -1}), 'error'),
+
+        (({'kelvin': 1, 'millimeter': 1, 'meter': -1}, {'degC': 1}), -273.14),
+        (({'kelvin': 1, 'millimeter': 1, 'meter': -1}, {'kelvin': 1}), 0.01),
+        (({'kelvin': 1, 'millimeter': 1, 'meter': -1}, {'degC': 1, 'millimeter': 1, 'meter': -1}), 'error'),
+        (({'kelvin': 1, 'millimeter': 1, 'meter': -1}, {'kelvin': 1, 'millimeter': 1, 'meter': -1}), 10),
+
+        (({'degC': 2}, {'kelvin': 2}), 'error'),
+        (({'degC': 1, 'degF': 1}, {'kelvin': 2}), 'error'),
+        (({'degC': 1, 'kelvin': 1}, {'kelvin': 2}), 'error'),
+        ]
+
+    @ParameterizedTestCase.parameterize(("input", "expected_output"),
+                                        convert_with_offset)
+    def test_to_and_from_offset_units(self, input_tuple, expected):
+        src, dst = input_tuple
+        src, dst = UnitsContainer(src), UnitsContainer(dst)
+        value = 10.
+        convert = self.ureg.convert
+        if isinstance(expected, string_types):
+            self.assertRaises(DimensionalityError, convert, value, src, dst)
+            if src != dst:
+                self.assertRaises(DimensionalityError, convert, value, dst, src)
+        else:
+            self.assertQuantityAlmostEqual(convert(value, src, dst),
+                                           expected, atol=0.001)
+            if src != dst:
+                self.assertQuantityAlmostEqual(convert(expected, dst, src),
+                                               value, atol=0.001)
diff --git a/pint/unit.py b/pint/unit.py
index 37cd63b..27f6b5e 100644
--- a/pint/unit.py
+++ b/pint/unit.py
@@ -87,7 +87,7 @@ class DimensionalityError(ValueError):
     """Raised when trying to convert between incompatible units.
     """
 
-    def __init__(self, units1, units2, dim1=None, dim2=None, extra_msg=None):
+    def __init__(self, units1, units2, dim1=None, dim2=None, extra_msg=''):
         super(DimensionalityError, self).__init__()
         self.units1 = units1
         self.units2 = units2
@@ -103,12 +103,25 @@ class DimensionalityError(ValueError):
             dim1 = ''
             dim2 = ''
 
-        msg = "Cannot convert from '{0}'{1} to '{2}'{3}".format(self.units1, dim1, self.units2, dim2)
+        msg = "Cannot convert from '{0}'{1} to '{2}'{3}" + self.extra_msg
 
-        if self.extra_msg:
-            return msg + self.extra_msg
+        return msg.format(self.units1, dim1, self.units2, dim2)
 
-        return msg
+
+class OffsetUnitCalculusError(ValueError):
+    """Raised on ambiguous operations with offset units.
+    """
+    def __init__(self, units1, units2='', extra_msg=''):
+        super(ValueError, self).__init__()
+        self.units1 = units1
+        self.units2 = units2
+        self.extra_msg = extra_msg
+
+    def __str__(self):
+        msg = ("Ambiguous operation with offset unit (%s)." %
+               ', '.join(['%s' % u for u in [self.units1, self.units2] if u])
+               + self.extra_msg)
+        return msg.format(self.units1, self.units2)
 
 
 class Converter(object):
@@ -289,7 +302,7 @@ class UnitDefinition(Definition):
                                  'Base units must be referenced only to dimensions. '
                                  'Derived units must be referenced only to units.')
             self.reference = UnitsContainer(converter.items())
-            if 'offset' in modifiers:
+            if modifiers.get('offset', 0.) != 0.:
                 converter = OffsetConverter(converter.scale, modifiers['offset'])
             else:
                 converter = ScaleConverter(converter.scale)
@@ -447,12 +460,17 @@ class UnitRegistry(object):
     :param force_ndarray: convert any input, scalar or not to a numpy.ndarray.
     :param default_as_delta: In the context of a multiplication of units, interpret
                              non-multiplicative units as their *delta* counterparts.
+    :autoconvert_offset_to_baseunit: If True converts offset units in quantites are
+                                     converted to their base units in multiplicative
+                                     context. If False no conversion happens.
     :param on_redefinition: action to take in case a unit is redefined.
                             'warn', 'raise', 'ignore'
     :type on_redefintion: str
     """
 
-    def __init__(self, filename='', force_ndarray=False, default_as_delta=True, on_redefinition='warn'):
+    def __init__(self, filename='', force_ndarray=False, default_as_delta=True,
+                 autoconvert_offset_to_baseunit=False,
+                 on_redefinition='warn'):
         self.Quantity = build_quantity_class(self, force_ndarray)
         self.Measurement = build_measurement_class(self, force_ndarray)
 
@@ -495,6 +513,10 @@ class UnitRegistry(object):
         #: non-multiplicative units as their *delta* counterparts.
         self.default_as_delta = default_as_delta
 
+        # Determines if quantities with offset units are converted to their
+        # base units on multiplication and division.
+        self.autoconvert_offset_to_baseunit = autoconvert_offset_to_baseunit
+
         if filename == '':
             self.load_definitions('default_en.txt', True)
         elif filename is not None:
@@ -697,7 +719,8 @@ class UnitRegistry(object):
 
             _adder(alias, definition)
 
-        if isinstance(definition.converter, OffsetConverter):
+        # define additional "delta_" units for units with an offset
+        if getattr(definition.converter, "offset", 0.0) != 0.0:
             d_name = 'delta_' + definition.name
             if definition.symbol:
                 d_symbol = 'Δ' + definition.symbol
@@ -710,7 +733,7 @@ class UnitRegistry(object):
                     return '[delta_' + _name[1:]
                 return 'delta_' + _name
 
-            d_reference = UnitsContainer(dict((prep(ref), value)
+            d_reference = UnitsContainer(dict((ref, value)
                                          for ref, value in definition.reference.items()))
 
             self.define(UnitDefinition(d_name, d_symbol, d_aliases,
@@ -966,9 +989,9 @@ class UnitRegistry(object):
         :return: converted value
         """
         if isinstance(src, string_types):
-            src = ParserHelper.from_string(src)
+            src = self.parse_units(src)
         if isinstance(dst, string_types):
-            dst = ParserHelper.from_string(dst)
+            dst = self.parse_units(dst)
         if src == dst:
             return value
 
@@ -996,36 +1019,59 @@ class UnitRegistry(object):
         if src_dim != dst_dim:
             raise DimensionalityError(src, dst, src_dim, dst_dim)
 
-        if len(src) == 1:
-            # If the source has a single element, it might be a non-multiplicative unit
-            # and therefore it is treated differently.
-            src_unit, src_value = list(src.items())[0]
-            src_unit = self._units[self.get_name(src_unit)]
-
-            # We only continue if is a ScaleConverter,
-            # if not just exit to use the standard src / dst.
-            # TODO: This will fail and should not degK * meter / nanometer -> degC
-            if not isinstance(src_unit.converter, ScaleConverter):
+        # Conversion needs to consider if non-multiplicative (AKA offset
+        # units) are involved. Conversion is only possible if src and dst
+        # have at most one offset unit per dimension.
+        src_offset_units = [(u, e) for u, e in src.items()
+                            if not self._units[u].is_multiplicative]
+        dst_offset_units = [(u, e) for u, e in dst.items()
+                            if not self._units[u].is_multiplicative]
+
+        # For offset units we need to check if the conversion is allowed.
+        if src_offset_units or dst_offset_units:
+
+            # Validate that not more than one offset unit is present
+            if len(src_offset_units) > 1 or len(dst_offset_units) > 1:
+                raise DimensionalityError(
+                    src, dst, src_dim, dst_dim,
+                    extra_msg=' - more than one offset unit.')
+
+            # validate that offset unit is not used in multiplicative context
+            if ((len(src_offset_units) == 1 and len(src) > 1)
+                    or (len(dst_offset_units) == 1 and len(dst) > 1)
+                    and not self.autoconvert_offset_to_baseunit):
+                raise DimensionalityError(
+                    src, dst, src_dim, dst_dim,
+                    extra_msg=' - offset unit used in multiplicative context.')
+
+            # Validate that order of offset unit is exactly one.
+            if src_offset_units:
+                if src_offset_units[0][1] != 1:
+                    raise DimensionalityError(
+                        src, dst, src_dim, dst_dim,
+                        extra_msg=' - offset units in higher order.')
+            else:
+                if dst_offset_units[0][1] != 1:
+                    raise DimensionalityError(
+                        src, dst, src_dim, dst_dim,
+                        extra_msg=' - offset units in higher order.')
 
-                if len(dst) > 1:
-                    # If the destination has more than one element,
-                    # then the conversion is not possible.
-                    # TODO: This will fail and should not degC -> degK * meter / nanometer
-                    raise DimensionalityError(src, dst, src_dim, dst_dim)
+        # Here we convert only the offset quantities. Any remaining scaled
+        # quantities will be converted later.
 
-                dst_unit, dst_value = list(dst.items())[0]
-                dst_unit = self._units[self.get_name(dst_unit)]
-                if not type(src_unit.converter) is type(dst_unit.converter):
-                    raise DimensionalityError(src, dst, src_dim, dst_dim)
+        # clean src from offset units by converting to reference
+        for u, e in src_offset_units:
+            value = self._units[u].converter.to_reference(value, inplace)
+            src.pop(u)
 
-                return dst_unit.converter.from_reference(src_unit.converter.to_reference(value, inplace), inplace)
+        # clean dst units from offset units
+        for u, e in dst_offset_units:
+            dst.pop(u)
 
+        # Here src and dst have only multiplicative units left. Thus we can
+        # convert with a factor.
         factor, units = self.get_base_units(src / dst)
 
-        if factor is None:
-            raise DimensionalityError(src, dst, src_dim, dst_dim,
-                                      'Non-multiplicative unit found.')
-
         # factor is type float and if our magnitude is type Decimal then
         # must first convert to Decimal before we can '*' the values
         if isinstance(value, Decimal):
@@ -1036,6 +1082,16 @@ class UnitRegistry(object):
         else:
             value = value * factor
 
+        # Finally convert to offset units specified in destination
+        for u, e in dst_offset_units:
+            value = self._units[u].converter.from_reference(value, inplace)
+            # add back offset units to dst
+            dst[u] = e
+
+        # restore offset conversion of src units
+        for u, e in src_offset_units:
+            src[u] = e
+
         return value
 
     def pi_theorem(self, quantities):
@@ -1116,7 +1172,7 @@ class UnitRegistry(object):
             cname = self.get_name(name)
             if not cname:
                 continue
-            if as_delta and (many or (not many and abs(value) != 1)):
+            if as_delta and (many or (not many and value != 1)):
                 definition = self._units[cname]
                 if not definition.is_multiplicative:
                     cname = 'delta_' + cname
@@ -1140,7 +1196,7 @@ class UnitRegistry(object):
         unknown = set()
         for toknum, tokval, _, _, _ in gen:
             if toknum == NAME:
-                # TODO: Integrate math better, Replace eval
+                # TODO: Integrate math better, Replace eval, make as_delta-aware
                 if tokval == 'pi' or tokval in values:
                     result.append((toknum, tokval))
                     continue