pcb-stator-coil-generator/simulations/magnet_field.ipynb

{
 "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
}