solenoidBNrz.py

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

# --- 1. Physical Constants & Parameters ---
mu0 = 4 * np.pi * 1e-7  # Permeability of free space
mu_r = 2000             # Relative permeability
n = 1000                # Turns per meter
I = 1                   # Current in A
R = 0.01                # Radius (1 cm)
L_half = 0.05           # Half-length (5 cm), Total length = 10 cm

# --- 2. Demagnetization Factor N (Joseph/Chen exact formula) ---
# We use the analytic formula to get the scalar N for the bulk magnetization.
# gamma is the aspect ratio (Length / Diameter)
gamma = (2 * L_half) / (2 * R) 

# Exact N_z for a cylinder (Chen et al. 1991 / Joseph 1966)
# This formulation is numerically stable and dimensionally correct.
def calculate_N_cylinder(gamma):
    return 1 - (4 / (3 * np.pi * gamma)) * (
        np.sqrt(1 + gamma**2) * ( (ellipk(1/(1+gamma**2)) - ellipe(1/(1+gamma**2)))/gamma**2 + ellipe(1/(1+gamma**2)) ) - 1
    )
    # Note: ellipk(m) in scipy takes m = k^2. Here m = 1/(1+gamma^2).

N = calculate_N_cylinder(gamma)

# Calculate Effective Magnetization M
# H_in = H_applied - N * M
# M = chi * H_in  =>  M = chi * (H_app - N * M)
# M (1 + chi * N) = chi * H_app
# M = (chi * H_app) / (1 + chi * N)
chi = mu_r - 1
H_applied = n * I
M_eff = (chi * H_applied) / (1 + chi * N)

print(f"Aspect Ratio (L_total/D): {gamma}")
print(f"Demagnetization Factor N: {N:.5f}")
print(f"Effective Magnetization M: {M_eff:.2f} A/m")

# --- 3. Exact Field Calculation (Off-Axis) ---

def Bz_integrand(z_prime, r, z, R):
    """
    The integrand representing the contribution of a current loop at z_prime 
    to the field at (r, z).
    Derived from the off-axis field of a current loop using Elliptic Integrals.
    """
    dZ = z - z_prime
    sqrt_term = np.sqrt((R + r)**2 + dZ**2)
    m = (4 * R * r) / ((R + r)**2 + dZ**2)
    
    # Elliptic integrals (scipy uses parameter m = k^2)
    K = ellipk(m)
    E = ellipe(m)
    
    # Singularity handling: if r=R and z=z_prime (on the wire), E diverges.
    # The integration routine usually handles integrable singularities, 
    # but we avoid r=R exactly in the plotting lists.
    
    term1 = K
    term2 = ((R**2 - r**2 - dZ**2) / ((R - r)**2 + dZ**2)) * E
    
    return (1.0 / sqrt_term) * (term1 + term2)

def calculate_Bz_total(r_obs, z_points, M_total, R, L_half):
    """
    Integrates the loop contributions over the length of the cylinder.
    Total 'current' comes from two sources:
    1. Real current n*I (Solenoid)
    2. Magnetization current M_eff (Core)
    We treat them as a single equivalent surface current sheet K_tot = n*I + M_eff
    """
    K_total = n * I + M_eff
    factor = (mu0 * K_total) / (2 * np.pi)
    
    B_values = []
    for z in z_points:
        # Integrate from -L to +L
        result, error = quad(Bz_integrand, -L_half, L_half, args=(r_obs, z, R))
        B_values.append(factor * result)
    return np.array(B_values)

# --- 4. Plotting ---
z_plot = np.linspace(0, 0.15, 200) # Plot from center (0) to past the end (15cm)
radii = [0, 0.005, 0.009, 0.0099]  # r = 0, 5mm, 9mm, 9.9mm (near surface)

plt.figure(figsize=(10, 7))

for r in radii:
    # Calculate B array
    B_profile = calculate_Bz_total(r, z_plot, M_eff, R, L_half)
    
    # Labeling
    r_label = f"r = {r*1000:.1f} mm"
    if r == 0: r_label += " (Axis)"
    elif r == 0.0099: r_label += " (Near Surface)"
        
    plt.plot(z_plot * 100, B_profile, label=r_label, linewidth=2)

# Add geometric markers
plt.axvline(x=L_half * 100, color='k', linestyle='--', label='Cylinder End')
plt.text(L_half * 100 + 0.5, plt.ylim()[1]*0.9, 'End Face', rotation=90)

plt.title(f'Magnetic Field $B_z$ Distribution (Total Length = {2*L_half*100} cm)')
plt.xlabel('Axial Position z (cm)')
plt.ylabel('Magnetic Field $B_z$ (Tesla)')
plt.grid(True, which='both', alpha=0.3)
plt.legend(loc='lower left')
plt.tight_layout()
plt.show()