Sample a contour in a 2D analytic function. Arguments:

    f
      A function, mapping (x, y) -> z.
    df
      The partial derivative: df (x, y) -> [dz/dx, dz/dy]. If None,
      the derivative of f is approximated numerically with
      scipy.derivative.
    x0
      Initial x value. Should be of "typical" size for the problem;
      avoid 0.
    y0
      Initial y value. Should be of "typical" size for the problem;
      avoid 0.

    Optional arguments:

    maxiters
      Maximum number of points to create. Default 5000.
    defeta
      Initially offset by distances of defeta*[df/dx, df/dy]
      Default 0.05.
    netastep
      Number of steps between defeta and the machine resolution
      in which we test eta values for goodness. (OMG FIXME doc).
      Default 12.
    vtol1
      Tolerance for constancy in the value of the function in the
      initial offset step. The value is only allowed to vary by
      ``f(x0,y0) * vtol1``. Default 1e-3.
    vtol2
      Tolerance for constancy in the value of the function in the
      along the contour. The value is only allowed to vary by
      ``f(x0,y0) * vtol2``. Default 1e-8.
    maxnewt
      Maximum number of Newton's method steps to take when
      attempting to hone in on the desired function value. Default 20.
    dorder
      Number of function evaluations to perform when evaluating
      the derivative of f numerically. Must be an odd integer greater
      than 1. Default 7.
    goright
      If True, trace the contour rightward (as looking uphill),
      rather than leftward (the default).

