Example is requested for python

[Steffen-W]

Hi @shredEngineer ,

your project looks good,
can you give me an example of how I can calculate an air coil without the gui but only via python?

I would really appreciate it

[shredEngineer]

Hey Steffen,

thanks for the kind words.

MagnetiCalc is very monolithic in that regard; putting all needed classes (SamplingVolume, Wire, BiotSavart, ...) together in a minimal example would take some work which I have not done yet.

All I can do right now is to give you MagnetiCalc v1.0 which is at least all in one file, but it also has some extra stuff in there, and it does NOT use any optimization, so it all runs on one core, without any speedups:

# (C) <mailto:anfrage@paulwilhelm.de>
import numpy as np
import matplotlib.pyplot as plt

# Parametrisierung und Anzeige
animation = False       # Automatische 3D-Rotation
azim, elev = 270, 90    # Anfängliche 3D-Rotation
loops = 2               # Anzahl der Schleifen
radius = 2              # Radius der Gesamtschleife
height = 3              # Höhe einer Schleife
points = False          # Zeige nummerierte Schleifenpunkte
components = True       # Zeige die Komponenten der Parametrisierung der einzelnen Schleife
slicing = 16            # Elemente der Gesamtschleife in noch kleinere Segmente zerlegen
grid = True             # Koordinatensystem anzeigen

# Berechnung der magnetischen Flussdichte
I_DC = .25              # Stromstärke
B_box = 5               # "Radius" (halbe Seitenlänge) des zu berechnenden (kubischen) Volumens (zentriert um Nullpunkt)
B_res = 151             # Auflösung des Gitters
B_normalize = True      # Vektoren normalisieren
B_scale = .15           # Skalierung der Vektoren
B_alpha = 1.0           # Sichtbarkeit der Vektoren
B_limit = 100           # Vektoren mit größerem Betrag werden darauf begrenzt

torus = True            # Begrenze Feldberechnung auf einen Torus:
torus_shell_min = 1.5   # Radius der inneren Kugelschale
torus_shell_max = 2.5   # Radius der äußeren Kugelschale
torus_z_min = 0         # Untere Z-Ebene
torus_z_max = 0         # Obere Z-Ebene

# Parametrisierung der einzelnen Schleife
#              0  1  2  3  4  5  6  7  8  9 10 11
x = -np.array([0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3]) / 3 + 1/2
y = -np.array([0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0]) + 1/2
z = +np.array([0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0]) - 1/2
p = np.linspace(0, 1, len(x))

# Erzeuge die gewünschte Anzahl Schleifen mit vorgegebenem Radius (Gesamtschleife), schließe die Gesamtschleife
X, Y, Z = [], [], []
for a in np.linspace(0, 2 * np.pi, loops, endpoint=False):
    X = np.append(X, x * np.cos(a) - (y + radius) * np.sin(a))
    Y = np.append(Y, x * np.sin(a) + (y + radius) * np.cos(a))
    Z = np.append(Z, z * height)
X, Y, Z = np.append(X, X[0]), np.append(Y, Y[0]), np.append(Z, Z[0])
P = np.linspace(0, 1, len(X))

# Linienelemente der Gesamtschleife in noch kleinere Segmente zerlegen, schließe die Gesamtschleife
nX, nY, nZ = [], [], []
for i in range(len(P) - 1):
    p1 = np.array([X[i], Y[i], Z[i]])               # Startpunkt dieses Elements
    p2 = np.array([X[i + 1], Y[i + 1], Z[i + 1]])   # Endpunkt dieses Elements
    dI = p2 - p1                                    # Linienelement
    for j in range(slicing):
        nX = np.append(nX, X[i] + dI[0] * j / slicing)
        nY = np.append(nY, Y[i] + dI[1] * j / slicing)
        nZ = np.append(nZ, Z[i] + dI[2] * j / slicing)
X, Y, Z = np.append(nX, nX[0]), np.append(nY, nY[0]), np.append(nZ, nZ[0])
P = np.linspace(0, 1, len(X))

# Zeige Matplotlib-Fenster
fig = plt.figure(
    "MagnetiCalc-v1.0: Magnetic flux density calculation example <mailto:anfrage@paulwilhelm.de>  ::  Drücke [q] zum Beenden",
    figsize=(12 if components else 6, 6)
)

