Add Beta1D model

astropy/modeling/functional_models.py

@@ -18,10 +18,10 @@
 from ..units import Quantity, UnitsError
 from ..utils.decorators import deprecated
 
-__all__ = ['AiryDisk2D', 'Moffat1D', 'Moffat2D', 'Box1D', 'Box2D', 'Const1D',
-           'Const2D', 'Ellipse2D', 'Disk2D', 'BaseGaussian1D', 'Gaussian1D',
-           'GaussianAbsorption1D', 'Gaussian2D', 'Linear1D', 'Lorentz1D',
-           'MexicanHat1D', 'MexicanHat2D', 'RedshiftScaleFactor',
+__all__ = ['AiryDisk2D', 'Moffat1D', 'Moffat2D', 'Beta1D', 'Box1D', 'Box2D',
+           'Const1D', 'Const2D', 'Ellipse2D', 'Disk2D', 'BaseGaussian1D',
+           'Gaussian1D', 'GaussianAbsorption1D', 'Gaussian2D', 'Linear1D',
+           'Lorentz1D', 'MexicanHat1D', 'MexicanHat2D', 'RedshiftScaleFactor',
            'Scale', 'Sersic1D', 'Sersic2D', 'Shift', 'Sine1D', 'Trapezoid1D',
            'TrapezoidDisk2D', 'Ring2D', 'Voigt1D']
 
@@ -862,6 +862,61 @@ def _parameter_units_for_data_units(self, inputs_unit, outputs_unit):
                             ('amplitude', outputs_unit['y'])])
 
 
+class Beta1D(Fittable1DModel):
+    """
+    1-D Lorentz model with a varying power law, a.k.a, beta model.
+
+    Parameters
+    ----------
+    amplitude : float
+        Amplitude at x=xpos
+    beta : float
+        Beta index
+    r0 : float
+        Core radius
+    xpos : float
+        Offset from x=0
+
+    Notes
+    -----
+    Model formula:
+
+        .. math:: f(x) = a * (1 + [(x - x_{pos})/r_0] ^ 2)^(-3 * \\beta + 1/2)
+
+    See Also
+    --------
+    Lorentz1D
+    """
+
+    amplitude = Parameter(default=1)
+    beta = Parameter(default=1)
+    r0 = Parameter(default=1)
+    xpos = Parameter(default=0)
+
+    @staticmethod
+    def evaluate(x, amplitude, beta, r0, xpos):
+        return amplitude * (1 + ((x - xpos) / r0) ** 2) ** (-3 * beta + .5)
+
+    @staticmethod
+    def fit_deriv(x, amplitude, beta, r0, xpos):
+        d_r0 = (-2 * amplitude * ((x - xpos) / r0) ** 2 / r0 * (-3 * beta + .5)
+                * (1 + ((x - xpos) / r0) ** 2) ** (-3 * beta - .5))
+        d_beta = (- 3 * amplitude * (1 + ((x - xpos) / r0) ** 2) ** (-3 * beta
+                  + .5) * np.log(1 + ((x - xpos) / r0) ** 2))
+        d_xpos = (-2 * amplitude * (-3 * beta + .5) * (1 + ((x - xpos) / r0)
+                  ** 2) ** (-3 * beta - .5) * (x - xpos) / r0 ** 2)
+        d_amplitude = (1 + ((x - xpos) / r0) ** 2) ** (-3 * beta + .5)
+
+        return [d_amplitude, d_beta, d_r0, d_xpos]
+
+    @property
+    def input_units(self):
+        if self.xpos.unit is None:
+            return None
+        else:
+            return {'x': self.xpos.unit}
+
+
 class Linear1D(Fittable1DModel):
     """
     One dimensional Line model.

astropy/modeling/tests/test_functional_models.py

@@ -20,6 +20,33 @@
     HAS_SCIPY = False
 
 
+def test_Beta1D():
+    model = models.Beta1D()
+    x = np.linspace(0, 4, 10)
+    y = model(x)
+    assert y[0] == model.amplitude.value
+    assert (y[:len(y) - 1] > y[1:]).all()
+    assert_allclose(0, model(20), atol=1e-6)
+
+    model = models.Beta1D(beta=1/6)
+    y = model(x)
+    assert (y == model.amplitude.value).all()
+
+    model = models.Beta1D(r0=1e20)
+    y = model(x)
+    assert (y == model.amplitude.value).all()
+
+    model = models.Beta1D(beta=1e20)
+    y = model(x)
+    assert y[0] == model.amplitude.value
+    assert (y[1:] == 0).all()
+
+    model = models.Beta1D(beta=-1e10)
+    y = model(x)
+    assert y[0] == model.amplitude.value
+    assert (y[1:] == np.inf).all()
+
+
 def test_Trapezoid1D():
     """Regression test for https://github.com/astropy/astropy/issues/1721"""