shredEngineer
commented
1 year ago 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 ,
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