import math


def vortex_feldtkeller(beta, center, right_handed=True, polarity=+1):
    """
    Returns a function f: (x,y,z) -> m representing a vortex
    magnetisation pattern where the vortex lies in the x/y-plane (i.e.
    the magnetisation is constant along the z-direction), the vortex
    core is centered around the point `center` and the m_z component
    falls off exponentially as given by the following formula, which
    is a result by Feldtkeller and Thomas [1].:

        m_z = exp(-2*r^2/beta^2).

    Here `r` is the distance from the vortex core centre and `beta` is
    a parameter, whose value is taken from the first argument to this
    function.

    [1] E. Feldtkeller and H. Thomas, "Struktur und Energie von
        Blochlinien in Duennen Ferromagnetischen Schichten", Phys.
        Kondens. Materie 8, 8 (1965).

    """
    beta_sq = beta ** 2

    def f(pt):
        x = pt[0]
        y = pt[1]
        # To start with, create a right-handed vortex with polarity 1.
        xc = x - center[0]
        yc = y - center[1]
        phi = math.atan2(yc, xc)
        r_sq = xc ** 2 + yc ** 2
        mz = math.exp(-2.0 * r_sq / beta_sq)
        mx = -math.sqrt(1 - mz * mz) * math.sin(phi)
        my = math.sqrt(1 - mz * mz) * math.cos(phi)

        # If we actually want a different polarity, flip the z-coordinates
        if polarity < 0:
            mz = -mz

        # Adapt the chirality accordingly
        if ((polarity > 0) and (not right_handed)) or\
           ((polarity < 0) and right_handed):
            mx = -mx
            my = -my

        return (mx, my, mz)

    return f