Given an 2D Gaussian expressed as an ellipse (major, minor, pa), compute a
    "squared distance parameter" such that

       z = exp (-0.5 * d2)

    Inputs:

    * x0: position of Gaussian center on x axis
    * y0: position of Gaussian center on y axis
    * mjr: major axis (sigma not FWHM) of the Gaussian
    * mnr: minor axis (sigma not FWHM) of the Gaussian
    * pa: position angle (from +x to +y) of the Gaussian, radians
    * x: x coordinates of the locations for which to evaluate d2
    * y: y coordinates of the locations for which to evaluate d2

    Returns: d2, distance parameter defined as above.

    x0, y0, mjr, and mnr may be in any units so long as they're consistent. x
    and y may be arrays (of the same shape), in which case d2 will be an array
    as well.

def elld2 (x0, y0, mjr, mnr, pa, x, y):
    """Given an 2D Gaussian expressed as an ellipse (major, minor, pa), compute a
    "squared distance parameter" such that

       z = exp (-0.5 * d2)

    Inputs:

    * x0: position of Gaussian center on x axis
    * y0: position of Gaussian center on y axis
    * mjr: major axis (sigma not FWHM) of the Gaussian
    * mnr: minor axis (sigma not FWHM) of the Gaussian
    * pa: position angle (from +x to +y) of the Gaussian, radians
    * x: x coordinates of the locations for which to evaluate d2
    * y: y coordinates of the locations for which to evaluate d2

    Returns: d2, distance parameter defined as above.

    x0, y0, mjr, and mnr may be in any units so long as they're consistent. x
    and y may be arrays (of the same shape), in which case d2 will be an array
    as well.

    """
    _ellcheck (mjr, mnr, pa)

    dx, dy = x - x0, y - y0
    c, s = np.cos (pa), np.sin (pa)
    a = c * dx + s * dy
    b = -s * dx + c * dy
    return (a / mjr)**2 + (b / mnr)**2