moffat.py

class EllipticalMoffat2D(Fittable2DModel):
    """
    Two dimensional elliptical Moffat model.

    Parameters
    ----------
    amplitude : float
        Amplitude of the model.
    x_0 : float
        x position of the maximum of the Moffat model.
    y_0 : float
        y position of the maximum of the Moffat model.
    x_gamma : float
        x core width of the Moffat model.
    y_gamma : float
        y core width of the Moffat model.
    alpha : float
        Power index of the Moffat model.
    theta : float
        Rotation.

    See Also
    --------
    Gaussian2D, Box2D

    Notes
    -----
    Model formula:

    .. math::

        f(x, y) = A \\left(1 + \\frac{\\left(x - x_{0}\\right)^{2} +
        \\left(y - y_{0}\\right)^{2}}{\\gamma^{2}}\\right)^{- \\alpha}
    """

    amplitude = Parameter(default=1)
    x_0 = Parameter(default=0)
    y_0 = Parameter(default=0)
    x_gamma = Parameter(default=1)
    y_gamma = Parameter(default=1)
    alpha = Parameter(default=1)
    theta = Parameter(default=0)

    @property
    def x_fwhm(self):
        """Moffat full width at half maximum."""
        return 2.0 * self.x_gamma * np.sqrt(2.0 ** (1.0 / self.alpha) - 1.0)

    @property
    def y_fwhm(self):
        """Moffat full width at half maximum."""
        return 2.0 * self.y_gamma * np.sqrt(2.0 ** (1.0 / self.alpha) - 1.0)

    @staticmethod
    def evaluate(x, y, amplitude, x_0, y_0, x_gamma, y_gamma, alpha, theta):
        """Two dimensional Moffat model function"""

        cost2 = np.cos(theta) ** 2
        sint2 = np.sin(theta) ** 2
        sin2t = np.sin(2. * theta)
        xstd2 = x_gamma ** 2
        ystd2 = y_gamma ** 2
        xdiff = x - x_0
        ydiff = y - y_0

        a = ((cost2 / xstd2) + (sint2 / ystd2))
        b = ((sin2t / xstd2) - (sin2t / ystd2))
        c = ((sint2 / xstd2) + (cost2 / ystd2))

        rr_gg = (a * xdiff ** 2) + (b * xdiff * ydiff) + (c * ydiff ** 2)
        return amplitude * (1 + rr_gg) ** (-alpha)