<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 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