astropy
diff --git a/‎CHANGES.rst
Lines changed: 5 additions & 0 deletions b/‎CHANGES.rst
Lines changed: 5 additions & 0 deletions
diff --git a/‎astropy/coordinates/angles.py
Lines changed: 20 additions & 8 deletions b/‎astropy/coordinates/angles.py
Lines changed: 20 additions & 8 deletions
diff --git a/‎astropy/modeling/functional_models.py
Lines changed: 7 additions & 5 deletions b/‎astropy/modeling/functional_models.py
Lines changed: 7 additions & 5 deletions
@@ -187,6 +187,11 @@ New Features
   - ``Quantity`` now supports the ``@`` operator for matrix multiplication that
     was introduced in Python 3.5, for all supported versions of numpy. [#6144]
 
+  - ``Quantity`` supports the new ``__array_ufunc__`` protocol introduced in
+    numpy 1.13.  As a result, operations that involve unit conversion will be
+    sped up considerably (by up to a factor of two for costly operations such
+    as trigonometric ones). [#2583]
+
 - ``astropy.utils``
 
   - Added a new ``dataurl_mirror`` configuration item in ``astropy.utils.data``
 
@@ -444,6 +444,18 @@ def _repr_latex_(self):
                                    formatter={'str_kind': lambda x: x})
 
 
+def _no_angle_subclass(obj):
+    """Return any Angle subclass objects as an Angle objects.
+
+    This is used to ensure that Latitute and Longitude change to Angle
+    objects when they are used in calculations (such as lon/2.)
+    """
+    if isinstance(obj, tuple):
+        return tuple(_no_angle_subclass(_obj) for _obj in obj)
+
+    return obj.view(Angle) if isinstance(obj, Angle) else obj
+
+
 class Latitude(Angle):
     """
     Latitude-like angle(s) which must be in the range -90 to +90 deg.
@@ -525,11 +537,11 @@ def __setitem__(self, item, value):
     # Any calculation should drop to Angle
     def __array_wrap__(self, obj, context=None):
         obj = super(Angle, self).__array_wrap__(obj, context=context)
+        return _no_angle_subclass(obj)
 
-        if isinstance(obj, Angle):
-            return obj.view(Angle)
-
-        return obj
+    def __array_ufunc__(self, *args, **kwargs):
+        results = super(Latitude, self).__array_ufunc__(*args, **kwargs)
+        return _no_angle_subclass(results)
 
 
 class Longitude(Angle):
@@ -644,8 +656,8 @@ def __array_finalize__(self, obj):
     # Any calculation should drop to Angle
     def __array_wrap__(self, obj, context=None):
         obj = super(Angle, self).__array_wrap__(obj, context=context)
+        return _no_angle_subclass(obj)
 
-        if isinstance(obj, Angle):
-            return obj.view(Angle)
-
-        return obj
+    def __array_ufunc__(self, *args, **kwargs):
+        results = super(Longitude, self).__array_ufunc__(*args, **kwargs)
+        return _no_angle_subclass(results)
@@ -2185,19 +2185,21 @@ def evaluate(cls, x, y, amplitude, x_0, y_0, radius):
         r = np.sqrt((x - x_0) ** 2 + (y - y_0) ** 2) / (radius / cls._rz)
 
         if isinstance(r, Quantity):
-            r = r.decompose().value
+            # scipy function cannot handle Quantity, so turn into array.
+            r = r.to_value(u.dimensionless_unscaled)
 
         # Since r can be zero, we have to take care to treat that case
         # separately so as not to raise a numpy warning
         z = np.ones(r.shape)
         rt = np.pi * r[r > 0]
         z[r > 0] = (2.0 * cls._j1(rt) / rt) ** 2
-        z *= amplitude
 
         if isinstance(amplitude, Quantity):
-            return Quantity(z, unit=amplitude.unit, copy=False)
-        else:
-            return z
+            # make z quantity too, otherwise in-place multiplication fails.
+            z = Quantity(z, u.dimensionless_unscaled, copy=False)
+
+        z *= amplitude
+        return z
 
     @property
     def input_units(self):