nCoils.py

import numpy as np
from scipy.special import ellipk, ellipe
import matplotlib.pyplot as plt

# Constants (normalized: mu0=1, I=1, a=1 for simplicity)
mu0 = 1.0
I = 1.0
a = 1.0
d = 0.01  # Small d to see changes with large n

# Function to compute B_rho and B_z at (rho, z) for a single coil centered at z_coil
def get_B(rho, z, R=a, I=I, mu0=mu0, z_coil=0):
    z_rel = z - z_coil
    if rho == 0:
        B_rho = 0
        B_z = (mu0 * I * R**2) / (2 * (R**2 + z_rel**2)**(3/2))
        return B_rho, B_z
    sqrt_term = np.sqrt((R + rho)**2 + z_rel**2)
    k2 = 4 * R * rho / ((R + rho)**2 + z_rel**2)
    if k2 >= 1 or k2 < 0:
        return 0, 0
    K = ellipk(k2)
    E = ellipe(k2)
    denom = (R - rho)**2 + z_rel**2
    if denom < 1e-10:
        return np.nan, np.nan
    factor = mu0 * I / (2 * np.pi * sqrt_term)
    term_z = K + E * (R**2 - rho**2 - z_rel**2) / denom
    B_z = factor * term_z
    term_rho = -K + E * (R**2 + rho**2 + z_rel**2) / denom
    B_rho = factor * (z_rel / rho) * term_rho
    return B_rho, B_z

# Create x values (from -3a to 3a)
x_vals = np.linspace(-3*a, 3*a, 500)  # Reduced resolution for speed with large n
z_point = 0.0  # Evaluate in the plane of the first coil (z=0)

# List of n values
n_values = [1, 10, 100, 1000]

# Prepare plot
fig, ax = plt.subplots(figsize=(10, 6))
colors = ['b', 'g', 'r', 'm']  # Colors for each n

for idx, n in enumerate(n_values):
    # Coils placed at z = 0, d, 2d, ..., n*d
    z_coils = [k * d for k in range(n + 1)]
    
    B_mag = np.zeros_like(x_vals)
    
    for i, x in enumerate(x_vals):
        rho = np.abs(x)
        B_rho_tot = 0.0
        B_z_tot = 0.0
        for zc in z_coils:
            br, bz = get_B(rho, z_point, z_coil=zc)
            if not np.isnan(br):
                B_rho_tot += br
            if not np.isnan(bz):
                B_z_tot += bz
        B_mag[i] = np.sqrt(B_rho_tot**2 + B_z_tot**2)
    
    ax.plot(x_vals, B_mag, color=colors[idx], linewidth=1.5, label=f'n={n} (total {n+1} coils)')

ax.set_xlabel(r'$x$')
ax.set_ylabel(r'Magnetic Field Magnitude $|B|$')
ax.set_title('Magnetic Field Magnitude vs x in the Plane of the First Coil for Various n (d=0.01)')
ax.grid(True)
ax.legend()
plt.tight_layout()
plt.show()