Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
A quick demonstration of the magnetic field of a circular current loop
import numpy as np
import matplotlib.pyplot as plt
import scipy.special as sp
a = 1.55
current = 1
mu = 4 * np.pi / 10000000
# equal to sqrt(4ra(z^2 + (a+r)^2)^(-1))
def k_val (R, Z):
denominator = np.square(Z) + np.square(np.add(a, R))
K = np.sqrt(np.divide(4 * R * a, denominator))
return K
# equal to sqrt(z^2 + (a+r)^2)
def gamma_val (R, Z):
gamma = np.sqrt(np.square(Z) + np.square(np.add(a, R)))
return gamma
# equal to sqrt(z^2 + (r-a)^2)
def epsilon_val (R, Z):
epsilon = np.square(Z) + np.square(np.subtract(a, R))
return epsilon
R, Z = np.meshgrid(np.arange(0.1, 3, 0.1), np.arange(-3, 3, 0.1))
K = k_val(R, Z)
gamma = gamma_val (R, Z)
epsilon = epsilon_val (R, Z)
E1 = sp.ellipk(np.square(K))
E2 = sp.ellipe(np.square(K))
Ua = np.divide((mu * current) / (2 * np.pi), gamma)
Va = np.multiply(np.divide(np.divide((mu * current) / (2 * np.pi), gamma), R), Z)
U = np.multiply(
Ua,
np.add(
np.multiply(
np.divide(
(np.square(a) - np.square(Z) - np.square(R)),
epsilon),
E2),
E1))
V = np.multiply(
Va,
np.subtract(
np.multiply(
np.divide(
(np.square(a) + np.square(Z) + np.square(R)),
epsilon),
E2),
E1))
plt.figure(1)
ax =plt.gca()
mag = np.sqrt(V**2 + U**2)
Us = np.divide(U, mag)
Vs = np.divide(V, mag)
ax.quiver(R, Z, Vs, Us, angles='xy', scale_units='xy', scale=10 )
ax.set_xlim([0,3])
ax.set_ylim([-3,3])
plt.draw()
plt.figure(2)
ax = plt.gca()
ax.plot (Z[:,10], U[:,10])
ax.plot (Z[:,20], U[:,20])
ax.set_xlim([-3,3])
ax.set_ylim([-0.0000003,0.0000003])
plt.draw()
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment