solenoid_example.py

import numpy as np
import matplotlib.pyplot as plt
import isosurface as isosurface


def create_magnetic_field_solenoid(grid_size=50, radius=1.5, length=6.0,
                                   n_turns=50, current=10.0):
    """
    Returns:
    --------
    field_data : 3D scalar field (magnitude of B field)
    bounds : spatial bounds
    """
    bounds = {'x': [-5, 5], 'y': [-15, 15], 'z': [-5, 5]}
    x = np.linspace(bounds['x'][0], bounds['x'][1], grid_size)
    y = np.linspace(bounds['y'][0], bounds['y'][1], grid_size)
    z = np.linspace(bounds['z'][0], bounds['z'][1], grid_size)
    X, Y, Z = np.meshgrid(x, y, z, indexing='ij')
    Bx = np.zeros_like(X)
    By = np.zeros_like(Y)
    Bz = np.zeros_like(Z)

    n_segments_per_turn = 40
    mu_0_over_4pi = 1.0

    # i switched y and x, so solenoid is parallel with y axis
    for turn_idx in range(n_turns):
        if turn_idx % 10 == 0:
            print(f"  Progress: {turn_idx}/{n_turns} turns")
        z_turn = -length / 2 + (turn_idx / (n_turns - 1)) * length if n_turns > 1 else 0
        theta = np.linspace(0, 2 * np.pi, n_segments_per_turn, endpoint=False)
        dtheta = 2 * np.pi / n_segments_per_turn
        for th in theta:
            dl_x = -radius * np.sin(th) * dtheta
            dl_y = 0.0
            dl_z = radius * np.cos(th) * dtheta
            r_source_x = radius * np.cos(th)
            r_source_y = z_turn
            r_source_z = radius * np.sin(th)
            rx = X - r_source_x
            ry = Y - r_source_y
            rz = Z - r_source_z
            r_mag = np.sqrt(rx ** 2 + ry ** 2 + rz ** 2)
            r_mag = np.maximum(r_mag, 0.01)
            r_mag_cubed = r_mag ** 3
            cross_x = dl_y * rz - dl_z * ry
            cross_y = dl_z * rx - dl_x * rz
            cross_z = dl_x * ry - dl_y * rx
            factor = mu_0_over_4pi * current / r_mag_cubed
            Bx += factor * cross_x
            By += factor * cross_y
            Bz += factor * cross_z
    B_magnitude = np.sqrt(Bx ** 2 + By ** 2 + Bz ** 2)
    # i found log makes for better visual
    field_data = np.log10(B_magnitude + 1e-10)
    return field_data, bounds


field_data, bounds = create_magnetic_field_solenoid(grid_size=40)
print("creating isosurfaces...")
iso = isosurface.Isosurface(field_data, bounds)

fig, ax = iso.plot_isosurfaces(num_surfaces=11, opacity=0.2, fade_rate=1.0)
ax.set_box_aspect((1, 3, 1))
ax.set_title('solenoid')


# copy solenoid wire (i switched y and x)
def solenoid_wire_xyz():
    wire_x = []
    wire_y = []
    wire_z = []
    radius = 1.45
    length = 6.0
    n_turns = 50
    n_segments_per_turn = 20
    for turn_idx in range(n_turns):
        z_turn = -length / 2 + (turn_idx / (n_turns - 1)) * length if n_turns > 1 else 0
        theta = np.linspace(0, 2 * np.pi, n_segments_per_turn, endpoint=False)
        for th in theta:
            wire_x.append(radius * np.cos(th))
            wire_y.append(z_turn)
            wire_z.append(radius * np.sin(th))
    return wire_x, wire_y, wire_z


ax.plot3D(*solenoid_wire_xyz(), color='black', linewidth=2)
plt.show()