(r, n, r0, R)
return B
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.

    Parameters
    ----------
        n is normal vector to the plane of the loop at the center, current
            is oriented by the right-hand-rule.
        r is a position vector where the Bfield is evaluated:
            [x1 y2 z3 ; x2 y2 z2 ; ... ]
        r is in units of d
        r0 is the location of the center of the loop in units of d: [x y z]
        R is the radius of the loop

    Returns
    -------
        B is a vector for the B field at point 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
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.
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 .
def B_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.

    Parameters
    ----------
        n is normal vector to the plane of the loop at the center, current
            is oriented by the right-hand-rule.
        r is a position vector where the Bfield is evaluated:
            [x1 y2 z3 ; x2 y2 z2 ; ... ]
        r is in units of d
        r0 is the location of the center of the loop in units of d: [x y z]
        R is the radius of the loop

    Returns
    -------
        B is a vector for the B field at point 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 = linalg.inv(trans)

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

    #### 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)
    # NaNs are generated where y is zero.
    theta[y == 0] = np.pi / 2

    E = special.ellipe((4 * R * rho)/( (R + rho)**2 + z**2))
    K = special.ellipk((4 * R * rho)/( (R + rho)**2 + z**2))
    dist = ((R - rho)**2 + z**2)
    Bz = 1 / np.sqrt((R + rho)**2 + z**2) * (
                K
              + E * (R**2 - rho**2 - z**2) / dist
              )
    Brho = z / (rho*np.sqrt((R + rho)**2 + z**2)) * (
               -K
              + E * (R**2 + rho**2 + z**2)/ dist
              )
    # On the axis of the coil we get a divided by zero here. This returns a
    # NaN, where the field is actually zero :
    Brho[dist == 0] = 0
    Brho[rho == 0] = 0
    Bz[dist == 0] = 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
def B_field ( r , n , r0 , R ) : ### 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 = linalg . inv ( trans ) r = r - r0 #point location from center of coil r = np . dot ( r , inv_trans ) #transform vector to coil frame #### 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 ) # NaNs are generated where y is zero. theta [ y == 0 ] = np . pi / 2 E = special . ellipe ( ( 4 * R * rho ) / ( ( R + rho ) ** 2 + z ** 2 ) ) K = special . ellipk ( ( 4 * R * rho ) / ( ( R + rho ) ** 2 + z ** 2 ) ) dist = ( ( R - rho ) ** 2 + z ** 2 ) Bz = 1 / np . sqrt ( ( R + rho ) ** 2 + z ** 2 ) * ( K + E * ( R ** 2 - rho ** 2 - z ** 2 ) / dist ) Brho = z / ( rho * np . sqrt ( ( R + rho ) ** 2 + z ** 2 ) ) * ( - K + E * ( R ** 2 + rho ** 2 + z ** 2 ) / dist ) # On the axis of the coil we get a divided by zero here. This returns a # NaN, where the field is actually zero : Brho [ dist == 0 ] = 0 Brho [ rho == 0 ] = 0 Bz [ dist == 0 ] = 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