N52.ipynb

n52.ipynb

second.py

# Computing and plotting the magnetic field just above the surface of an N52 cylindrical magnet
# using a surface magnetic-charge model (end-face magnetic charges).
# Geometry: diameter 20 mm, thickness 10 mm, Br = 1.45 T (N52)
# Observation plane: z = +L/2 + epsilon (just above the top face), epsilon = 0.1 mm
# We compute B using B = mu0/(4*pi) * integral(sigma * (r_obs - r')/|r_obs - r'|^3 dA')
# where sigma = +M on top face (z = +L/2) and sigma = -M on bottom face (z = -L/2),
# and M = Br/mu0 (approximate uniform magnetization).
#
# This is a direct numerical surface integral (discretized in polar coordinates).
# We plot B_z (axial component) and |B| vs radius (0 -> 2*R).

import numpy as np
import matplotlib.pyplot as plt
from math import pi

# Parameters (meters)
diameter_mm = 20.0
thickness_mm = 10.0
R = (diameter_mm/2) * 1e-3
L = thickness_mm * 1e-3
Br = 1.45  # Tesla, remanence for N52
mu0 = 4 * np.pi * 1e-7
M = Br / mu0  # approximate magnetization (A/m) 

# Integration settings
Nr = 120       # radial subdivisions on each face
Nphi = 180     # angular subdivisions
r_src = np.linspace(0, R, Nr)
phi_src = np.linspace(0, 2*np.pi, Nphi, endpoint=False)
dr = r_src[1] - r_src[0]
dphi = phi_src[1] - phi_src[0]

# Precompute grid of source points (r', phi')
r_grid, phi_grid = np.meshgrid(r_src, phi_src, indexing='xy')  # shape (Nphi, Nr)
x_src = (r_grid * np.cos(phi_grid)).ravel()
y_src = (r_grid * np.sin(phi_grid)).ravel()
area_src = (r_grid.ravel() * dr * dphi)  # area element r' dr dphi

# Face z positions and sigma
z_top = +L/2
z_bot = -L/2
sigma_top = M
sigma_bot = -M

# Source coordinates for top and bottom faces (stacked)
x_top = x_src.copy()
y_top = y_src.copy()
z_top_arr = np.full_like(x_top, z_top)
area_top = area_src.copy()

x_bot = x_src.copy()
y_bot = y_src.copy()
z_bot_arr = np.full_like(x_bot, z_bot)
area_bot = area_src.copy()

# Observation radii (we'll compute from 0 to 2R)
r_obs = np.linspace(0.0, 2*R, 200)
z_obs = z_top + 1e-4  # 0.1 mm above the top face (in meters)

# Prepare result arrays
Bz_vals = np.zeros_like(r_obs)
Bmag_vals = np.zeros_like(r_obs)

const = mu0 / (4*np.pi)

# Loop over observation radii (vectorized over sources)
for i, ro in enumerate(r_obs):
    # place observation at (ro, 0, z_obs) (use symmetry axis phi=0)
    x_o = ro
    y_o = 0.0
    z_o = z_obs
    
    # Top face contributions
    dx = x_o - x_top
    dy = y_o - y_top
    dz = z_o - z_top_arr
    rvec_sq = dx*dx + dy*dy + dz*dz
    rmag = np.sqrt(rvec_sq)
    # avoid singularity by ensuring minimum distance
    rmag[rmag < 1e-12] = 1e-12
    inv_r3 = 1.0 / (rmag**3)
    # dB = const * sigma * dA * rvec * inv_r3
    dBx_top = const * sigma_top * area_top * dx * inv_r3
    dBy_top = const * sigma_top * area_top * dy * inv_r3
    dBz_top = const * sigma_top * area_top * dz * inv_r3
    
    # Bottom face contributions
    dx = x_o - x_bot
    dy = y_o - y_bot
    dz = z_o - z_bot_arr
    rvec_sq = dx*dx + dy*dy + dz*dz
    rmag = np.sqrt(rvec_sq)
    rmag[rmag < 1e-12] = 1e-12
    inv_r3 = 1.0 / (rmag**3)
    dBx_bot = const * sigma_bot * area_bot * dx * inv_r3
    dBy_bot = const * sigma_bot * area_bot * dy * inv_r3
    dBz_bot = const * sigma_bot * area_bot * dz * inv_r3
    
    # Sum contributions
    Bx = dBx_top.sum() + dBx_bot.sum()
    By = dBy_top.sum() + dBy_bot.sum()
    Bz = dBz_top.sum() + dBz_bot.sum()
    
    Bz_vals[i] = Bz
    Bmag_vals[i] = np.sqrt(Bx*Bx + By*By + Bz*Bz)

# Convert radii back to mm and fields to mT for plotting convenience
r_obs_mm = r_obs * 1e3
Bz_mT = Bz_vals * 1e3
Bmag_mT = Bmag_vals * 1e3

# Plot
plt.figure(figsize=(7,4.5))
plt.plot(r_obs_mm, Bz_mT, label='B_z (mT)')
plt.plot(r_obs_mm, Bmag_mT, label='|B| (mT)')
plt.axvline(x=diameter_mm/2, linestyle='--', linewidth=0.8, label='magnet edge')
plt.xlabel('Radial distance from center (mm)')
plt.ylabel('Magnetic field (mT)')
plt.title('N52 disk (D=20 mm, L=10 mm) — field 0.1 mm above top face')
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

# Print a short summary numeric values at a few radii
for rr, bz, bmag in zip([0.0, R*1e3/2, R*1e3, 1.5*R*1e3, 2*R*1e3],
                        np.interp([0.0, R*1e3/2, R*1e3, 1.5*R*1e3, 2*R*1e3], r_obs_mm, Bz_mT),
                        np.interp([0.0, R*1e3/2, R*1e3, 1.5*R*1e3, 2*R*1e3], r_obs_mm, Bmag_mT)):
    print(f"r = {rr:.1f} mm: B_z ≈ {bz:.1f} mT, |B| ≈ {bmag:.1f} mT")

https://youtu.be/dCeycaEkmh4?si=bRFiDqhXg2o1iICR

https://youtu.be/vrPwx9sHSbM?si=UpY_oucalEkW5mrs

https://youtu.be/DgwwWY4Ji2Q?si=d6eVvDLtYEGaCX17

https://www.kjmagnetics.com/magnetic-field-calculator.asp