def analytic_2d (f, df, x0, y0,
                 maxiters=5000,
                 defeta=0.05,
                 netastep=12,
                 vtol1=1e-3,
                 vtol2=1e-8,
                 maxnewt=20,
                 dorder=7,
                 goright=False):
    """Sample a contour in a 2D analytic function. Arguments:

    f
      A function, mapping (x, y) -> z.
    df
      The partial derivative: df (x, y) -> [dz/dx, dz/dy]. If None,
      the derivative of f is approximated numerically with
      scipy.derivative.
    x0
      Initial x value. Should be of "typical" size for the problem;
      avoid 0.
    y0
      Initial y value. Should be of "typical" size for the problem;
      avoid 0.

    Optional arguments:

    maxiters
      Maximum number of points to create. Default 5000.
    defeta
      Initially offset by distances of defeta*[df/dx, df/dy]
      Default 0.05.
    netastep
      Number of steps between defeta and the machine resolution
      in which we test eta values for goodness. (OMG FIXME doc).
      Default 12.
    vtol1
      Tolerance for constancy in the value of the function in the
      initial offset step. The value is only allowed to vary by
      ``f(x0,y0) * vtol1``. Default 1e-3.
    vtol2
      Tolerance for constancy in the value of the function in the
      along the contour. The value is only allowed to vary by
      ``f(x0,y0) * vtol2``. Default 1e-8.
    maxnewt
      Maximum number of Newton's method steps to take when
      attempting to hone in on the desired function value. Default 20.
    dorder
      Number of function evaluations to perform when evaluating
      the derivative of f numerically. Must be an odd integer greater
      than 1. Default 7.
    goright
      If True, trace the contour rightward (as looking uphill),
      rather than leftward (the default).

    """
    # Coerce argument types.

    if not callable (f):
        raise ValueError ('f')

    if df is not None and not callable (df):
        raise ValueError ('df')

    x0 = float (x0)
    if x0 == 0.:
        raise ValueError ('x0')

    y0 = float (y0)
    if y0 == 0.:
        raise ValueError ('y0')

    maxiters = int (maxiters)
    if maxiters < 3:
        raise ValueError ('maxiters')

    defeta = float (defeta)
    if defeta <= 0:
        raise ValueError ('defeta')

    netastep = int (netastep)
    if netastep < 2:
        raise ValueError ('netastep')

    vtol1 = float (vtol1)
    if vtol1 <= 0:
        raise ValueError ('vtol1')

    vtol2 = float (vtol2)
    if vtol2 >= vtol1:
        raise ValueError ('vtol2')

    maxnewt = int (maxnewt)
    if maxnewt < 1:
        raise ValueError ('maxnewt')

    # What value are we contouring?
    v = f (x0, y0)

    # If no derivative is given, use a numerical approximation.

    if df is None:
        derivx = abs (x0 * 0.025)
        derivy = abs (y0 * 0.025)
        from scipy import derivative

        if dorder == 2:
            # simple derivative
            def df (x1, y1):
                z0 = f (x1, y1)
                dx = max (abs (x1) * 1e-5, 1e-8)
                dy = max (abs (y1) * 1e-5, 1e-8)
                dzdx = (f (x1 + dx, y1) - z0) / dx
                dzdy = (f (x1, y1 + dy) - z0) / dy
                return [dzdx, dzdy]
        else:
            def df (x1, y1):
                dx = derivative (lambda x: f (x, y1), x1, derivx, order=dorder)
                dy = derivative (lambda y: f (x1, y), y1, derivy, order=dorder)
                return [dx, dy]

    # Init eta progression.
    rez = np.finfo (np.double).resolution
    if rez > defeta:
        raise PKError ('defeta below resolution!')
    eta_scale = np.exp ((np.log (rez) - np.log (defeta)) / netastep)

    # Init data storage
    n = 1
    pts = np.empty ((maxiters, 2))
    pts[0] = (x0, y0)
    x = x0
    y = y0

    # Quitflag: 0 if first iteration
    #           1 if inited but not yet ok to quit (definition of this below)
    #           2 if ok to quit
    # initquad: 0 if x > 0, y > 0
    #           1 if x < 0, y > 0
    #           2 if x < 0, y < 0
    #           3 if x > 0, y < 0
    # We invert these senses in the in-loop test to make comparison easy.

    quitflag = 0
    initquad = -1

    # Start finding contours.

    while n < maxiters:
        dfdx, dfdy = df (x, y)

        # If we're booting up, remember the quadrant that df/dx points in.
        # Once we've rotated around to the other direction, it is safe to quit
        # once we return close to the original point, since we must have
        # completed a circle.

        if quitflag == 0:
            if dfdx > 0:
                if dfdy > 0:
                    initquad = 0
                else:
                    initquad = 3
            else:
                if dfdy > 0:
                    initquad = 1
                else:
                    initquad = 2
            quitflag = 1
        elif quitflag == 1:
            if dfdx > 0:
                if dfdy > 0:
                    curquad = 2
                else:
                    curquad = 1
            else:
                if dfdy > 0:
                    curquad = 3
                else:
                    curquad = 0

            if curquad == initquad:
                quitflag = 2

        # We will move perpendicular to [df/dx, df/dy], rotating to the left
        # (arbitrarily) from that direction. We need to figure out how far we
        # can safely move in this direction.

        if goright:
            dx = dfdy * defeta
            dy = -dfdx * defeta
        else:
            dx = -dfdy * defeta
            dy = dfdx * defeta

        i = 0

        while i < netastep:
            nx = x + dx
            ny = y + dy
            nv = f (nx, ny)

            # Is the value of the function sufficently close to what
            # we're aiming for?
            if abs (nv / v - 1) < vtol1:
                break

            # No. Try a smaller dx/dy.
            dx *= eta_scale
            dy *= eta_scale
            i += 1
        else:
            # Failed to find a sufficiently small eta (did not break out of
            # loop)
            raise PKError ('failed to find sufficiently small eta: xy %g,%g; '
                           'dv %g; df %g,%g; dxy %g,%g; defeta %g; eta_scale '
                           '%g' % (x, y, nv - v, dfdx, dfdy, dx, dy, defeta,
                                   eta_scale))

        # Now compute a new [df/dx, df/dy], and move along it, finding our way
        # back to the desired value, 'v'. Newton's method should suffice. This
        # loop usually exits after one iteration.

        i = 0

        while i < maxnewt:
            dfdx, dfdy = df (nx, ny)
            df2 = dfdx**2 + dfdy**2
            dv = nv - v

            nx -= dv * dfdx / df2
            ny -= dv * dfdy / df2
            nv = f (nx, ny)

            if abs (nv/v - 1) < vtol2:
                break

            i += 1
        else:
            # Did not break out of loop.
            raise PKError ('failed to converge with Newton\'s method')

        # Ok, we found our next value.
        pts[n] = (nx, ny)
        x = nx
        y = ny
        n += 1

        # Time to stop? Make sure we've gone at least a half-turn so that we
        # don't just exit on the first iteration.
        if quitflag == 2:
            dist2 = (x/x0 - 1)**2 + (y/y0 - 1)**2
            if dist2 < 3 * (dx**2 + dy**2):
                break
    else:
        raise PKError ('needed too many points to close contour')

    # Woohoo! All done.
    return pts[:n]