mirror of
https://github.com/atomic14/kicad-coil-plugins.git
synced 2024-10-18 09:06:57 +00:00
410 lines
12 KiB
Text
410 lines
12 KiB
Text
{
|
|
"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, 0.01)\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, 0.01)\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",
|
|
" l = l * 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, 10):\n",
|
|
" angle = np.deg2rad(degrees)\n",
|
|
" x_, y_, z_ = dipole(\n",
|
|
" x + d / 2 * np.cos(angle),\n",
|
|
" y + d / 2 * np.sin(angle),\n",
|
|
" z - l / 2,\n",
|
|
" m / (360 / 10),\n",
|
|
" )\n",
|
|
" bx += x_\n",
|
|
" by += y_\n",
|
|
" bz += z_\n",
|
|
" x_, y_, z_ = dipole(\n",
|
|
" x + d / 2 * np.cos(angle),\n",
|
|
" y + d / 2 * np.sin(angle),\n",
|
|
" z + l / 2,\n",
|
|
" m / (360 / 10),\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": [
|
|
"import multiprocess as mp\n",
|
|
"\n",
|
|
"\n",
|
|
"def calculate_forces_on_coil(coil, x, y, z):\n",
|
|
" points = coil.copy()\n",
|
|
" for i in range(len(points)):\n",
|
|
" points[i][0] = points[i][0] / 100 + x\n",
|
|
" points[i][1] = points[i][1] / 100 + y\n",
|
|
" points[i][2] = points[i][2] / 100 + z\n",
|
|
" # feel the force\n",
|
|
" Fx_, Fy_, Fz_ = calculate_forces_on_wire_points(points)\n",
|
|
" return sum(Fx_), sum(Fy_), sum(Fz_)\n",
|
|
"\n",
|
|
"\n",
|
|
"# 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",
|
|
"\n",
|
|
" # run this in parallel\n",
|
|
" with mp.Pool(mp.cpu_count()) as pool:\n",
|
|
" params = [(coil, X[p], Y[p], Z) for p in range(len(theta))]\n",
|
|
" results = [\n",
|
|
" pool.apply_async(calculate_forces_on_coil, params) for params in params\n",
|
|
" ]\n",
|
|
"\n",
|
|
" for result in results:\n",
|
|
" fx, fy, fz = result.get()\n",
|
|
" Fx.append(fx)\n",
|
|
" Fy.append(fy)\n",
|
|
" Fz.append(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=\"coil2\",\n",
|
|
" color=\"blue\",\n",
|
|
")\n",
|
|
"points = coil1.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=\"coil1\",\n",
|
|
" color=\"red\",\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, 20.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
|
|
}
|