def magnetic_field(r, n, r0, R):
    """
    Returns the magnetic field from an arbitrary current loop calculated from
    eqns (1) and (2) in Phys Rev A Vol. 35, N 4, pp. 1535-1546; 1987.

    Arguments
    ----------
        n: ndarray, shape (3, )
            The normal vector to the plane of the loop at the center,
            current is oriented by the right-hand-rule.
        r: ndarray, shape (m, 3)
            A position vector where the magnetic field is evaluated:
            [x1 y2 z3 ; x2 y2 z2 ; ... ]
            r is in units of d
        r0: ndarray, shape (3, )
            The location of the center of the loop in units of d: [x y z]
        R: float
            The radius of the current loop

    Returns
    --------
    B: ndarray, shape (m, 3)
        a vector for the B field at each position specified in r
        in inverse units of (mu I) / (2 pi d)
        for I in amps and d in meters and mu = 4 pi * 10^-7 we get Tesla
    """
    ### Translate the coordinates in the coil's frame
    n, l, m = base_vectors(n)

    # transformation matrix coil frame to lab frame
    trans = np.vstack((l, m, n))
    # transformation matrix to lab frame to coil frame
    inv_trans = np.linalg.inv(trans)

    # point location from center of coil
    r = r - r0
    # transform vector to coil frame
    r = np.dot(r, inv_trans)

    #### calculate field

    # express the coordinates in polar form
    x = r[:, 0]
    y = r[:, 1]
    z = r[:, 2]
    rho = np.sqrt(x**2 + y**2)
    theta = np.arctan(x/y)
    theta[y==0] = 0

    E = special.ellipe((4 * R * rho)/( (R + rho)**2 + z**2))
    K = special.ellipk((4 * R * rho)/( (R + rho)**2 + z**2))
    Bz =  1/np.sqrt((R + rho)**2 + z**2) * (
                K
              + E * (R**2 - rho**2 - z**2)/((R - rho)**2 + z**2)
              )
    Brho = z/(rho*np.sqrt((R + rho)**2 + z**2)) * (
               -K
              + E * (R**2 + rho**2 + z**2)/((R - rho)**2 + z**2)
              )
    # On the axis of the coil we get a divided by zero here. This returns a
    # NaN, where the field is actually zero :
    Brho[np.isnan(Brho)] = 0
    Brho[np.isinf(Brho)] = 0
    Bz[np.isnan(Bz)]     = 0
    Bz[np.isinf(Bz)]     = 0

    B = np.c_[np.cos(theta)*Brho, np.sin(theta)*Brho, Bz ]

    # Rotate the field back in the lab's frame
    B = np.dot(B, trans)
    return B
dataset/ETHPy150Open enthought/mayavi/docs/source/mayavi/auto/magnetic_field.py/magnetic_field/original