Skip to content

Instantly share code, notes, and snippets.

@GrantTrebbin
Created November 28, 2015 10:31
Show Gist options
  • Save GrantTrebbin/e9ca33585c36a6174e32 to your computer and use it in GitHub Desktop.
Save GrantTrebbin/e9ca33585c36a6174e32 to your computer and use it in GitHub Desktop.
In plane magnetic field of a current loop of radius a, distance r from the loop axis. Field will be perpendicular to the plane.
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 z^2 + (r-a)^2
def epsilon_val (R, Z):
epsilon = np.square(Z) + np.square(np.subtract(a, R))
return epsilon
# equal to a^2 - z^2 - r^2
def delta_val (R, Z):
delta = np.square(a) - np.square(Z) - np.square(R)
return delta
R = np.arange(0.0, 5, 0.01)
K = k_val(R, 0)
gamma = gamma_val (R, 0)
epsilon = epsilon_val (R, 0)
delta = delta_val (R, 0)
E1 = sp.ellipk(np.square(K))
E2 = sp.ellipe(np.square(K))
Bza = np.divide((mu * current) / (2 * np.pi), gamma)
Bzb = np.multiply(np.divide(delta, epsilon), E2) + E1
Bz = Bza * Bzb
plt.figure(1)
ax = plt.gca()
ax.plot (R, Bz)
ax.set_ylim([-0.000003,0.000003])
plt.draw()
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment