def gaussian_deconvolve (smaj, smin, spa, bmaj, bmin, bpa):
    """Deconvolve two Gaussians analytically.

    Given the shapes of 2-dimensional “source” and “beam” Gaussians, this
    returns a deconvolved “result” Gaussian such that the convolution of
    “beam” and “result” is “source”.

    Arguments:

    smaj
      Major axis of source Gaussian.
    smin
      Minor axis of source Gaussian.
    spa
      Orientation angle of source Gaussian, in radians.
    bmaj
      Major axis of beam Gaussian.
    bmin
      Minor axis of beam Gaussian.
    bpa
      Orientation angle of beam Gaussian, in radians.

    The return value is ``(rmaj, rmin, rpa, status)``. The first three values
    have the same format as the input arguments. The *status* result is one of
    "ok", "pointlike", or "fail". A "pointlike" status indicates that the
    source and beam shapes are difficult to distinguish; a "fail" status
    indicates that the two shapes seem to be mutually incompatible (e.g.,
    source and beam are very narrow and orthogonal).

    The axes can be measured in any units, so long as they're consistent.

    Ideally if::

      rmaj, rmin, rpa, status = gaussian_deconvolve (smaj, smin, spa, bmaj, bmin, bpa)

    then::

      smaj, smin, spa = gaussian_convolve (rmaj, rmin, rpa, bmaj, bmin, bpa)

    Implementation derived from MIRIAD’s ``gaudfac``. This function currently
    doesn't do a great job of dealing with pointlike sources, i.e. ones where
    “source” and “beam” are nearly indistinguishable.

    """
    # I've added extra code to ensure ``smaj >= bmaj``, ``smin >= bmin``, and
    # increased the coefficient in front of "limit" from 0.1 to 0.5. Feel a
    # little wary about that first change.

    from numpy import cos, sin, sqrt, min, abs, arctan2

    if smaj < bmaj:
        smaj = bmaj
    if smin < bmin:
        smin = bmin

    alpha = ((smaj * cos (spa))**2 + (smin * sin (spa))**2 -
             (bmaj * cos (bpa))**2 - (bmin * sin (bpa))**2)
    beta = ((smaj * sin (spa))**2 + (smin * cos (spa))**2 -
            (bmaj * sin (bpa))**2 - (bmin * cos (bpa))**2)
    gamma = 2 * ((smin**2 - smaj**2) * sin (spa) * cos (spa) -
                 (bmin**2 - bmaj**2) * sin (bpa) * cos (bpa))

    s = alpha + beta
    t = sqrt ((alpha - beta)**2 + gamma**2)
    limit = 0.5 * min ([smaj, smin, bmaj, bmin])**2
    status = 'ok'

    if alpha < 0 or beta < 0 or s < t:
        dmaj = dmin = dpa = 0

        if 0.5 * (s - t) < limit and alpha > -limit and beta > -limit:
            status = 'pointlike'
        else:
            status = 'fail'
    else:
        dmaj = sqrt (0.5 * (s + t))
        dmin = sqrt (0.5 * (s - t))

        if abs (gamma) + abs (alpha - beta) == 0:
            dpa = 0
        else:
            dpa = 0.5 * arctan2 (-gamma, alpha - beta)

    return dmaj, dmin, dpa, status
Deconvolve two Gaussians analytically.

    Given the shapes of 2-dimensional “source” and “beam” Gaussians, this
    returns a deconvolved “result” Gaussian such that the convolution of
    “beam” and “result” is “source”.

    Arguments:

    smaj
      Major axis of source Gaussian.
    smin
      Minor axis of source Gaussian.
    spa
      Orientation angle of source Gaussian, in radians.
    bmaj
      Major axis of beam Gaussian.
    bmin
      Minor axis of beam Gaussian.
    bpa
      Orientation angle of beam Gaussian, in radians.

    The return value is ``(rmaj, rmin, rpa, status)``. The first three values
    have the same format as the input arguments. The *status* result is one of
    "ok", "pointlike", or "fail". A "pointlike" status indicates that the
    source and beam shapes are difficult to distinguish; a "fail" status
    indicates that the two shapes seem to be mutually incompatible (e.g.,
    source and beam are very narrow and orthogonal).

    The axes can be measured in any units, so long as they're consistent.

    Ideally if::

      rmaj, rmin, rpa, status = gaussian_deconvolve (smaj, smin, spa, bmaj, bmin, bpa)

    then::

      smaj, smin, spa = gaussian_convolve (rmaj, rmin, rpa, bmaj, bmin, bpa)

    Implementation derived from MIRIAD’s ``gaudfac``. This function currently
    doesn't do a great job of dealing with pointlike sources, i.e. ones where
    “source” and “beam” are nearly indistinguishable.