From fcc9c6727c4c3c2efc9caa61da70debc02897ed3 Mon Sep 17 00:00:00 2001 From: Chris Greening Date: Tue, 22 Nov 2022 11:58:21 +0000 Subject: [PATCH] Cleaned up simulation code --- coil_generator-12.ipynb | 2 +- coil_generator-6.ipynb | 2 +- requirements.txt | 5 - simulations/coil_6_spiral_coil.svg | 1 - ...ields.ipynb => coil_magnetic_fields.ipynb} | 26 +- simulations/helpers.py | 111 ---- simulations/magnet_field.ipynb | 607 ------------------ simulations/magnet_field_single_coil.ipynb | 525 --------------- simulations/magnetic_force_on_coils.ipynb | 389 +++++++++++ simulations/simulate_coil.ipynb | 162 ----- 10 files changed, 411 insertions(+), 1419 deletions(-) delete mode 100644 simulations/coil_6_spiral_coil.svg rename simulations/{compare_coil_fields.ipynb => coil_magnetic_fields.ipynb} (92%) delete mode 100644 simulations/helpers.py delete mode 100644 simulations/magnet_field.ipynb delete mode 100644 simulations/magnet_field_single_coil.ipynb create mode 100644 simulations/magnetic_force_on_coils.ipynb delete mode 100644 simulations/simulate_coil.ipynb diff --git a/coil_generator-12.ipynb b/coil_generator-12.ipynb index 8eb098c..39e726f 100644 --- a/coil_generator-12.ipynb +++ b/coil_generator-12.ipynb @@ -714,7 +714,7 @@ }, "vscode": { "interpreter": { - "hash": "1ce20143987840b9786ebb5907032c9c3a8efacbb887dbb0ebc4934f2ad26cb3" + "hash": "fc384f9db26c31784edfba3761ba3d2c7b2f9b8a63e03a9eb0778fc35334efe1" } } }, diff --git a/coil_generator-6.ipynb b/coil_generator-6.ipynb index daf459a..0f99dba 100644 --- a/coil_generator-6.ipynb +++ b/coil_generator-6.ipynb @@ -693,7 +693,7 @@ }, "vscode": { "interpreter": { - "hash": "1ce20143987840b9786ebb5907032c9c3a8efacbb887dbb0ebc4934f2ad26cb3" + "hash": "fc384f9db26c31784edfba3761ba3d2c7b2f9b8a63e03a9eb0778fc35334efe1" } } }, diff --git a/requirements.txt b/requirements.txt index e693d95..d76c92f 100644 --- a/requirements.txt +++ b/requirements.txt @@ -4,9 +4,4 @@ numpy pandas scikit-spatial pre-commit -tqdm plotly -magpylib -pyvista -mayavi -PySide \ No newline at end of file diff --git a/simulations/coil_6_spiral_coil.svg b/simulations/coil_6_spiral_coil.svg deleted file mode 100644 index 5df132c..0000000 --- a/simulations/coil_6_spiral_coil.svg +++ /dev/null @@ -1 +0,0 @@ - \ No newline at end of file diff --git a/simulations/compare_coil_fields.ipynb b/simulations/coil_magnetic_fields.ipynb similarity index 92% rename from simulations/compare_coil_fields.ipynb rename to simulations/coil_magnetic_fields.ipynb index 8a095df..3acd0e6 100644 --- a/simulations/compare_coil_fields.ipynb +++ b/simulations/coil_magnetic_fields.ipynb @@ -68,7 +68,7 @@ "# read in the results\n", "total_volume = bs.read_target_volume(COIL_NAME + \"volume.vol\")\n", "\n", - "total_volume = total_volume\n", + "total_volume = -total_volume\n", "\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)" @@ -141,7 +141,7 @@ "xG, yG = np.meshgrid(x, y)\n", "\n", "# calclulate the magnitude of the field\n", - "mag = np.sqrt(x_contour**2 + y_contour**2 + z_contour**2)\n", + "mag = np.sqrt(x_contour ** 2 + y_contour ** 2 + z_contour ** 2)\n", "\n", "# color = 2 * np.log(np.hypot(x_contour, z_contour))\n", "color = np.log(mag)\n", @@ -204,16 +204,30 @@ "# make the plot larger\n", "plt.rcParams[\"figure.figsize\"] = (10, 10)\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.contourf(\n", - " x_contour,\n", + " y_contour,\n", " vmin=vmin,\n", " vmax=vmax,\n", " cmap=cm.magma,\n", " levels=50,\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", "print(vmin, vmax)\n", "\n", @@ -255,7 +269,7 @@ "y_contour = slice[1]\n", "z_contour = slice[2]\n", "\n", - "mag = np.sqrt(x_contour**2 + y_contour**2 + z_contour**2)\n", + "mag = np.sqrt(x_contour ** 2 + y_contour ** 2 + z_contour ** 2)\n", "\n", "print(np.max(mag), np.min(mag))\n", "\n", @@ -350,7 +364,7 @@ }, "vscode": { "interpreter": { - "hash": "1ce20143987840b9786ebb5907032c9c3a8efacbb887dbb0ebc4934f2ad26cb3" + "hash": "fc384f9db26c31784edfba3761ba3d2c7b2f9b8a63e03a9eb0778fc35334efe1" } } }, diff --git a/simulations/helpers.py b/simulations/helpers.py deleted file mode 100644 index 234576b..0000000 --- a/simulations/helpers.py +++ /dev/null @@ -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) diff --git a/simulations/magnet_field.ipynb b/simulations/magnet_field.ipynb deleted file mode 100644 index c10f20e..0000000 --- a/simulations/magnet_field.ipynb +++ /dev/null @@ -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 -} diff --git a/simulations/magnet_field_single_coil.ipynb b/simulations/magnet_field_single_coil.ipynb deleted file mode 100644 index a41c36c..0000000 --- a/simulations/magnet_field_single_coil.ipynb +++ /dev/null @@ -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 -} diff --git a/simulations/magnetic_force_on_coils.ipynb b/simulations/magnetic_force_on_coils.ipynb new file mode 100644 index 0000000..d484f97 --- /dev/null +++ b/simulations/magnetic_force_on_coils.ipynb @@ -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 +} diff --git a/simulations/simulate_coil.ipynb b/simulations/simulate_coil.ipynb deleted file mode 100644 index 31ffbe3..0000000 --- a/simulations/simulate_coil.ipynb +++ /dev/null @@ -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 -}