def B(t, r):
    """Returns magnetic field of a current loop in point r = (r1, r2, r3) and time t
    The expressions for the magnatic field of the current loop are taken from 
    [J. Simpson et al. Simple Analytic Expressions for the Magnetic Field of a Circular Current Loop. 
    Transactions on Magnetics (2001)]
    
    Arguments:
        r -- vector of position
        t -- time
        
    Returns:
        B -- vector of magnetic field components
    """
    x,y,z = r
    
    a = 0.03 # diameter of the loop in meters
    mu0 = 4*np.pi*1e-7 # vacuum permeability
    I = 10 # current in Amps
    #D = 0.1 # distance of two current loops
    C = mu0 * I / np.pi
    
    rho2 = np.power(x,2) + np.power(y,2)
    rho = np.sqrt(rho2)
    r2 = np.power(x,2) + np.power(y,2) + np.power(z,2)
    a2 = np.power(a, 2)
    alpha2 = a2 + r2 - 2*a*rho
    beta2 = a2 + r2 + 2*a*rho
    beta = np.sqrt(beta2)
    k2 = 1 - alpha2/beta2
    
    #print(ellipk(k2))
    #print(2 * alpha2 * beta * rho2)
    
    Bx = (C * x * z) / (2 * alpha2 * beta * rho2) * ( (a2 + r2) * ellipe(k2) - alpha2 * ellipk(k2) )
    By = Bx * y / x
    Bz = C / ( 2 * alpha2 * beta ) * ( (a2 - r2) * ellipe(k2) + alpha2 * ellipk(k2) )
    
    # Function diverges if x = y = 0, in which case we fix it.
    if x == 0 or y == 0:
        Bx, By = (0, 0)
    
    return (Bx, By, Bz)