Add Beta1D model

CHANGES.rst

@@ -40,6 +40,8 @@ astropy.io.votable
 astropy.modeling
 ^^^^^^^^^^^^^^^^
 
+- Add Beta1D to models. [#6374]
+
 astropy.nddata
 ^^^^^^^^^^^^^^
 

astropy/modeling/functional_models.py

@@ -21,8 +21,8 @@
 __all__ = ['AiryDisk2D', 'Moffat1D', 'Moffat2D', 'Box1D', 'Box2D', 'Const1D',
            'Const2D', 'Ellipse2D', 'Disk2D', 'BaseGaussian1D', 'Gaussian1D',
            'GaussianAbsorption1D', 'Gaussian2D', 'Linear1D', 'Lorentz1D',
-           'MexicanHat1D', 'MexicanHat2D', 'RedshiftScaleFactor',
-           'Scale', 'Sersic1D', 'Sersic2D', 'Shift', 'Sine1D', 'Trapezoid1D',
+           'MexicanHat1D', 'MexicanHat2D', 'RedshiftScaleFactor', 'Scale',
+           'Sersic1D', 'Sersic2D', 'Shift', 'Sine1D', 'Trapezoid1D',
            'TrapezoidDisk2D', 'Ring2D', 'Voigt1D']
 
 TWOPI = 2 * np.pi

astropy/modeling/powerlaws.py

@@ -15,8 +15,65 @@
 from .parameters import Parameter, InputParameterError
 from ..units import Quantity
 
-__all__ = ['PowerLaw1D', 'BrokenPowerLaw1D', 'SmoothlyBrokenPowerLaw1D',
-           'ExponentialCutoffPowerLaw1D', 'LogParabola1D']
+__all__ = ['PowerLaw1D', 'Beta1D', 'BrokenPowerLaw1D',
+           'SmoothlyBrokenPowerLaw1D', 'ExponentialCutoffPowerLaw1D',
+           'LogParabola1D']
+
+
+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):
+        aux = ((x - xpos) / r0) ** 2 + 1
+
+        d_r0 = (-2 * amplitude * (aux - 1) / r0 * (-3 * beta + .5)
+                * aux ** (-3 * beta - .5))
+        d_beta = (- 3 * amplitude * aux ** (-3 * beta + .5) * np.log(aux))
+        d_xpos = (-2 * amplitude * (-3 * beta + .5) * aux ** (-3 * beta - .5)
+                  * (x - xpos) / r0 ** 2)
+        d_amplitude = aux ** (-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 PowerLaw1D(Fittable1DModel):

astropy/modeling/tests/example_models.py

@@ -60,7 +60,7 @@
     Ring2D, Sersic1D, Sersic2D, Voigt1D, Planar2D)
 from ..polynomial import Polynomial1D, Polynomial2D
 from ..powerlaws import (
-    PowerLaw1D, BrokenPowerLaw1D, SmoothlyBrokenPowerLaw1D, ExponentialCutoffPowerLaw1D,
+    PowerLaw1D, BrokenPowerLaw1D, Beta1D, SmoothlyBrokenPowerLaw1D, ExponentialCutoffPowerLaw1D,
     LogParabola1D)
 import numpy as np
 
@@ -90,6 +90,13 @@
         'integral': 10
     },
 
+    Beta1D: {
+        'parameters': [1, 1, 1, 0],
+        'x_values': [0, 0.25, 0.5, 1.0],
+        'y_values': [1., 0.85936498, 0.5724334, 0.1767767],
+        'x_lim': [-10, 10]
+    },
+
     Linear1D: {
         'parameters': [1, 0],
         'x_values': [0, np.pi, 42, -1],

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"""