mirror of
https://github.com/atomic14/kicad-coil-plugins.git
synced 2024-10-18 09:06:57 +00:00
Cleaned up simulation code
This commit is contained in:
Chris Greening
parent
1461c6799b
commit
fcc9c6727c
10 changed files with 411 additions and 1419 deletions
|
|
@ -714,7 +714,7 @@
|
||||||
},
|
},
|
||||||
"vscode": {
|
"vscode": {
|
||||||
"interpreter": {
|
"interpreter": {
|
||||||
"hash": "1ce20143987840b9786ebb5907032c9c3a8efacbb887dbb0ebc4934f2ad26cb3"
|
"hash": "fc384f9db26c31784edfba3761ba3d2c7b2f9b8a63e03a9eb0778fc35334efe1"
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
},
|
},
|
||||||
|
|
|
||||||
|
|
@ -693,7 +693,7 @@
|
||||||
},
|
},
|
||||||
"vscode": {
|
"vscode": {
|
||||||
"interpreter": {
|
"interpreter": {
|
||||||
"hash": "1ce20143987840b9786ebb5907032c9c3a8efacbb887dbb0ebc4934f2ad26cb3"
|
"hash": "fc384f9db26c31784edfba3761ba3d2c7b2f9b8a63e03a9eb0778fc35334efe1"
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
},
|
},
|
||||||
|
|
|
||||||
|
|
@ -4,9 +4,4 @@ numpy
|
||||||
pandas
|
pandas
|
||||||
scikit-spatial
|
scikit-spatial
|
||||||
pre-commit
|
pre-commit
|
||||||
tqdm
|
|
||||||
plotly
|
plotly
|
||||||
magpylib
|
|
||||||
pyvista
|
|
||||||
mayavi
|
|
||||||
PySide
|
|
||||||
File diff suppressed because one or more lines are too long
|
Before Width: | Height: | Size: 36 KiB |
|
|
@ -68,7 +68,7 @@
|
||||||
"# read in the results\n",
|
"# read in the results\n",
|
||||||
"total_volume = bs.read_target_volume(COIL_NAME + \"volume.vol\")\n",
|
"total_volume = bs.read_target_volume(COIL_NAME + \"volume.vol\")\n",
|
||||||
"\n",
|
"\n",
|
||||||
"total_volume = total_volume\n",
|
"total_volume = -total_volume\n",
|
||||||
"\n",
|
"\n",
|
||||||
"# reads the volume we created\n",
|
"# reads the volume we created\n",
|
||||||
"bs.plot_fields(total_volume, (4, 4, 4), (-2, -2, -2), 0.1, which_plane=\"z\", level=0.1)"
|
"bs.plot_fields(total_volume, (4, 4, 4), (-2, -2, -2), 0.1, which_plane=\"z\", level=0.1)"
|
||||||
|
|
@ -204,16 +204,30 @@
|
||||||
"# make the plot larger\n",
|
"# make the plot larger\n",
|
||||||
"plt.rcParams[\"figure.figsize\"] = (10, 10)\n",
|
"plt.rcParams[\"figure.figsize\"] = (10, 10)\n",
|
||||||
"\n",
|
"\n",
|
||||||
|
"# plot the coil\n",
|
||||||
|
"coil = bs.parse_coil(\"coils/\" + COIL_NAME + \".csv\")\n",
|
||||||
|
"coil = coil[0:3, :]\n",
|
||||||
|
"coil = coil * 10 + 20\n",
|
||||||
|
"plt.plot(coil[0, :], coil[1, :], color=\"black\", alpha=0.1)\n",
|
||||||
|
"\n",
|
||||||
"# plt.imshow(x_contour, cmap='hot', interpolation='nearest')\n",
|
"# plt.imshow(x_contour, cmap='hot', interpolation='nearest')\n",
|
||||||
"plt.contourf(\n",
|
"plt.contourf(\n",
|
||||||
" x_contour,\n",
|
" y_contour,\n",
|
||||||
" vmin=vmin,\n",
|
" vmin=vmin,\n",
|
||||||
" vmax=vmax,\n",
|
" vmax=vmax,\n",
|
||||||
" cmap=cm.magma,\n",
|
" cmap=cm.magma,\n",
|
||||||
" levels=50,\n",
|
" levels=50,\n",
|
||||||
")\n",
|
")\n",
|
||||||
"\n",
|
"\n",
|
||||||
"plt.savefig(COIL_NAME + \"_field_x.png\")\n",
|
"# hide the axis\n",
|
||||||
|
"plt.axis(\"off\")\n",
|
||||||
|
"\n",
|
||||||
|
"# remove any padding\n",
|
||||||
|
"plt.tight_layout()\n",
|
||||||
|
"plt.margins(x=0)\n",
|
||||||
|
"plt.margins(x=0)\n",
|
||||||
|
"\n",
|
||||||
|
"plt.savefig(COIL_NAME + \"_field_y.png\", bbox_inches=\"tight\", pad_inches=0)\n",
|
||||||
"\n",
|
"\n",
|
||||||
"print(vmin, vmax)\n",
|
"print(vmin, vmax)\n",
|
||||||
"\n",
|
"\n",
|
||||||
|
|
@ -350,7 +364,7 @@
|
||||||
},
|
},
|
||||||
"vscode": {
|
"vscode": {
|
||||||
"interpreter": {
|
"interpreter": {
|
||||||
"hash": "1ce20143987840b9786ebb5907032c9c3a8efacbb887dbb0ebc4934f2ad26cb3"
|
"hash": "fc384f9db26c31784edfba3761ba3d2c7b2f9b8a63e03a9eb0778fc35334efe1"
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
},
|
},
|
||||||
|
|
@ -1,111 +0,0 @@
|
||||||
import numpy as np
|
|
||||||
|
|
||||||
|
|
||||||
# get the point on an arc at the given angle
|
|
||||||
def get_arc_point(angle, radius):
|
|
||||||
return (
|
|
||||||
radius * np.cos(np.deg2rad(angle)),
|
|
||||||
radius * np.sin(np.deg2rad(angle)),
|
|
||||||
)
|
|
||||||
|
|
||||||
|
|
||||||
# draw an arc
|
|
||||||
def draw_arc(start_angle, end_angle, radius, step=5):
|
|
||||||
# make sure start_angle is less then end_angle
|
|
||||||
if start_angle > end_angle:
|
|
||||||
start_angle, end_angle = end_angle, start_angle
|
|
||||||
|
|
||||||
points = []
|
|
||||||
for angle in np.arange(start_angle, end_angle, step):
|
|
||||||
x = radius * np.cos(np.deg2rad(angle))
|
|
||||||
y = radius * np.sin(np.deg2rad(angle))
|
|
||||||
points.append((x, y))
|
|
||||||
if angle != end_angle:
|
|
||||||
x = radius * np.cos(np.deg2rad(end_angle))
|
|
||||||
y = radius * np.sin(np.deg2rad(end_angle))
|
|
||||||
points.append((x, y))
|
|
||||||
return points
|
|
||||||
|
|
||||||
|
|
||||||
# roate the points by the required angle
|
|
||||||
def rotate(points, angle):
|
|
||||||
return [
|
|
||||||
[
|
|
||||||
x * np.cos(np.deg2rad(angle)) - y * np.sin(np.deg2rad(angle)),
|
|
||||||
x * np.sin(np.deg2rad(angle)) + y * np.cos(np.deg2rad(angle)),
|
|
||||||
]
|
|
||||||
for x, y in points
|
|
||||||
]
|
|
||||||
|
|
||||||
def scale(points, scale):
|
|
||||||
return [[x * scale, y * scale] for x, y in points]
|
|
||||||
|
|
||||||
# rotate a point
|
|
||||||
def rotate_point(x, y, angle, ox=0, oy=0):
|
|
||||||
x -= ox
|
|
||||||
y -= oy
|
|
||||||
qx = x * np.cos(np.deg2rad(angle)) - y * np.sin(np.deg2rad(angle))
|
|
||||||
qy = x * np.sin(np.deg2rad(angle)) + y * np.cos(np.deg2rad(angle))
|
|
||||||
qx += ox
|
|
||||||
qy += oy
|
|
||||||
return qx, qy
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
# move the points out to the distance at the requited angle
|
|
||||||
def translate(points, distance, angle):
|
|
||||||
return [
|
|
||||||
[
|
|
||||||
x + distance * np.cos(np.deg2rad(angle)),
|
|
||||||
y + distance * np.sin(np.deg2rad(angle)),
|
|
||||||
]
|
|
||||||
for x, y in points
|
|
||||||
]
|
|
||||||
|
|
||||||
|
|
||||||
# flip the y coordinate
|
|
||||||
def flip_y(points):
|
|
||||||
return [[x, -y] for x, y in points]
|
|
||||||
|
|
||||||
|
|
||||||
def flip_x(points):
|
|
||||||
return [[-x, y] for x, y in points]
|
|
||||||
|
|
||||||
|
|
||||||
def optimize_points(points):
|
|
||||||
# follow the line and remove points that are in the same direction as the previous poin
|
|
||||||
# keep doing this until the direction changes significantly
|
|
||||||
# this is a very simple optimization that removes a lot of points
|
|
||||||
# it's not perfect but it's a good start
|
|
||||||
optimized_points = []
|
|
||||||
for i in range(len(points)):
|
|
||||||
if i == 0:
|
|
||||||
optimized_points.append(points[i])
|
|
||||||
else:
|
|
||||||
vector1 = np.array(points[i]) - np.array(points[i - 1])
|
|
||||||
vector2 = np.array(points[(i + 1) % len(points)]) - np.array(points[i])
|
|
||||||
length1 = np.linalg.norm(vector1)
|
|
||||||
length2 = np.linalg.norm(vector2)
|
|
||||||
if length1 > 0 and length2 > 0:
|
|
||||||
dot = np.dot(vector1, vector2) / (length1 * length2)
|
|
||||||
# clamp dot between -1 and 1
|
|
||||||
dot = max(-1, min(1, dot))
|
|
||||||
angle = np.arccos(dot)
|
|
||||||
if angle > np.deg2rad(5):
|
|
||||||
optimized_points.append(points[i])
|
|
||||||
print("Optimised from {} to {} points".format(len(points), len(optimized_points)))
|
|
||||||
return optimized_points
|
|
||||||
|
|
||||||
|
|
||||||
def chaikin(points, iterations):
|
|
||||||
if iterations == 0:
|
|
||||||
return points
|
|
||||||
l = len(points)
|
|
||||||
smoothed = []
|
|
||||||
for i in range(l - 1):
|
|
||||||
x1, y1 = points[i]
|
|
||||||
x2, y2 = points[i + 1]
|
|
||||||
smoothed.append([0.95 * x1 + 0.05 * x2, 0.95 * y1 + 0.05 * y2])
|
|
||||||
smoothed.append([0.05 * x1 + 0.95 * x2, 0.05 * y1 + 0.95 * y2])
|
|
||||||
smoothed.append(points[l - 1])
|
|
||||||
return chaikin(smoothed, iterations - 1)
|
|
||||||
|
|
@ -1,607 +0,0 @@
|
||||||
{
|
|
||||||
"cells": [
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"import numpy as np\n",
|
|
||||||
"import matplotlib.pyplot as plt\n",
|
|
||||||
"from mpl_toolkits.mplot3d import Axes3D\n",
|
|
||||||
"from biot_savart_v4_3 import parse_coil, plot_coil, slice_coil\n",
|
|
||||||
"from tqdm.notebook import trange, tqdm"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"# load up the simple spiral coil\n",
|
|
||||||
"coil1 = parse_coil(\"coils/coil_12_spiral.csv\")\n",
|
|
||||||
"plot_coil(\"coils/coil_12_spiral.csv\")\n",
|
|
||||||
"coil1 = slice_coil(coil1, 1)\n",
|
|
||||||
"coil1 = coil1.T\n",
|
|
||||||
"print(coil1.shape)\n",
|
|
||||||
"\n",
|
|
||||||
"coil2 = parse_coil(\"coils/coil_12_custom.csv\")\n",
|
|
||||||
"plot_coil(\"coils/coil_12_custom.csv\")\n",
|
|
||||||
"coil2 = slice_coil(coil2, 1)\n",
|
|
||||||
"coil2 = coil2.T\n",
|
|
||||||
"print(coil2.shape)"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "markdown",
|
|
||||||
"metadata": {},
|
|
||||||
"source": [
|
|
||||||
"# Simple Simulation of a dipole magnet"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"# magnetic field at a point x,y,z of a dipole magnet with moment m in the z direction\n",
|
|
||||||
"def B_(x, y, z, m=0.185):\n",
|
|
||||||
" mu0 = 4 * np.pi * 1e-7\n",
|
|
||||||
" r = np.sqrt(x**2 + y**2 + z**2)\n",
|
|
||||||
" return (\n",
|
|
||||||
" np.array(\n",
|
|
||||||
" [3 * x * z / r**5, 3 * y * z / r**5, (3 * z**2 / r**5 - 1 / r**3)]\n",
|
|
||||||
" )\n",
|
|
||||||
" * m\n",
|
|
||||||
" * mu0\n",
|
|
||||||
" )\n",
|
|
||||||
"\n",
|
|
||||||
"\n",
|
|
||||||
"def B(x, y, z, m=0.185, l=0.003, d=0.01):\n",
|
|
||||||
" d = d * 0.75\n",
|
|
||||||
" # simulate multiple points in the cylinder and add them together to create a disk of field\n",
|
|
||||||
" bx = 0\n",
|
|
||||||
" by = 0\n",
|
|
||||||
" bz = 0\n",
|
|
||||||
" for degrees in range(0, 360, 30):\n",
|
|
||||||
" angle = np.deg2rad(degrees)\n",
|
|
||||||
" x_, y_, z_ = B_(\n",
|
|
||||||
" x + d / 2 * np.cos(angle),\n",
|
|
||||||
" y + d / 2 * np.sin(angle),\n",
|
|
||||||
" z - l / 2,\n",
|
|
||||||
" m / (360 / 60),\n",
|
|
||||||
" )\n",
|
|
||||||
" bx += x_\n",
|
|
||||||
" by += y_\n",
|
|
||||||
" bz += z_\n",
|
|
||||||
" x_, y_, z_ = B_(\n",
|
|
||||||
" x + d / 2 * np.cos(angle),\n",
|
|
||||||
" y + d / 2 * np.sin(angle),\n",
|
|
||||||
" z + l / 2,\n",
|
|
||||||
" m / (360 / 60),\n",
|
|
||||||
" )\n",
|
|
||||||
" bx += x_\n",
|
|
||||||
" by += y_\n",
|
|
||||||
" bz += z_\n",
|
|
||||||
" return np.array([bx, by, bz])\n",
|
|
||||||
"\n",
|
|
||||||
"\n",
|
|
||||||
"def plot_field_slice(x, y, bx, by, mag, name=\"magnetic_field.png\"):\n",
|
|
||||||
" # plot the magnetic field\n",
|
|
||||||
" fig = plt.figure()\n",
|
|
||||||
" ax = fig.add_subplot(111)\n",
|
|
||||||
" ax.streamplot(\n",
|
|
||||||
" x,\n",
|
|
||||||
" y,\n",
|
|
||||||
" bx,\n",
|
|
||||||
" by,\n",
|
|
||||||
" linewidth=1,\n",
|
|
||||||
" cmap=plt.cm.inferno,\n",
|
|
||||||
" density=2,\n",
|
|
||||||
" arrowstyle=\"->\",\n",
|
|
||||||
" arrowsize=1.5,\n",
|
|
||||||
" )\n",
|
|
||||||
"\n",
|
|
||||||
" ax.set_xlabel(\"$x$\")\n",
|
|
||||||
" ax.set_ylabel(\"$y$\")\n",
|
|
||||||
" ax.set_xlim(-0.1, 0.1)\n",
|
|
||||||
" ax.set_ylim(-0.1, 0.1)\n",
|
|
||||||
" ax.set_aspect(\"equal\")\n",
|
|
||||||
"\n",
|
|
||||||
" # plot the magniture of the field as an image\n",
|
|
||||||
" im = ax.imshow(\n",
|
|
||||||
" mag, extent=[-0.1, 0.1, -0.1, 0.1], origin=\"lower\", cmap=plt.cm.inferno\n",
|
|
||||||
" )\n",
|
|
||||||
"\n",
|
|
||||||
" # draw the magnet\n",
|
|
||||||
" ax.add_patch(plt.Rectangle((-0.005, -0.0015), 0.01, 0.003, fc=\"w\", ec=\"k\", lw=1))\n",
|
|
||||||
"\n",
|
|
||||||
" fig.show()\n",
|
|
||||||
" # save the figure\n",
|
|
||||||
" fig.savefig(name)\n",
|
|
||||||
"\n",
|
|
||||||
"\n",
|
|
||||||
"# # calculate the magnetic field at y = 0, over z = -1, 1 and x = -1, 1\n",
|
|
||||||
"x = np.linspace(-0.1, 0.1, 100)\n",
|
|
||||||
"z = np.linspace(-0.1, 0.1, 100)\n",
|
|
||||||
"X, Z = np.meshgrid(x, z)\n",
|
|
||||||
"Bx, By, Bz = B(X, 0, Z)\n",
|
|
||||||
"\n",
|
|
||||||
"print(Bx.shape, By.shape, Bz.shape)\n",
|
|
||||||
"\n",
|
|
||||||
"plot_field_slice(\n",
|
|
||||||
" X,\n",
|
|
||||||
" Z,\n",
|
|
||||||
" Bx,\n",
|
|
||||||
" Bz,\n",
|
|
||||||
" np.log(np.sqrt(Bx**2 + By**2 + Bz**2)),\n",
|
|
||||||
" \"magnetic_field_side.png\",\n",
|
|
||||||
")\n",
|
|
||||||
"\n",
|
|
||||||
"# # calculate the magnetic field at z = 1, over y = -1, 1 and x = -1, 1\n",
|
|
||||||
"x = np.linspace(-0.1, 0.1, 100)\n",
|
|
||||||
"y = np.linspace(-0.1, 0.1, 100)\n",
|
|
||||||
"X, Y = np.meshgrid(x, y)\n",
|
|
||||||
"Bx, By, Bz = B(X, Y, 0.01)\n",
|
|
||||||
"\n",
|
|
||||||
"plot_field_slice(\n",
|
|
||||||
" X,\n",
|
|
||||||
" Y,\n",
|
|
||||||
" Bx,\n",
|
|
||||||
" By,\n",
|
|
||||||
" np.log(np.sqrt(Bx**2 + By**2 + Bz**2)),\n",
|
|
||||||
" \"magnetic_field_bottom.png\",\n",
|
|
||||||
")\n",
|
|
||||||
"\n",
|
|
||||||
"# calculate the magnetic field in a 3d volume\n",
|
|
||||||
"x = np.linspace(-0.1, 0.1, 100)\n",
|
|
||||||
"y = np.linspace(-0.1, 0.1, 100)\n",
|
|
||||||
"z = np.linspace(-0.1, 0.1, 100)\n",
|
|
||||||
"X, Y, Z = np.meshgrid(x, y, z)\n",
|
|
||||||
"Bx, By, Bz = B(X, Y, Z)"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"# do a 3d quivwer plot of the magnetic field\n",
|
|
||||||
"fig = plt.figure()\n",
|
|
||||||
"ax = fig.add_subplot(111, projection=\"3d\")\n",
|
|
||||||
"# down sample the results\n",
|
|
||||||
"ax.quiver(\n",
|
|
||||||
" X[::10, ::10, ::10],\n",
|
|
||||||
" Y[::10, ::10, ::10],\n",
|
|
||||||
" Z[::10, ::10, ::10],\n",
|
|
||||||
" Bx[::10, ::10, ::10],\n",
|
|
||||||
" By[::10, ::10, ::10],\n",
|
|
||||||
" Bz[::10, ::10, ::10],\n",
|
|
||||||
" length=0.01,\n",
|
|
||||||
" normalize=True,\n",
|
|
||||||
")\n",
|
|
||||||
"ax.set_xlabel(\"$x$\")\n",
|
|
||||||
"ax.set_ylabel(\"$y$\")\n",
|
|
||||||
"ax.set_zlabel(\"$z$\")\n",
|
|
||||||
"ax.set_xlim(-0.1, 0.1)\n",
|
|
||||||
"ax.set_ylim(-0.1, 0.1)\n",
|
|
||||||
"ax.set_zlim(-0.1, 0.1)\n",
|
|
||||||
"ax.set_aspect(\"equal\")\n",
|
|
||||||
"# make the plot larger\n",
|
|
||||||
"fig.set_size_inches(10, 10)\n",
|
|
||||||
"fig.show()"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"# use mayavi to plot the magnetic field\n",
|
|
||||||
"from mayavi import mlab\n",
|
|
||||||
"\n",
|
|
||||||
"mlab.figure(bgcolor=(1, 1, 1), size=(800, 800))\n",
|
|
||||||
"mlab.quiver3d(\n",
|
|
||||||
" X,\n",
|
|
||||||
" Y,\n",
|
|
||||||
" Z,\n",
|
|
||||||
" Bx,\n",
|
|
||||||
" By,\n",
|
|
||||||
" Bz,\n",
|
|
||||||
" mode=\"2ddash\",\n",
|
|
||||||
" scale_factor=1,\n",
|
|
||||||
" scale_mode=\"none\",\n",
|
|
||||||
" color=(0, 0, 0),\n",
|
|
||||||
")\n",
|
|
||||||
"mlab.axes(\n",
|
|
||||||
" xlabel=\"$x$\", ylabel=\"$y$\", zlabel=\"$z$\", ranges=[-0.1, 0.1, -0.1, 0.1, -0.1, 0.1]\n",
|
|
||||||
")\n",
|
|
||||||
"mlab.show()"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"# use myavi to plot a contour plot of the magnetic field\n",
|
|
||||||
"from mayavi import mlab\n",
|
|
||||||
"\n",
|
|
||||||
"mlab.figure(bgcolor=(1, 1, 1), size=(800, 800))\n",
|
|
||||||
"mlab.contour3d(\n",
|
|
||||||
" X, Y, Z, np.sqrt(Bx**2 + By**2 + Bz**2), contours=20, transparent=True\n",
|
|
||||||
")\n",
|
|
||||||
"mlab.axes(\n",
|
|
||||||
" xlabel=\"$x$\", ylabel=\"$y$\", zlabel=\"$z$\", ranges=[-0.1, 0.1, -0.1, 0.1, -0.1, 0.1]\n",
|
|
||||||
")\n",
|
|
||||||
"mlab.show()"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"# calculate the force on a wire of length l carrying current I at a point x,y,z with a direction vector d\n",
|
|
||||||
"def F(p, d, I, l):\n",
|
|
||||||
" return I * l * np.cross(d, B(p[0], p[1], p[2]))\n",
|
|
||||||
"\n",
|
|
||||||
"\n",
|
|
||||||
"def calculate_forces_on_wire_points(points):\n",
|
|
||||||
" Fx = []\n",
|
|
||||||
" Fy = []\n",
|
|
||||||
" Fz = []\n",
|
|
||||||
" # calculate the force on each point\n",
|
|
||||||
" for i in range(len(points)):\n",
|
|
||||||
" # calculate the direction vector\n",
|
|
||||||
" dx = points[i][0] - points[(i + 1) % len(points)][0]\n",
|
|
||||||
" dy = points[i][1] - points[(i + 1) % len(points)][1]\n",
|
|
||||||
" dz = points[i][2] - points[(i + 1) % len(points)][2]\n",
|
|
||||||
" d = np.array([dx, dy, dz])\n",
|
|
||||||
" # get the length of d\n",
|
|
||||||
" l = np.sqrt(dx**2 + dy**2 + dz**2)\n",
|
|
||||||
" if l > 0:\n",
|
|
||||||
" # normalise d\n",
|
|
||||||
" d = d / l\n",
|
|
||||||
" # calculate the force\n",
|
|
||||||
" fx, fy, fz = F(points[i], d, points[i][3], l)\n",
|
|
||||||
" Fx.append(fx)\n",
|
|
||||||
" Fy.append(fy)\n",
|
|
||||||
" Fz.append(fz)\n",
|
|
||||||
" else:\n",
|
|
||||||
" Fx.append(0)\n",
|
|
||||||
" Fy.append(0)\n",
|
|
||||||
" Fz.append(0)\n",
|
|
||||||
" return Fx, Fy, Fz\n",
|
|
||||||
"\n",
|
|
||||||
"\n",
|
|
||||||
"# locate the loop of wire directly below the magnet\n",
|
|
||||||
"x = 0.01\n",
|
|
||||||
"y = 0\n",
|
|
||||||
"z = -0.002\n",
|
|
||||||
"r1 = 0.001\n",
|
|
||||||
"r2 = 0.01\n",
|
|
||||||
"\n",
|
|
||||||
"\n",
|
|
||||||
"points = coil2.copy()\n",
|
|
||||||
"# scale the points from mm to m\n",
|
|
||||||
"# shift the coil to the correct position\n",
|
|
||||||
"for i in range(len(points)):\n",
|
|
||||||
" points[i][0] = points[i][0] / 1000 + x\n",
|
|
||||||
" points[i][1] = points[i][1] / 1000 + y\n",
|
|
||||||
" points[i][2] = points[i][2] / 1000 + z\n",
|
|
||||||
"\n",
|
|
||||||
"\n",
|
|
||||||
"# plot the points in 2D x,y\n",
|
|
||||||
"fig = plt.figure()\n",
|
|
||||||
"ax = fig.add_subplot(111)\n",
|
|
||||||
"ax.plot([p[0] for p in points], [p[1] for p in points])\n",
|
|
||||||
"ax.set_xlabel(\"$x$\")\n",
|
|
||||||
"ax.set_ylabel(\"$y$\")\n",
|
|
||||||
"ax.set_xlim(-0.01 + x, 0.01 + x)\n",
|
|
||||||
"ax.set_ylim(-0.01 + y, 0.01 + y)\n",
|
|
||||||
"ax.set_aspect(\"equal\")\n",
|
|
||||||
"fig.show()\n",
|
|
||||||
"\n",
|
|
||||||
"Fx, Fy, Fz = calculate_forces_on_wire_points(points)\n",
|
|
||||||
"\n",
|
|
||||||
"\n",
|
|
||||||
"# plot the wire along with arrows showing the force\n",
|
|
||||||
"fig = plt.figure()\n",
|
|
||||||
"ax = fig.add_subplot(111, projection=\"3d\")\n",
|
|
||||||
"ax.plot([p[0] for p in points], [p[1] for p in points], [p[2] for p in points])\n",
|
|
||||||
"\n",
|
|
||||||
"print(sum(Fx) / 9.8)\n",
|
|
||||||
"\n",
|
|
||||||
"# subsample the points and force vectors to make the plot clearer\n",
|
|
||||||
"points = points[::500]\n",
|
|
||||||
"Fx = Fx[::500]\n",
|
|
||||||
"Fy = Fy[::500]\n",
|
|
||||||
"Fz = Fz[::500]\n",
|
|
||||||
"\n",
|
|
||||||
"ax.quiver(\n",
|
|
||||||
" [p[0] for p in points],\n",
|
|
||||||
" [p[1] for p in points],\n",
|
|
||||||
" [p[2] for p in points],\n",
|
|
||||||
" np.sqrt(Fx),\n",
|
|
||||||
" 0,\n",
|
|
||||||
" 0,\n",
|
|
||||||
" length=0.01,\n",
|
|
||||||
" normalize=True,\n",
|
|
||||||
")\n",
|
|
||||||
"ax.set_xlabel(\"$x$\")\n",
|
|
||||||
"ax.set_ylabel(\"$y$\")\n",
|
|
||||||
"ax.set_zlabel(\"$z$\")\n",
|
|
||||||
"ax.set_xlim(x - 0.02, x + 0.02)\n",
|
|
||||||
"ax.set_ylim(y - 0.02, y + 0.02)\n",
|
|
||||||
"ax.set_zlim(z - 0.02, z + 0.02)\n",
|
|
||||||
"ax.set_aspect(\"equal\")\n",
|
|
||||||
"\n",
|
|
||||||
"# change the figure size\n",
|
|
||||||
"fig.set_size_inches(10, 10)\n",
|
|
||||||
"\n",
|
|
||||||
"fig.show()\n",
|
|
||||||
"# coil2 = 0.0375037573536258\n",
|
|
||||||
"# coil1 = 0.010254238165764389"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"def sweep_coil(coil, X):\n",
|
|
||||||
" Y = np.zeros(len(X))\n",
|
|
||||||
" Z = -0.01 * np.ones(len(X))\n",
|
|
||||||
"\n",
|
|
||||||
" # loop through the locations and calculate the forces, sum up the force in the X direction for each location\n",
|
|
||||||
" Fx = []\n",
|
|
||||||
" Fy = []\n",
|
|
||||||
" Fz = []\n",
|
|
||||||
" for p in trange(len(X)):\n",
|
|
||||||
" points = coil.copy()\n",
|
|
||||||
" # scale the points from mm to m\n",
|
|
||||||
" # shift the coil to the correct position\n",
|
|
||||||
" for i in range(len(points)):\n",
|
|
||||||
" points[i][0] = points[i][0] / 1000 + X[p]\n",
|
|
||||||
" points[i][1] = points[i][1] / 1000 + Y[p]\n",
|
|
||||||
" points[i][2] = points[i][2] / 1000 + Z[p]\n",
|
|
||||||
" Fx_, Fy_, Fz_ = calculate_forces_on_wire_points(points)\n",
|
|
||||||
" Fx.append(sum(Fx_))\n",
|
|
||||||
" Fy.append(sum(Fy_))\n",
|
|
||||||
" Fz.append(sum(Fz_))\n",
|
|
||||||
" return Fx, Fy, Fz\n",
|
|
||||||
"\n",
|
|
||||||
"\n",
|
|
||||||
"# sweep the coild from -3cm to 3cm in 0.01m steps\n",
|
|
||||||
"X = np.linspace(-0.03, 0.03, 100)\n",
|
|
||||||
"Fx_1_straight, Fy_1_straight, Fz_1_straight = sweep_coil(coil1, X)\n",
|
|
||||||
"Fx_2_straight, Fy_2_straight, Fz_2_straight = sweep_coil(coil2, X)"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"# plot the force as a function of x\n",
|
|
||||||
"plt.plot(X, -np.array(Fx_1_straight) / 9.8, label=\"coil1\", color=\"red\")\n",
|
|
||||||
"plt.plot(X, np.array(Fx_2_straight) / 9.8, label=\"coil2\", color=\"blue\")\n",
|
|
||||||
"# plot a dotted line along y = 0\n",
|
|
||||||
"plt.plot([X[0], X[-1]], [0, 0], \"--\")"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"# plot the force as a function of x\n",
|
|
||||||
"plt.plot(X, -np.array(Fz_1_straight) / 9.8, label=\"coil1\", color=\"red\")\n",
|
|
||||||
"plt.plot(X, np.array(Fz_2_straight) / 9.8, label=\"coil2\", color=\"blue\")\n",
|
|
||||||
"# plot a dotted line along y = 0\n",
|
|
||||||
"plt.plot([X[0], X[-1]], [0, 0], \"--\")"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"# instead of sweeping horizontally, we'll sweep the coils around a circle\n",
|
|
||||||
"def sweep_coil_cirlc(coil, coil_center_radius, theta):\n",
|
|
||||||
" X = coil_center_radius * np.cos(np.deg2rad(theta))\n",
|
|
||||||
" Y = coil_center_radius * np.sin(np.deg2rad(theta))\n",
|
|
||||||
" Z = -0.01 * np.ones(100)\n",
|
|
||||||
"\n",
|
|
||||||
" # loop through the locations and calculate the forces, sum up the force in the X direction for each location\n",
|
|
||||||
" Torque = []\n",
|
|
||||||
" Fx = []\n",
|
|
||||||
" Fy = []\n",
|
|
||||||
" Fz = []\n",
|
|
||||||
" for p in trange(len(theta)):\n",
|
|
||||||
" angle = np.deg2rad(theta[p]) - np.pi / 2\n",
|
|
||||||
" x = X[p]\n",
|
|
||||||
" y = Y[p]\n",
|
|
||||||
" z = Z[p]\n",
|
|
||||||
"\n",
|
|
||||||
" points = coil.copy()\n",
|
|
||||||
" for i in range(len(points)):\n",
|
|
||||||
" px = points[i][0] / 1000\n",
|
|
||||||
" py = points[i][1] / 1000\n",
|
|
||||||
" pz = points[i][2] / 1000\n",
|
|
||||||
" # rotate the points so the coil is correctly oriented\n",
|
|
||||||
" points[i][0] = px * np.cos(angle) - py * np.sin(angle) + x\n",
|
|
||||||
" points[i][1] = (\n",
|
|
||||||
" px * np.sin(angle) + py * np.cos(angle) + y - coil_center_radius\n",
|
|
||||||
" )\n",
|
|
||||||
" points[i][2] = pz + z\n",
|
|
||||||
" # feel the force\n",
|
|
||||||
" Fx_, Fy_, Fz_ = calculate_forces_on_wire_points(points)\n",
|
|
||||||
" Fx.append(sum(Fx_))\n",
|
|
||||||
" Fy.append(sum(Fy_))\n",
|
|
||||||
" Fz.append(sum(Fz_))\n",
|
|
||||||
" # calculate the torque - which should be 90 degress to the angle\n",
|
|
||||||
" torque_angle = np.deg2rad(theta[p] - 90)\n",
|
|
||||||
" Torque.append(sum(Fx_) * np.cos(torque_angle) + sum(Fy_) * np.sin(torque_angle))\n",
|
|
||||||
"\n",
|
|
||||||
" return Fx, Fy, Fz, Torque\n",
|
|
||||||
"\n",
|
|
||||||
"\n",
|
|
||||||
"# sweep the coils from -45 to 45 degrees in 1 degree steps\n",
|
|
||||||
"theta = np.linspace(0, 180, 100)\n",
|
|
||||||
"Fx_1_curve, Fy_1_curve, Fz_1_curve, Torque_1 = sweep_coil_cirlc(\n",
|
|
||||||
" coil1, 20.5 / 1000, theta\n",
|
|
||||||
")\n",
|
|
||||||
"Fx_2_curve, Fy_2_curve, Fz_2_curve, Torque_2 = sweep_coil_cirlc(\n",
|
|
||||||
" coil2, 19.5 / 1000, theta\n",
|
|
||||||
")"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"plt.plot(\n",
|
|
||||||
" theta, -np.array(Fx_1_curve) + np.array(Fx_1_curve), label=\"coil1\", color=\"red\"\n",
|
|
||||||
")\n",
|
|
||||||
"plt.plot(theta, np.array(Torque_2) - np.array(Fx_2_curve), label=\"coil2\", color=\"blue\")\n",
|
|
||||||
"# plot a dotted line along y = 0\n",
|
|
||||||
"plt.plot([theta[0], theta[-1]], [0, 0], \"--\")"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"plt.plot(theta, -np.array(Torque_1) / 9.8, label=\"coil1\", color=\"red\")\n",
|
|
||||||
"plt.plot(theta, np.array(Torque_2) / 9.8, label=\"coil2\", color=\"blue\")\n",
|
|
||||||
"# plot a dotted line along y = 0\n",
|
|
||||||
"plt.plot([theta[0], theta[-1]], [0, 0], \"--\")"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"# plot arrows for the Fx and Fy components\n",
|
|
||||||
"X = 19.5 * np.cos(np.deg2rad(theta)) / 1000\n",
|
|
||||||
"Y = 19.5 * np.sin(np.deg2rad(theta)) / 1000\n",
|
|
||||||
"plt.quiver(X[::5], Y[::5], Fx_1_curve[::5], Fy_1_curve[::5], color=\"red\")\n",
|
|
||||||
"# make the axis equal so the arrows are not stretched\n",
|
|
||||||
"plt.axis(\"equal\")"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"plt.plot(theta, -np.array(Fz_1_curve) / 9.8, label=\"coil1\", color=\"red\")\n",
|
|
||||||
"plt.plot(theta, np.array(Fz_2_curve) / 9.8, label=\"coil2\", color=\"blue\")\n",
|
|
||||||
"# plot a dotted line along y = 0\n",
|
|
||||||
"plt.plot([theta[0], theta[-1]], [0, 0], \"--\")"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"# instead of sweeping horizontally, we'll sweep the coils around a circle\n",
|
|
||||||
"def sweep_coil_cirlc(coil, coil_center_radius):\n",
|
|
||||||
" # sweep the coils from -45 to 45 degrees in 1 degree steps\n",
|
|
||||||
" theta = np.linspace(0, 180, 20)\n",
|
|
||||||
" X = coil_center_radius * np.cos(np.deg2rad(theta))\n",
|
|
||||||
" Y = coil_center_radius * np.sin(np.deg2rad(theta))\n",
|
|
||||||
" Z = 1 * np.ones(100)\n",
|
|
||||||
"\n",
|
|
||||||
" # loop through the locations and calculate the forces, sum up the force in the X direction for each location\n",
|
|
||||||
" Fx = []\n",
|
|
||||||
" Fy = []\n",
|
|
||||||
" Fz = []\n",
|
|
||||||
" for p in trange(len(theta)):\n",
|
|
||||||
" angle = np.deg2rad(theta[p]) - np.pi / 2\n",
|
|
||||||
" x = X[p]\n",
|
|
||||||
" y = Y[p]\n",
|
|
||||||
" z = Z[p]\n",
|
|
||||||
"\n",
|
|
||||||
" points = coil.copy()\n",
|
|
||||||
" for i in range(len(points)):\n",
|
|
||||||
" px = points[i][0] / 1000\n",
|
|
||||||
" py = points[i][1] / 1000\n",
|
|
||||||
" pz = points[i][2] / 1000\n",
|
|
||||||
" # rotate the points so the coil is correctly oriented\n",
|
|
||||||
" points[i][0] = px * np.cos(angle) - py * np.sin(angle) + x\n",
|
|
||||||
" points[i][1] = (\n",
|
|
||||||
" px * np.sin(angle) + py * np.cos(angle) + y - coil_center_radius\n",
|
|
||||||
" )\n",
|
|
||||||
" points[i][2] = pz + z\n",
|
|
||||||
" plt.plot([p[0] for p in points], [p[1] for p in points], linewidth=0.5)\n",
|
|
||||||
" # add the torque arrow to the plot\n",
|
|
||||||
" torque_angle = np.deg2rad(theta[p] - 90)\n",
|
|
||||||
" torque_x1 = x\n",
|
|
||||||
" torque_y1 = y\n",
|
|
||||||
" torque_x2 = x + 0.01 * np.cos(torque_angle)\n",
|
|
||||||
" torque_y2 = y + 0.01 * np.sin(torque_angle)\n",
|
|
||||||
" plt.arrow(\n",
|
|
||||||
" torque_x1,\n",
|
|
||||||
" torque_y1,\n",
|
|
||||||
" torque_x2 - torque_x1,\n",
|
|
||||||
" torque_y2 - torque_y1,\n",
|
|
||||||
" head_width=0.001,\n",
|
|
||||||
" head_length=0.002,\n",
|
|
||||||
" fc=\"k\",\n",
|
|
||||||
" ec=\"k\",\n",
|
|
||||||
" )\n",
|
|
||||||
"\n",
|
|
||||||
"\n",
|
|
||||||
"sweep_coil_cirlc(coil2, 19.5 / 1000)"
|
|
||||||
]
|
|
||||||
}
|
|
||||||
],
|
|
||||||
"metadata": {
|
|
||||||
"kernelspec": {
|
|
||||||
"display_name": "Python 3.10.7 ('venv': venv)",
|
|
||||||
"language": "python",
|
|
||||||
"name": "python3"
|
|
||||||
},
|
|
||||||
"language_info": {
|
|
||||||
"codemirror_mode": {
|
|
||||||
"name": "ipython",
|
|
||||||
"version": 3
|
|
||||||
},
|
|
||||||
"file_extension": ".py",
|
|
||||||
"mimetype": "text/x-python",
|
|
||||||
"name": "python",
|
|
||||||
"nbconvert_exporter": "python",
|
|
||||||
"pygments_lexer": "ipython3",
|
|
||||||
"version": "3.10.7"
|
|
||||||
},
|
|
||||||
"vscode": {
|
|
||||||
"interpreter": {
|
|
||||||
"hash": "1ce20143987840b9786ebb5907032c9c3a8efacbb887dbb0ebc4934f2ad26cb3"
|
|
||||||
}
|
|
||||||
}
|
|
||||||
},
|
|
||||||
"nbformat": 4,
|
|
||||||
"nbformat_minor": 2
|
|
||||||
}
|
|
||||||
|
|
@ -1,525 +0,0 @@
|
||||||
{
|
|
||||||
"cells": [
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"import numpy as np\n",
|
|
||||||
"import matplotlib.pyplot as plt\n",
|
|
||||||
"from mpl_toolkits.mplot3d import Axes3D\n",
|
|
||||||
"from biot_savart_v4_3 import parse_coil, plot_coil, slice_coil, plot_coil2\n",
|
|
||||||
"from tqdm.notebook import trange, tqdm\n",
|
|
||||||
"from helpers import get_arc_point, draw_arc, rotate, scale, rotate_point, translate"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"# Track width and spacing\n",
|
|
||||||
"TRACK_WIDTH = 0.127\n",
|
|
||||||
"TRACK_SPACING = 0.127\n",
|
|
||||||
"\n",
|
|
||||||
"# via defaults\n",
|
|
||||||
"VIA_DIAM = 0.8\n",
|
|
||||||
"VIA_DRILL = 0.4\n",
|
|
||||||
"\n",
|
|
||||||
"# this is for a 1.27mm pitch pin\n",
|
|
||||||
"PIN_DIAM = 1.0\n",
|
|
||||||
"PIN_DRILL = 0.65\n",
|
|
||||||
"\n",
|
|
||||||
"# this is for the PCB connector - see https://www.farnell.com/datasheets/2003059.pdf\n",
|
|
||||||
"PAD_WIDTH = 3\n",
|
|
||||||
"PAD_HEIGHT = 2\n",
|
|
||||||
"PAD_PITCH = 2.5\n",
|
|
||||||
"\n",
|
|
||||||
"STATOR_HOLE_RADIUS = 5.5\n",
|
|
||||||
"HOLE_SPACING = 0.25\n",
|
|
||||||
"\n",
|
|
||||||
"# PCB Edge size\n",
|
|
||||||
"STATOR_RADIUS = 30\n",
|
|
||||||
"\n",
|
|
||||||
"# where to puth the mounting pins\n",
|
|
||||||
"SCREW_HOLE_DRILL_DIAM = 2.3 # 2.3mm drill for a 2mm screw\n",
|
|
||||||
"SCREW_HOLE_RADIUS = STATOR_RADIUS\n",
|
|
||||||
"\n",
|
|
||||||
"# Coil params\n",
|
|
||||||
"TURNS = 16\n",
|
|
||||||
"COIL_CENTER_RADIUS = 19.95\n",
|
|
||||||
"COIL_VIA_RADIUS = 20.95\n",
|
|
||||||
"\n",
|
|
||||||
"\n",
|
|
||||||
"def get_points(spacing, inner_radius, outer_radius, start_angle, end_angle):\n",
|
|
||||||
" # first calculate the angle step size from the spacing and the inner_radius\n",
|
|
||||||
" spacing_angle = np.rad2deg(np.arctan2(spacing, inner_radius))\n",
|
|
||||||
" print(spacing_angle)\n",
|
|
||||||
" # now calculate the points be iterating from start_angle to end_angle with the spacing_angle\n",
|
|
||||||
" points = []\n",
|
|
||||||
" for angle in np.arange(start_angle, end_angle, spacing_angle * 2):\n",
|
|
||||||
" p1 = get_arc_point(angle, inner_radius)\n",
|
|
||||||
" points.append([p1[0], p1[1], 0, 0.5])\n",
|
|
||||||
" p2 = get_arc_point(angle + spacing_angle, outer_radius)\n",
|
|
||||||
" points.append([p2[0], p2[1], 0, 0.5])\n",
|
|
||||||
" points.append([p2[0], p2[1], -0.8, 0.5])\n",
|
|
||||||
" points.append([p1[0], p1[1], -0.8, 0.5])\n",
|
|
||||||
" return points\n",
|
|
||||||
"\n",
|
|
||||||
"\n",
|
|
||||||
"coil1 = get_points(\n",
|
|
||||||
" TRACK_SPACING + TRACK_WIDTH,\n",
|
|
||||||
" (COIL_CENTER_RADIUS - 10),\n",
|
|
||||||
" COIL_CENTER_RADIUS + 10,\n",
|
|
||||||
" -15,\n",
|
|
||||||
" 15,\n",
|
|
||||||
")\n",
|
|
||||||
"# move all the points in by COIL_CENTER_RADIUS\n",
|
|
||||||
"for p in coil1:\n",
|
|
||||||
" p[0] -= COIL_CENTER_RADIUS\n",
|
|
||||||
"# rotate the points by 90 degrees\n",
|
|
||||||
"for p in coil1:\n",
|
|
||||||
" p[0], p[1] = rotate_point(p[0], p[1], 90)\n",
|
|
||||||
"coil1 = np.array(coil1).T\n",
|
|
||||||
"plot_coil2(coil1)\n",
|
|
||||||
"coil1 = slice_coil(coil1, 0.1)\n",
|
|
||||||
"coil1 = coil1.T\n",
|
|
||||||
"print(coil1.shape)\n",
|
|
||||||
"print(len(coil1))"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "markdown",
|
|
||||||
"metadata": {},
|
|
||||||
"source": [
|
|
||||||
"# Simple Simulation of a dipole magnet"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"# magnetic field at a point x,y,z of a dipole magnet with moment m in the z direction\n",
|
|
||||||
"def B(x, y, z, m=0.185):\n",
|
|
||||||
" mu0 = 4 * np.pi * 1e-7\n",
|
|
||||||
" r = np.sqrt(x**2 + y**2 + z**2)\n",
|
|
||||||
" return (\n",
|
|
||||||
" np.array(\n",
|
|
||||||
" [3 * x * z / r**5, 3 * y * z / r**5, (3 * z**2 / r**5 - 1 / r**3)]\n",
|
|
||||||
" )\n",
|
|
||||||
" * m\n",
|
|
||||||
" * mu0\n",
|
|
||||||
" )\n",
|
|
||||||
"\n",
|
|
||||||
"\n",
|
|
||||||
"def plot_field_slice(x, y, bx, by, mag, name=\"magnetic_field.png\"):\n",
|
|
||||||
" # plot the magnetic field\n",
|
|
||||||
" fig = plt.figure()\n",
|
|
||||||
" ax = fig.add_subplot(111)\n",
|
|
||||||
" ax.streamplot(\n",
|
|
||||||
" x,\n",
|
|
||||||
" y,\n",
|
|
||||||
" bx,\n",
|
|
||||||
" by,\n",
|
|
||||||
" linewidth=1,\n",
|
|
||||||
" cmap=plt.cm.inferno,\n",
|
|
||||||
" density=2,\n",
|
|
||||||
" arrowstyle=\"->\",\n",
|
|
||||||
" arrowsize=1.5,\n",
|
|
||||||
" )\n",
|
|
||||||
"\n",
|
|
||||||
" ax.set_xlabel(\"$x$\")\n",
|
|
||||||
" ax.set_ylabel(\"$y$\")\n",
|
|
||||||
" ax.set_xlim(-0.1, 0.1)\n",
|
|
||||||
" ax.set_ylim(-0.1, 0.1)\n",
|
|
||||||
" ax.set_aspect(\"equal\")\n",
|
|
||||||
"\n",
|
|
||||||
" # plot the magniture of the field as an image\n",
|
|
||||||
" im = ax.imshow(\n",
|
|
||||||
" mag, extent=[-0.1, 0.1, -0.1, 0.1], origin=\"lower\", cmap=plt.cm.inferno\n",
|
|
||||||
" )\n",
|
|
||||||
"\n",
|
|
||||||
" fig.show()\n",
|
|
||||||
" # save the figure\n",
|
|
||||||
" fig.savefig(name)\n",
|
|
||||||
"\n",
|
|
||||||
"\n",
|
|
||||||
"# # calculate the magnetic field at y = 0, over z = -1, 1 and x = -1, 1\n",
|
|
||||||
"x = np.linspace(-0.1, 0.1, 100)\n",
|
|
||||||
"z = np.linspace(-0.1, 0.1, 100)\n",
|
|
||||||
"X, Z = np.meshgrid(x, z)\n",
|
|
||||||
"Bx, By, Bz = B(X, 0, Z)\n",
|
|
||||||
"\n",
|
|
||||||
"plot_field_slice(\n",
|
|
||||||
" X,\n",
|
|
||||||
" Z,\n",
|
|
||||||
" Bx,\n",
|
|
||||||
" Bz,\n",
|
|
||||||
" np.log(np.sqrt(Bx**2 + By**2 + Bz**2)),\n",
|
|
||||||
" \"magnetic_field_side.png\",\n",
|
|
||||||
")\n",
|
|
||||||
"\n",
|
|
||||||
"# # calculate the magnetic field at z = 1, over y = -1, 1 and x = -1, 1\n",
|
|
||||||
"x = np.linspace(-0.1, 0.1, 100)\n",
|
|
||||||
"y = np.linspace(-0.1, 0.1, 100)\n",
|
|
||||||
"X, Y = np.meshgrid(x, y)\n",
|
|
||||||
"Bx, By, Bz = B(X, Y, 0.01)\n",
|
|
||||||
"\n",
|
|
||||||
"plot_field_slice(\n",
|
|
||||||
" X, Y, Bx, By, np.sqrt(Bx**2 + By**2 + Bz**2), \"magnetic_field_bottom.png\"\n",
|
|
||||||
")\n",
|
|
||||||
"\n",
|
|
||||||
"# calculate the magnetic field in a 3d volume\n",
|
|
||||||
"x = np.linspace(-0.1, 0.1, 100)\n",
|
|
||||||
"y = np.linspace(-0.1, 0.1, 100)\n",
|
|
||||||
"z = np.linspace(-0.1, 0.1, 100)\n",
|
|
||||||
"X, Y, Z = np.meshgrid(x, y, z)\n",
|
|
||||||
"Bx, By, Bz = B(X, Y, Z)"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"# calculate the force on a wire of length l carrying current I at a point x,y,z with a direction vector d\n",
|
|
||||||
"def F(p, d, I, l):\n",
|
|
||||||
" return I * l * np.cross(d, B(p[0], p[1], p[2]))\n",
|
|
||||||
"\n",
|
|
||||||
"\n",
|
|
||||||
"def calculate_forces_on_wire_points(points):\n",
|
|
||||||
" Fx = []\n",
|
|
||||||
" Fy = []\n",
|
|
||||||
" Fz = []\n",
|
|
||||||
" # calculate the force on each point\n",
|
|
||||||
" for i in range(len(points)):\n",
|
|
||||||
" # calculate the direction vector\n",
|
|
||||||
" dx = points[i][0] - points[(i + 1) % len(points)][0]\n",
|
|
||||||
" dy = points[i][1] - points[(i + 1) % len(points)][1]\n",
|
|
||||||
" dz = points[i][2] - points[(i + 1) % len(points)][2]\n",
|
|
||||||
" d = np.array([dx, dy, dz])\n",
|
|
||||||
" # get the length of d\n",
|
|
||||||
" l = np.sqrt(dx**2 + dy**2 + dz**2)\n",
|
|
||||||
" if l > 0:\n",
|
|
||||||
" # normalise d\n",
|
|
||||||
" d = d / l\n",
|
|
||||||
" # calculate the force\n",
|
|
||||||
" fx, fy, fz = F(points[i], d, points[i][3], l)\n",
|
|
||||||
" Fx.append(fx)\n",
|
|
||||||
" Fy.append(fy)\n",
|
|
||||||
" Fz.append(fz)\n",
|
|
||||||
" else:\n",
|
|
||||||
" Fx.append(0)\n",
|
|
||||||
" Fy.append(0)\n",
|
|
||||||
" Fz.append(0)\n",
|
|
||||||
" return Fx, Fy, Fz\n",
|
|
||||||
"\n",
|
|
||||||
"\n",
|
|
||||||
"# locate the loop of wire directly below the magnet\n",
|
|
||||||
"x = 0.01\n",
|
|
||||||
"y = 0\n",
|
|
||||||
"z = -0.002\n",
|
|
||||||
"r1 = 0.001\n",
|
|
||||||
"r2 = 0.01\n",
|
|
||||||
"\n",
|
|
||||||
"\n",
|
|
||||||
"points = coil1.copy()\n",
|
|
||||||
"# scale the points from mm to m\n",
|
|
||||||
"# shift the coil to the correct position\n",
|
|
||||||
"for i in range(len(points)):\n",
|
|
||||||
" points[i][0] = points[i][0] / 1000 + x\n",
|
|
||||||
" points[i][1] = points[i][1] / 1000 + y\n",
|
|
||||||
" points[i][2] = points[i][2] / 1000 + z\n",
|
|
||||||
"\n",
|
|
||||||
"\n",
|
|
||||||
"# plot the points in 2D x,y\n",
|
|
||||||
"fig = plt.figure()\n",
|
|
||||||
"ax = fig.add_subplot(111)\n",
|
|
||||||
"ax.plot([p[0] for p in points], [p[1] for p in points])\n",
|
|
||||||
"ax.set_xlabel(\"$x$\")\n",
|
|
||||||
"ax.set_ylabel(\"$y$\")\n",
|
|
||||||
"ax.set_xlim(-0.01 + x, 0.01 + x)\n",
|
|
||||||
"ax.set_ylim(-0.01 + y, 0.01 + y)\n",
|
|
||||||
"ax.set_aspect(\"equal\")\n",
|
|
||||||
"fig.show()\n",
|
|
||||||
"\n",
|
|
||||||
"Fx, Fy, Fz = calculate_forces_on_wire_points(points)\n",
|
|
||||||
"\n",
|
|
||||||
"\n",
|
|
||||||
"# plot the wire along with arrows showing the force\n",
|
|
||||||
"fig = plt.figure()\n",
|
|
||||||
"ax = fig.add_subplot(111, projection=\"3d\")\n",
|
|
||||||
"ax.plot([p[0] for p in points], [p[1] for p in points], [p[2] for p in points])\n",
|
|
||||||
"\n",
|
|
||||||
"print(sum(Fx) / 9.8)\n",
|
|
||||||
"\n",
|
|
||||||
"# subsample the points and force vectors to make the plot clearer\n",
|
|
||||||
"points = points[::50]\n",
|
|
||||||
"Fx = Fx[::50]\n",
|
|
||||||
"Fy = Fy[::50]\n",
|
|
||||||
"Fz = Fz[::50]\n",
|
|
||||||
"\n",
|
|
||||||
"ax.quiver(\n",
|
|
||||||
" [p[0] for p in points],\n",
|
|
||||||
" [p[1] for p in points],\n",
|
|
||||||
" [p[2] for p in points],\n",
|
|
||||||
" np.sqrt(Fx),\n",
|
|
||||||
" 0,\n",
|
|
||||||
" 0,\n",
|
|
||||||
" length=0.01,\n",
|
|
||||||
" normalize=True,\n",
|
|
||||||
")\n",
|
|
||||||
"ax.set_xlabel(\"$x$\")\n",
|
|
||||||
"ax.set_ylabel(\"$y$\")\n",
|
|
||||||
"ax.set_zlabel(\"$z$\")\n",
|
|
||||||
"ax.set_xlim(x - 0.02, x + 0.02)\n",
|
|
||||||
"ax.set_ylim(y - 0.02, y + 0.02)\n",
|
|
||||||
"ax.set_zlim(z - 0.02, z + 0.02)\n",
|
|
||||||
"ax.set_aspect(\"equal\")\n",
|
|
||||||
"\n",
|
|
||||||
"# change the figure size\n",
|
|
||||||
"fig.set_size_inches(10, 10)\n",
|
|
||||||
"\n",
|
|
||||||
"fig.show()\n",
|
|
||||||
"# coil2 = -0.01311082557350831\n",
|
|
||||||
"# coil1 = 0.0088435517657609"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"def sweep_coil(coil, X):\n",
|
|
||||||
" Y = np.zeros(len(X))\n",
|
|
||||||
" Z = -0.01 * np.ones(len(X))\n",
|
|
||||||
"\n",
|
|
||||||
" # loop through the locations and calculate the forces, sum up the force in the X direction for each location\n",
|
|
||||||
" Fx = []\n",
|
|
||||||
" Fy = []\n",
|
|
||||||
" Fz = []\n",
|
|
||||||
" for p in trange(len(X)):\n",
|
|
||||||
" points = coil.copy()\n",
|
|
||||||
" # scale the points from mm to m\n",
|
|
||||||
" # shift the coil to the correct position\n",
|
|
||||||
" for i in range(len(points)):\n",
|
|
||||||
" points[i][0] = points[i][0] / 1000 + X[p]\n",
|
|
||||||
" points[i][1] = points[i][1] / 1000 + Y[p]\n",
|
|
||||||
" points[i][2] = points[i][2] / 1000 + Z[p]\n",
|
|
||||||
" Fx_, Fy_, Fz_ = calculate_forces_on_wire_points(points)\n",
|
|
||||||
" Fx.append(sum(Fx_))\n",
|
|
||||||
" Fy.append(sum(Fy_))\n",
|
|
||||||
" Fz.append(sum(Fz_))\n",
|
|
||||||
" return Fx, Fy, Fz\n",
|
|
||||||
"\n",
|
|
||||||
"\n",
|
|
||||||
"# sweep the coild from -3cm to 3cm in 0.01m steps\n",
|
|
||||||
"X = np.linspace(-0.03, 0.03, 100)\n",
|
|
||||||
"Fx_1_straight, Fy_1_straight, Fz_1_straight = sweep_coil(coil1, X)"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"# plot the force as a function of x\n",
|
|
||||||
"plt.plot(X, -np.array(Fx_1_straight) / 9.8, label=\"coil1\", color=\"red\")\n",
|
|
||||||
"# plot a dotted line along y = 0\n",
|
|
||||||
"plt.plot([X[0], X[-1]], [0, 0], \"--\")"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"# plot the force as a function of x\n",
|
|
||||||
"plt.plot(X, -np.array(Fz_1_straight) / 9.8, label=\"coil1\", color=\"red\")\n",
|
|
||||||
"# plot a dotted line along y = 0\n",
|
|
||||||
"plt.plot([X[0], X[-1]], [0, 0], \"--\")"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"# instead of sweeping horizontally, we'll sweep the coils around a circle\n",
|
|
||||||
"def sweep_coil_cirlc(coil, coil_center_radius, theta):\n",
|
|
||||||
" X = coil_center_radius * np.cos(np.deg2rad(theta))\n",
|
|
||||||
" Y = coil_center_radius * np.sin(np.deg2rad(theta))\n",
|
|
||||||
" Z = -0.01 * np.ones(100)\n",
|
|
||||||
"\n",
|
|
||||||
" # loop through the locations and calculate the forces, sum up the force in the X direction for each location\n",
|
|
||||||
" Torque = []\n",
|
|
||||||
" Fx = []\n",
|
|
||||||
" Fy = []\n",
|
|
||||||
" Fz = []\n",
|
|
||||||
" for p in trange(len(theta)):\n",
|
|
||||||
" angle = np.deg2rad(theta[p]) - np.pi / 2\n",
|
|
||||||
" x = X[p]\n",
|
|
||||||
" y = Y[p]\n",
|
|
||||||
" z = Z[p]\n",
|
|
||||||
"\n",
|
|
||||||
" points = coil.copy()\n",
|
|
||||||
" for i in range(len(points)):\n",
|
|
||||||
" px = points[i][0] / 1000\n",
|
|
||||||
" py = points[i][1] / 1000\n",
|
|
||||||
" pz = points[i][2] / 1000\n",
|
|
||||||
" # rotate the points so the coil is correctly oriented\n",
|
|
||||||
" points[i][0] = px * np.cos(angle) - py * np.sin(angle) + x\n",
|
|
||||||
" points[i][1] = (\n",
|
|
||||||
" px * np.sin(angle) + py * np.cos(angle) + y - coil_center_radius\n",
|
|
||||||
" )\n",
|
|
||||||
" points[i][2] = pz + z\n",
|
|
||||||
" # feel the force\n",
|
|
||||||
" Fx_, Fy_, Fz_ = calculate_forces_on_wire_points(points)\n",
|
|
||||||
" Fx.append(sum(Fx_))\n",
|
|
||||||
" Fy.append(sum(Fy_))\n",
|
|
||||||
" Fz.append(sum(Fz_))\n",
|
|
||||||
" # calculate the torque - which should be 90 degress to the angle\n",
|
|
||||||
" torque_angle = np.deg2rad(theta[p] - 90)\n",
|
|
||||||
" Torque.append(sum(Fx_) * np.cos(torque_angle) + sum(Fy_) * np.sin(torque_angle))\n",
|
|
||||||
"\n",
|
|
||||||
" return Fx, Fy, Fz, Torque\n",
|
|
||||||
"\n",
|
|
||||||
"\n",
|
|
||||||
"# sweep the coils from -45 to 45 degrees in 1 degree steps\n",
|
|
||||||
"theta = np.linspace(0, 180, 100)\n",
|
|
||||||
"Fx_1_curve, Fy_1_curve, Fz_1_curve, Torque_1 = sweep_coil_cirlc(\n",
|
|
||||||
" coil1, 19.5 / 1000, theta\n",
|
|
||||||
")"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"plt.plot(theta, -np.array(Torque_1) / 9.8, label=\"coil1\", color=\"red\")\n",
|
|
||||||
"# plot a dotted line along y = 0\n",
|
|
||||||
"plt.plot([theta[0], theta[-1]], [0, 0], \"--\")"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"# plot arrows for the Fx and Fy components\n",
|
|
||||||
"X = 19.5 * np.cos(np.deg2rad(theta)) / 1000\n",
|
|
||||||
"Y = 19.5 * np.sin(np.deg2rad(theta)) / 1000\n",
|
|
||||||
"plt.quiver(X[::5], Y[::5], Fx_1_curve[::5], Fy_1_curve[::5], color=\"red\")\n",
|
|
||||||
"# make the axis equal so the arrows are not stretched\n",
|
|
||||||
"plt.axis(\"equal\")"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"plt.plot(theta, -np.array(Fz_1_curve) / 9.8, label=\"coil1\", color=\"red\")\n",
|
|
||||||
"plt.plot(theta, np.array(Fz_2_curve) / 9.8, label=\"coil2\", color=\"blue\")\n",
|
|
||||||
"# plot a dotted line along y = 0\n",
|
|
||||||
"plt.plot([theta[0], theta[-1]], [0, 0], \"--\")"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"# instead of sweeping horizontally, we'll sweep the coils around a circle\n",
|
|
||||||
"def sweep_coil_cirlc(coil, coil_center_radius):\n",
|
|
||||||
" # sweep the coils from -45 to 45 degrees in 1 degree steps\n",
|
|
||||||
" theta = np.linspace(0, 180, 20)\n",
|
|
||||||
" X = coil_center_radius * np.cos(np.deg2rad(theta))\n",
|
|
||||||
" Y = coil_center_radius * np.sin(np.deg2rad(theta))\n",
|
|
||||||
" Z = 1 * np.ones(100)\n",
|
|
||||||
"\n",
|
|
||||||
" # loop through the locations and calculate the forces, sum up the force in the X direction for each location\n",
|
|
||||||
" Fx = []\n",
|
|
||||||
" Fy = []\n",
|
|
||||||
" Fz = []\n",
|
|
||||||
" for p in trange(len(theta)):\n",
|
|
||||||
" angle = np.deg2rad(theta[p]) - np.pi / 2\n",
|
|
||||||
" x = X[p]\n",
|
|
||||||
" y = Y[p]\n",
|
|
||||||
" z = Z[p]\n",
|
|
||||||
"\n",
|
|
||||||
" points = coil.copy()\n",
|
|
||||||
" for i in range(len(points)):\n",
|
|
||||||
" px = points[i][0] / 1000\n",
|
|
||||||
" py = points[i][1] / 1000\n",
|
|
||||||
" pz = points[i][2] / 1000\n",
|
|
||||||
" # rotate the points so the coil is correctly oriented\n",
|
|
||||||
" points[i][0] = px * np.cos(angle) - py * np.sin(angle) + x\n",
|
|
||||||
" points[i][1] = (\n",
|
|
||||||
" px * np.sin(angle) + py * np.cos(angle) + y - coil_center_radius\n",
|
|
||||||
" )\n",
|
|
||||||
" points[i][2] = pz + z\n",
|
|
||||||
" plt.plot([p[0] for p in points], [p[1] for p in points], linewidth=0.5)\n",
|
|
||||||
" # add the torque arrow to the plot\n",
|
|
||||||
" torque_angle = np.deg2rad(theta[p] - 90)\n",
|
|
||||||
" torque_x1 = x\n",
|
|
||||||
" torque_y1 = y\n",
|
|
||||||
" torque_x2 = x + 0.01 * np.cos(torque_angle)\n",
|
|
||||||
" torque_y2 = y + 0.01 * np.sin(torque_angle)\n",
|
|
||||||
" plt.arrow(\n",
|
|
||||||
" torque_x1,\n",
|
|
||||||
" torque_y1,\n",
|
|
||||||
" torque_x2 - torque_x1,\n",
|
|
||||||
" torque_y2 - torque_y1,\n",
|
|
||||||
" head_width=0.001,\n",
|
|
||||||
" head_length=0.002,\n",
|
|
||||||
" fc=\"k\",\n",
|
|
||||||
" ec=\"k\",\n",
|
|
||||||
" )\n",
|
|
||||||
"\n",
|
|
||||||
"\n",
|
|
||||||
"sweep_coil_cirlc(coil2, 19.5 / 1000)"
|
|
||||||
]
|
|
||||||
}
|
|
||||||
],
|
|
||||||
"metadata": {
|
|
||||||
"kernelspec": {
|
|
||||||
"display_name": "Python 3.10.7 ('venv': venv)",
|
|
||||||
"language": "python",
|
|
||||||
"name": "python3"
|
|
||||||
},
|
|
||||||
"language_info": {
|
|
||||||
"codemirror_mode": {
|
|
||||||
"name": "ipython",
|
|
||||||
"version": 3
|
|
||||||
},
|
|
||||||
"file_extension": ".py",
|
|
||||||
"mimetype": "text/x-python",
|
|
||||||
"name": "python",
|
|
||||||
"nbconvert_exporter": "python",
|
|
||||||
"pygments_lexer": "ipython3",
|
|
||||||
"version": "3.10.7"
|
|
||||||
},
|
|
||||||
"vscode": {
|
|
||||||
"interpreter": {
|
|
||||||
"hash": "1ce20143987840b9786ebb5907032c9c3a8efacbb887dbb0ebc4934f2ad26cb3"
|
|
||||||
}
|
|
||||||
}
|
|
||||||
},
|
|
||||||
"nbformat": 4,
|
|
||||||
"nbformat_minor": 2
|
|
||||||
}
|
|
||||||
389
simulations/magnetic_force_on_coils.ipynb
Normal file
389
simulations/magnetic_force_on_coils.ipynb
Normal file
|
|
@ -0,0 +1,389 @@
|
||||||
|
{
|
||||||
|
"cells": [
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": null,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [],
|
||||||
|
"source": [
|
||||||
|
"import numpy as np\n",
|
||||||
|
"import matplotlib.pyplot as plt\n",
|
||||||
|
"from mpl_toolkits.mplot3d import Axes3D\n",
|
||||||
|
"from biot_savart_v4_3 import parse_coil, plot_coil, slice_coil"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": null,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [],
|
||||||
|
"source": [
|
||||||
|
"# COIL = '6'\n",
|
||||||
|
"COIL = \"12\"\n",
|
||||||
|
"\n",
|
||||||
|
"# load up the simple spiral coil\n",
|
||||||
|
"coil1 = parse_coil(f\"coils/coil_{COIL}_spiral.csv\")\n",
|
||||||
|
"plot_coil(f\"coils/coil_{COIL}_spiral.csv\")\n",
|
||||||
|
"coil1 = slice_coil(coil1, 1)\n",
|
||||||
|
"coil1 = coil1.T\n",
|
||||||
|
"print(coil1.shape)\n",
|
||||||
|
"\n",
|
||||||
|
"coil2 = parse_coil(f\"coils/coil_{COIL}_custom.csv\")\n",
|
||||||
|
"plot_coil(f\"coils/coil_{COIL}_custom.csv\")\n",
|
||||||
|
"coil2 = slice_coil(coil2, 1)\n",
|
||||||
|
"coil2 = coil2.T\n",
|
||||||
|
"print(coil2.shape)"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": null,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [],
|
||||||
|
"source": [
|
||||||
|
"# show a simple dipole magnet\n",
|
||||||
|
"# magnetic field at a point x,y,z of a dipole magnet with moment m in the z direction\n",
|
||||||
|
"def dipole(x, y, z, m=0.185):\n",
|
||||||
|
" mu0 = 1e-7 / 4 * np.pi\n",
|
||||||
|
" r = np.sqrt(x ** 2 + y ** 2 + z ** 2)\n",
|
||||||
|
" return (\n",
|
||||||
|
" np.array(\n",
|
||||||
|
" [3 * x * z / r ** 5, 3 * y * z / r ** 5, (3 * z ** 2 / r ** 5 - 1 / r ** 3)]\n",
|
||||||
|
" )\n",
|
||||||
|
" * m\n",
|
||||||
|
" * mu0\n",
|
||||||
|
" )\n",
|
||||||
|
"\n",
|
||||||
|
"\n",
|
||||||
|
"def plot_field_slice(x, y, bx, by, mag, name=\"magnetic_field.png\", draw_magnet=True):\n",
|
||||||
|
" # plot the magnetic field\n",
|
||||||
|
" fig = plt.figure()\n",
|
||||||
|
" ax = fig.add_subplot(111)\n",
|
||||||
|
" ax.streamplot(\n",
|
||||||
|
" x,\n",
|
||||||
|
" y,\n",
|
||||||
|
" bx,\n",
|
||||||
|
" by,\n",
|
||||||
|
" linewidth=1,\n",
|
||||||
|
" cmap=plt.cm.inferno,\n",
|
||||||
|
" density=2,\n",
|
||||||
|
" arrowstyle=\"->\",\n",
|
||||||
|
" arrowsize=1.5,\n",
|
||||||
|
" )\n",
|
||||||
|
"\n",
|
||||||
|
" ax.set_xlabel(\"$x$\")\n",
|
||||||
|
" ax.set_ylabel(\"$y$\")\n",
|
||||||
|
" ax.set_xlim(-0.1, 0.1)\n",
|
||||||
|
" ax.set_ylim(-0.1, 0.1)\n",
|
||||||
|
" ax.set_aspect(\"equal\")\n",
|
||||||
|
"\n",
|
||||||
|
" # plot the magniture of the field as an image\n",
|
||||||
|
" im = ax.imshow(\n",
|
||||||
|
" mag, extent=[-0.1, 0.1, -0.1, 0.1], origin=\"lower\", cmap=plt.cm.inferno\n",
|
||||||
|
" )\n",
|
||||||
|
" if draw_magnet:\n",
|
||||||
|
" # draw the magnet\n",
|
||||||
|
" ax.add_patch(\n",
|
||||||
|
" plt.Rectangle((-0.005, -0.0015), 0.01, 0.003, fc=\"w\", ec=\"k\", lw=1)\n",
|
||||||
|
" )\n",
|
||||||
|
"\n",
|
||||||
|
" # make the figure bigger\n",
|
||||||
|
" fig.set_size_inches(10, 10)\n",
|
||||||
|
"\n",
|
||||||
|
" fig.show()\n",
|
||||||
|
" # save the figure\n",
|
||||||
|
" fig.savefig(name)\n",
|
||||||
|
"\n",
|
||||||
|
"\n",
|
||||||
|
"# # calculate the magnetic field at y = 0, over z = -1, 1 and x = -1, 1\n",
|
||||||
|
"x = np.linspace(-0.1, 0.1, 100)\n",
|
||||||
|
"z = np.linspace(-0.1, 0.1, 100)\n",
|
||||||
|
"X, Z = np.meshgrid(x, z)\n",
|
||||||
|
"Bx, By, Bz = dipole(X, 0, Z)\n",
|
||||||
|
"\n",
|
||||||
|
"print(Bx.shape, By.shape, Bz.shape)\n",
|
||||||
|
"\n",
|
||||||
|
"plot_field_slice(\n",
|
||||||
|
" X,\n",
|
||||||
|
" Z,\n",
|
||||||
|
" Bx,\n",
|
||||||
|
" Bz,\n",
|
||||||
|
" np.log(np.sqrt(Bx ** 2 + By ** 2 + Bz ** 2)),\n",
|
||||||
|
" \"dipole_field_side.png\",\n",
|
||||||
|
" False,\n",
|
||||||
|
")\n",
|
||||||
|
"\n",
|
||||||
|
"# # calculate the magnetic field at z = 1, over y = -1, 1 and x = -1, 1\n",
|
||||||
|
"x = np.linspace(-0.1, 0.1, 100)\n",
|
||||||
|
"y = np.linspace(-0.1, 0.1, 100)\n",
|
||||||
|
"X, Y = np.meshgrid(x, y)\n",
|
||||||
|
"Bx, By, Bz = dipole(X, Y, 0.01)\n",
|
||||||
|
"\n",
|
||||||
|
"plot_field_slice(\n",
|
||||||
|
" X,\n",
|
||||||
|
" Y,\n",
|
||||||
|
" Bx,\n",
|
||||||
|
" By,\n",
|
||||||
|
" np.log(np.sqrt(Bx ** 2 + By ** 2 + Bz ** 2)),\n",
|
||||||
|
" \"dipole_field_bottom.png\",\n",
|
||||||
|
" False,\n",
|
||||||
|
")"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "markdown",
|
||||||
|
"metadata": {},
|
||||||
|
"source": [
|
||||||
|
"# Simple Simulation of a dipole magnet"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": null,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [],
|
||||||
|
"source": [
|
||||||
|
"def B(x, y, z, m=0.185, l=0.003, d=0.01):\n",
|
||||||
|
" d = d * 0.75\n",
|
||||||
|
" # simulate multiple points in the cylinder and add them together to create a disk of field\n",
|
||||||
|
" bx = 0\n",
|
||||||
|
" by = 0\n",
|
||||||
|
" bz = 0\n",
|
||||||
|
" for degrees in range(0, 360, 30):\n",
|
||||||
|
" angle = np.deg2rad(degrees)\n",
|
||||||
|
" x_, y_, z_ = dipole(\n",
|
||||||
|
" x + 0.9 * d / 2 * np.cos(angle),\n",
|
||||||
|
" y + 0.9 * d / 2 * np.sin(angle),\n",
|
||||||
|
" z - 0.9 * l / 2,\n",
|
||||||
|
" m / (360 / 60),\n",
|
||||||
|
" )\n",
|
||||||
|
" bx += x_\n",
|
||||||
|
" by += y_\n",
|
||||||
|
" bz += z_\n",
|
||||||
|
" x_, y_, z_ = dipole(\n",
|
||||||
|
" x + 0.9 * d / 2 * np.cos(angle),\n",
|
||||||
|
" y + 0.9 * d / 2 * np.sin(angle),\n",
|
||||||
|
" z + 0.9 * l / 2,\n",
|
||||||
|
" m / (360 / 60),\n",
|
||||||
|
" )\n",
|
||||||
|
" bx += x_\n",
|
||||||
|
" by += y_\n",
|
||||||
|
" bz += z_\n",
|
||||||
|
" return np.array([bx, by, bz])\n",
|
||||||
|
"\n",
|
||||||
|
"\n",
|
||||||
|
"# # calculate the magnetic field at y = 0, over z = -1, 1 and x = -1, 1\n",
|
||||||
|
"x = np.linspace(-0.1, 0.1, 100)\n",
|
||||||
|
"z = np.linspace(-0.1, 0.1, 100)\n",
|
||||||
|
"X, Z = np.meshgrid(x, z)\n",
|
||||||
|
"Bx, By, Bz = B(X, 0, Z)\n",
|
||||||
|
"\n",
|
||||||
|
"print(Bx.shape, By.shape, Bz.shape)\n",
|
||||||
|
"\n",
|
||||||
|
"plot_field_slice(\n",
|
||||||
|
" X,\n",
|
||||||
|
" Z,\n",
|
||||||
|
" Bx,\n",
|
||||||
|
" Bz,\n",
|
||||||
|
" np.log(np.sqrt(Bx ** 2 + By ** 2 + Bz ** 2)),\n",
|
||||||
|
" \"magnetic_field_side.png\",\n",
|
||||||
|
")\n",
|
||||||
|
"\n",
|
||||||
|
"# # calculate the magnetic field at z = 1, over y = -1, 1 and x = -1, 1\n",
|
||||||
|
"x = np.linspace(-0.1, 0.1, 100)\n",
|
||||||
|
"y = np.linspace(-0.1, 0.1, 100)\n",
|
||||||
|
"X, Y = np.meshgrid(x, y)\n",
|
||||||
|
"Bx, By, Bz = B(X, Y, 0.01)\n",
|
||||||
|
"\n",
|
||||||
|
"plot_field_slice(\n",
|
||||||
|
" X,\n",
|
||||||
|
" Y,\n",
|
||||||
|
" Bx,\n",
|
||||||
|
" By,\n",
|
||||||
|
" np.log(np.sqrt(Bx ** 2 + By ** 2 + Bz ** 2)),\n",
|
||||||
|
" \"magnetic_field_bottom.png\",\n",
|
||||||
|
")\n",
|
||||||
|
"\n",
|
||||||
|
"# calculate the magnetic field in a 3d volume\n",
|
||||||
|
"x = np.linspace(-0.1, 0.1, 100)\n",
|
||||||
|
"y = np.linspace(-0.1, 0.1, 100)\n",
|
||||||
|
"z = np.linspace(-0.1, 0.1, 100)\n",
|
||||||
|
"X, Y, Z = np.meshgrid(x, y, z)\n",
|
||||||
|
"Bx, By, Bz = B(X, Y, Z)"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": null,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [],
|
||||||
|
"source": [
|
||||||
|
"# do a 3d quivwer plot of the magnetic field\n",
|
||||||
|
"fig = plt.figure()\n",
|
||||||
|
"ax = fig.add_subplot(111, projection=\"3d\")\n",
|
||||||
|
"# down sample the results\n",
|
||||||
|
"ax.quiver(\n",
|
||||||
|
" X[::10, ::10, ::10],\n",
|
||||||
|
" Y[::10, ::10, ::10],\n",
|
||||||
|
" Z[::10, ::10, ::10],\n",
|
||||||
|
" Bx[::10, ::10, ::10],\n",
|
||||||
|
" By[::10, ::10, ::10],\n",
|
||||||
|
" Bz[::10, ::10, ::10],\n",
|
||||||
|
" length=0.01,\n",
|
||||||
|
" normalize=True,\n",
|
||||||
|
")\n",
|
||||||
|
"ax.set_xlabel(\"$x$\")\n",
|
||||||
|
"ax.set_ylabel(\"$y$\")\n",
|
||||||
|
"ax.set_zlabel(\"$z$\")\n",
|
||||||
|
"ax.set_xlim(-0.1, 0.1)\n",
|
||||||
|
"ax.set_ylim(-0.1, 0.1)\n",
|
||||||
|
"ax.set_zlim(-0.1, 0.1)\n",
|
||||||
|
"ax.set_aspect(\"equal\")\n",
|
||||||
|
"# make the plot larger\n",
|
||||||
|
"fig.set_size_inches(10, 10)\n",
|
||||||
|
"fig.show()"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": null,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [],
|
||||||
|
"source": [
|
||||||
|
"# calculate the force on a wire of length l carrying current I at a point x,y,z with a direction vector d\n",
|
||||||
|
"def F(p, d, I, l):\n",
|
||||||
|
" return I * l * np.cross(d, B(p[0], p[1], p[2]))\n",
|
||||||
|
"\n",
|
||||||
|
"\n",
|
||||||
|
"def calculate_forces_on_wire_points(points):\n",
|
||||||
|
" Fx = []\n",
|
||||||
|
" Fy = []\n",
|
||||||
|
" Fz = []\n",
|
||||||
|
" # calculate the force on each point\n",
|
||||||
|
" for i in range(len(points) - 1):\n",
|
||||||
|
" # calculate the direction vector\n",
|
||||||
|
" dx = points[i][0] - points[(i + 1)][0]\n",
|
||||||
|
" dy = points[i][1] - points[(i + 1)][1]\n",
|
||||||
|
" dz = points[i][2] - points[(i + 1)][2]\n",
|
||||||
|
" d = np.array([dx, dy, dz])\n",
|
||||||
|
" # get the length of d\n",
|
||||||
|
" l = np.sqrt(dx ** 2 + dy ** 2 + dz ** 2)\n",
|
||||||
|
" if l > 0:\n",
|
||||||
|
" # normalise d\n",
|
||||||
|
" d = d / l\n",
|
||||||
|
" # calculate the force\n",
|
||||||
|
" fx, fy, fz = F(points[i], d, points[i][3], l)\n",
|
||||||
|
" Fx.append(fx)\n",
|
||||||
|
" Fy.append(fy)\n",
|
||||||
|
" Fz.append(fz)\n",
|
||||||
|
" else:\n",
|
||||||
|
" Fx.append(0)\n",
|
||||||
|
" Fy.append(0)\n",
|
||||||
|
" Fz.append(0)\n",
|
||||||
|
" return Fx, Fy, Fz"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": null,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [],
|
||||||
|
"source": [
|
||||||
|
"# instead of sweeping horizontally, we'll sweep the coils around a circle\n",
|
||||||
|
"def sweep_coil_circle(coil, coil_center_radius, theta):\n",
|
||||||
|
" X = coil_center_radius * np.cos(np.deg2rad(theta))\n",
|
||||||
|
" Y = coil_center_radius * np.sin(np.deg2rad(theta)) - coil_center_radius\n",
|
||||||
|
" Z = -0.01\n",
|
||||||
|
"\n",
|
||||||
|
" # loop through the locations and calculate the forces, sum up the force in the X direction for each location\n",
|
||||||
|
" Fx = []\n",
|
||||||
|
" Fy = []\n",
|
||||||
|
" Fz = []\n",
|
||||||
|
" for p in range(len(theta)):\n",
|
||||||
|
" x = X[p]\n",
|
||||||
|
" y = Y[p]\n",
|
||||||
|
" z = Z\n",
|
||||||
|
"\n",
|
||||||
|
" points = coil.copy()\n",
|
||||||
|
" for i in range(len(points)):\n",
|
||||||
|
" points[i][0] = points[i][0] / 1000 + x\n",
|
||||||
|
" points[i][1] = points[i][1] / 1000 + y\n",
|
||||||
|
" points[i][2] = points[i][2] / 1000 + z\n",
|
||||||
|
" # feel the force\n",
|
||||||
|
" Fx_, Fy_, Fz_ = calculate_forces_on_wire_points(points)\n",
|
||||||
|
" Fx.append(sum(Fx_))\n",
|
||||||
|
" Fy.append(sum(Fy_))\n",
|
||||||
|
" Fz.append(sum(Fz_))\n",
|
||||||
|
" return Fx, Fy, Fz\n",
|
||||||
|
"\n",
|
||||||
|
"\n",
|
||||||
|
"# sweep the coils from -45 to 45 degrees in 1 degree steps\n",
|
||||||
|
"theta = np.linspace(25, 175, 100)\n",
|
||||||
|
"\n",
|
||||||
|
"# do a quick sanity check\n",
|
||||||
|
"# coil center in meters\n",
|
||||||
|
"coil_center_radius = 20.5 / 1000\n",
|
||||||
|
"X = coil_center_radius * np.cos(np.deg2rad(theta))\n",
|
||||||
|
"Y = coil_center_radius * np.sin(np.deg2rad(theta)) - coil_center_radius\n",
|
||||||
|
"\n",
|
||||||
|
"# plot X and Y along with the coil\n",
|
||||||
|
"plt.plot(X, Y, color=\"red\")\n",
|
||||||
|
"points = coil2.copy()\n",
|
||||||
|
"plt.plot(\n",
|
||||||
|
" # points in meters (coil is in cm)\n",
|
||||||
|
" [x / 100 for x, y, z, _ in points],\n",
|
||||||
|
" [y / 100 for x, y, z, _ in points],\n",
|
||||||
|
" label=\"coil\",\n",
|
||||||
|
" color=\"blue\",\n",
|
||||||
|
")\n",
|
||||||
|
"# make the aspect ratio equal\n",
|
||||||
|
"plt.gca().set_aspect(\"equal\")\n",
|
||||||
|
"\n",
|
||||||
|
"Fx_1_curve, Fy_1_curve, Fz_1_curve = sweep_coil_circle(coil1, 20.5 / 1000, theta)\n",
|
||||||
|
"Fx_2_curve, Fy_2_curve, Fz_2_curve = sweep_coil_circle(coil2, 19.5 / 1000, theta)"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": null,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [],
|
||||||
|
"source": [
|
||||||
|
"plt.plot(theta, np.array(Fx_1_curve), label=\"coil1\", color=\"red\")\n",
|
||||||
|
"plt.plot(theta, -np.array(Fx_2_curve), label=\"coil2\", color=\"blue\")\n",
|
||||||
|
"# plot a dotted line along y = 0\n",
|
||||||
|
"plt.gcf().set_size_inches(16, 9)\n",
|
||||||
|
"plt.plot([theta[0], theta[-1]], [0, 0], \"--\")\n",
|
||||||
|
"print(np.max(np.abs(Fx_1_curve)))\n",
|
||||||
|
"# save the plot\n",
|
||||||
|
"plt.savefig(f\"torque_{COIL}_coil_sweep_circle_1mm_spiral.png\")"
|
||||||
|
]
|
||||||
|
}
|
||||||
|
],
|
||||||
|
"metadata": {
|
||||||
|
"kernelspec": {
|
||||||
|
"display_name": "Python 3.10.7 ('venv': venv)",
|
||||||
|
"language": "python",
|
||||||
|
"name": "python3"
|
||||||
|
},
|
||||||
|
"language_info": {
|
||||||
|
"codemirror_mode": {
|
||||||
|
"name": "ipython",
|
||||||
|
"version": 3
|
||||||
|
},
|
||||||
|
"file_extension": ".py",
|
||||||
|
"mimetype": "text/x-python",
|
||||||
|
"name": "python",
|
||||||
|
"nbconvert_exporter": "python",
|
||||||
|
"pygments_lexer": "ipython3",
|
||||||
|
"version": "3.10.7"
|
||||||
|
},
|
||||||
|
"vscode": {
|
||||||
|
"interpreter": {
|
||||||
|
"hash": "fc384f9db26c31784edfba3761ba3d2c7b2f9b8a63e03a9eb0778fc35334efe1"
|
||||||
|
}
|
||||||
|
}
|
||||||
|
},
|
||||||
|
"nbformat": 4,
|
||||||
|
"nbformat_minor": 2
|
||||||
|
}
|
||||||
|
|
@ -1,162 +0,0 @@
|
||||||
{
|
|
||||||
"cells": [
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"import numpy as np\n",
|
|
||||||
"import biot_savart_v4_3 as bs\n",
|
|
||||||
"import matplotlib.pyplot as plt"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"bs.plot_coil(\"coil-4-top-1.csv\")"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"bs.plot_coil(\"coil-4-bottom-1.csv\")"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"for i in range(0, 12):\n",
|
|
||||||
" bs.write_target_volume(\n",
|
|
||||||
" f\"coil-4-top-{i}.csv\",\n",
|
|
||||||
" f\"targetvol-top-{i}.vol\",\n",
|
|
||||||
" (6, 6, 1),\n",
|
|
||||||
" (-3, -3, -0.5),\n",
|
|
||||||
" 0.1,\n",
|
|
||||||
" 0.1,\n",
|
|
||||||
" )\n",
|
|
||||||
" bs.write_target_volume(\n",
|
|
||||||
" f\"coil-4-bottom-{i}.csv\",\n",
|
|
||||||
" f\"targetvol-bottom-{i}.vol\",\n",
|
|
||||||
" (6, 6, 1),\n",
|
|
||||||
" (-3, -3, -0.5),\n",
|
|
||||||
" 0.1,\n",
|
|
||||||
" 0.1,\n",
|
|
||||||
" )"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"# read in the results\n",
|
|
||||||
"top = [bs.read_target_volume(f\"targetvol-top-{i}.vol\") for i in range(0, 12)]\n",
|
|
||||||
"bottom = [bs.read_target_volume(f\"targetvol-bottom-{i}.vol\") for i in range(0, 12)]\n",
|
|
||||||
"\n",
|
|
||||||
"# add all the results together into one volume\n",
|
|
||||||
"volume = top[0] + bottom[0]\n",
|
|
||||||
"for i in range(1, 12):\n",
|
|
||||||
" volume += top[i]\n",
|
|
||||||
" volume += bottom[i]\n",
|
|
||||||
"\n",
|
|
||||||
"print(volume.shape)\n",
|
|
||||||
"\n",
|
|
||||||
"# reads the volume we created\n",
|
|
||||||
"bs.plot_fields(volume, (60, 60, 10), (0, 0, 0), 1, which_plane=\"z\", level=6)"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"import matplotlib as plt\n",
|
|
||||||
"\n",
|
|
||||||
"# get a slice of the volume at z=6\n",
|
|
||||||
"slice = [volume[50, :, :, i].T for i in range(3)]\n",
|
|
||||||
"x_contour = slice[0]\n",
|
|
||||||
"y_contour = slice[1]\n",
|
|
||||||
"z_contour = slice[2]\n",
|
|
||||||
"\n",
|
|
||||||
"print(x_contour.shape)\n",
|
|
||||||
"print(y_contour.shape)\n",
|
|
||||||
"print(z_contour.shape)\n",
|
|
||||||
"\n",
|
|
||||||
"x = np.arange(0, 61, 1)\n",
|
|
||||||
"y = np.arange(0, 11, 1)\n",
|
|
||||||
"xG, yG = np.meshgrid(x, y)\n",
|
|
||||||
"# do a quiver plot of the x and the z components - we won't worry about y\n",
|
|
||||||
"plt.pyplot.quiver(xG, yG, y_contour, z_contour)"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": null,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"fig = plt.pyplot.figure()\n",
|
|
||||||
"ax = fig.add_subplot(111)\n",
|
|
||||||
"\n",
|
|
||||||
"# Plot the streamlines with an appropriate colormap and arrow style\n",
|
|
||||||
"# color = 2 * np.log(np.hypot(x_contour, z_contour))\n",
|
|
||||||
"ax.streamplot(\n",
|
|
||||||
" x,\n",
|
|
||||||
" y,\n",
|
|
||||||
" y_contour,\n",
|
|
||||||
" z_contour,\n",
|
|
||||||
" linewidth=1,\n",
|
|
||||||
" cmap=plt.cm.inferno,\n",
|
|
||||||
" density=2,\n",
|
|
||||||
" arrowstyle=\"->\",\n",
|
|
||||||
" arrowsize=1.5,\n",
|
|
||||||
")\n",
|
|
||||||
"\n",
|
|
||||||
"ax.set_xlabel(\"$x$\")\n",
|
|
||||||
"ax.set_ylabel(\"$y$\")\n",
|
|
||||||
"# ax.set_xlim(-2,2)\n",
|
|
||||||
"# ax.set_ylim(-2,2)\n",
|
|
||||||
"ax.set_aspect(\"equal\")\n",
|
|
||||||
"fig.show()"
|
|
||||||
]
|
|
||||||
}
|
|
||||||
],
|
|
||||||
"metadata": {
|
|
||||||
"kernelspec": {
|
|
||||||
"display_name": "Python 3.10.7 ('venv': venv)",
|
|
||||||
"language": "python",
|
|
||||||
"name": "python3"
|
|
||||||
},
|
|
||||||
"language_info": {
|
|
||||||
"codemirror_mode": {
|
|
||||||
"name": "ipython",
|
|
||||||
"version": 3
|
|
||||||
},
|
|
||||||
"file_extension": ".py",
|
|
||||||
"mimetype": "text/x-python",
|
|
||||||
"name": "python",
|
|
||||||
"nbconvert_exporter": "python",
|
|
||||||
"pygments_lexer": "ipython3",
|
|
||||||
"version": "3.10.7"
|
|
||||||
},
|
|
||||||
"vscode": {
|
|
||||||
"interpreter": {
|
|
||||||
"hash": "1ce20143987840b9786ebb5907032c9c3a8efacbb887dbb0ebc4934f2ad26cb3"
|
|
||||||
}
|
|
||||||
}
|
|
||||||
},
|
|
||||||
"nbformat": 4,
|
|
||||||
"nbformat_minor": 2
|
|
||||||
}
|
|
||||||
Loading…
Reference in a new issue