if components:
    # Zeige die Komponenten der Parametrisierung der einzelnen Schleife
    ax_x = fig.add_subplot(322); ax_x.set(xlabel="$p$", ylabel="$x(p)$"); ax_x.plot(p, x); ax_x.scatter(p, x)
    ax_y = fig.add_subplot(324); ax_y.set(xlabel="$p$", ylabel="$y(p)$"); ax_y.plot(p, y); ax_y.scatter(p, y)
    ax_z = fig.add_subplot(326); ax_z.set(xlabel="$p$", ylabel="$z(p)$"); ax_z.plot(p, z); ax_z.scatter(p, z)

# Zeige die Gesamtschleife
ax_3d = fig.add_subplot(121 if components else 111, projection='3d', proj_type='ortho')
ax_3d.set(xlabel="$x$", ylabel="$y$", zlabel="$z$")
ax_3d.set_title(f"Radius = {radius}, Loops = {loops}", y=-0.086)
ax_3d.plot(X, Y, Z, c="r")

if not grid:
    # Verstecke das Koordinatensystem
    ax_3d.set_axis_off()

if points:
    # Zeige nummerierte Schleifenpunkte
    ax_3d.scatter(X, Y, Z)
    for i in range(len(P)):
        ax_3d.text(X[i], Y[i], Z[i], i)

# Berechne magnetische Flussdichte
side = np.linspace(-B_box, B_box, B_res)    # Koordinaten einer Gitterkante

# Erzeuge flaches Array der 3D-Gitterkoordinaten (Vektoren) und der Feldvektoren an jeder Gitterkoordinate
G, B = [], []
for gx in range(B_res):
    print(f"{100 * (gx + 1) / B_res:.2f} %")
    for gy in range(B_res):
        for gz in range(B_res):
            g = np.array([side[gx], side[gy], side[gz]])

            if torus:
                # Überspringe alle Punkte, die nicht auf Kreisscheiben in der XY-Ebene liegen
                g_dist = np.linalg.norm(np.array([side[gx], side[gy]]))
                if g_dist < torus_shell_min or g_dist > torus_shell_max:
                    continue

                # Überspringe alle Punkte, die nicht zwischen zwei Z-Ebenen liegen
                if g[2] < torus_z_min or g[2] > torus_z_max:
                    continue

            B_tot = np.zeros(3)  # Vektorieller Akkumulator

            # Iteriere über die Punkte der Gesamtschleife (allerletzter Punkt = allererster Punkt)
            for i in range(len(P) - 1):

                # Berechne Linienelement dI
                p1 = np.array([X[i], Y[i], Z[i]])               # Startpunkt dieses Elements
                p2 = np.array([X[i + 1], Y[i + 1], Z[i + 1]])   # Endpunkt dieses Elements
                dI = p2 - p1                                    # Linienelement
                m = p1 + dI / 2                                 # Mittelpunkt des Linienelements

                # Vektor und Abstand vom Mittelpunkt des Linienelements zur Gitterkoordinate
                R = g - m
                r = np.linalg.norm(R)

                # Berechne Flussdichteelement dB (Biot-Savart)
                dB = I_DC * np.cross(dI, R) / (r ** 3)

                # Vektorielle Superposition
                B_tot += dB

            # Betrag des gesamten Vektors
            B_mag = np.linalg.norm(B_tot)

            # Vektorlänge hart begrenzen
            if B_mag > B_limit:
                B_tot = B_tot / B_mag * B_limit

            # Speichere Koordinate und Feldvektor
            G.append(g)
            B.append(B_tot)

# Transponiere Gitter und Feld
Gx, Gy, Gz = [g[0] for g in G], [g[1] for g in G], [g[2] for g in G]
Bx, By, Bz = [b[0] for b in B], [b[1] for b in B], [b[2] for b in B]

# Zeige den Vektor der magnetischen Flussdichte an jeder Koordinate des 3D-Gitters
ax_3d.quiver(Gx, Gy, Gz, Bx, By, Bz, pivot="middle", length=B_scale, normalize=B_normalize, alpha=B_alpha)

# Erfrische Matplotlib-Fenster
plt.tight_layout()
figManager = plt.get_current_fig_manager()
figManager.window.showMaximized()

# Automatische 3D-Rotation
if animation:
    while plt.fignum_exists(fig.number):
        ax_3d.view_init(elev, azim)
        # azim += 1
        elev -= 1
        plt.draw()
        plt.pause(.001)
else:
    ax_3d.view_init(elev, azim)
    plt.show()

Let me know if this helps! :)

[Steffen-W]

Hi @shredEngineer,

Thanks a lot the code looks really good. I will have a closer look and try it out the day after tomorrow. I would like to determine the inductance but in principle I can calculate the integral with the help of the field strength in space. That is facilitated with the code in any case. Thanks a lot