diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000..61bed8f Binary files /dev/null and b/.DS_Store differ diff --git a/.gitignore b/.gitignore index e8f28be..676cc8e 100644 --- a/.gitignore +++ b/.gitignore @@ -160,4 +160,6 @@ cython_debug/ #.idea/ *.json -*.png \ No newline at end of file +*.png +*.csv +*.vol \ No newline at end of file diff --git a/README.md b/README.md index e370a19..43568fb 100644 --- a/README.md +++ b/README.md @@ -25,6 +25,10 @@ ln -s ${PWD}/coil_plugin.py ~/Documents/KiCad/6.0/scripting/plugins/coil_plugin. The project is still very experimental - so there are no guarantees that these will work or do anything useful... -* 6 circular coils - https://www.pcbway.com/project/shareproject/Proof_of_concept_6_coil_spiral_PCB_stator_23ae7370.html -* 6 wedge coils - https://www.pcbway.com/project/shareproject/Proof_of_concept_6_coil_wedge_PCB_stator_c231a379.html -* 12 wedge coils - https://www.pcbway.com/project/shareproject/Proof_of_concept_12_coil_wedge_PCB_stator_c54d9374.html +- 6 circular coils - https://www.pcbway.com/project/shareproject/Proof_of_concept_6_coil_spiral_PCB_stator_23ae7370.html +- 6 wedge coils - https://www.pcbway.com/project/shareproject/Proof_of_concept_6_coil_wedge_PCB_stator_c231a379.html +- 12 wedge coils - https://www.pcbway.com/project/shareproject/Proof_of_concept_12_coil_wedge_PCB_stator_c54d9374.html + +## Building python to take advantage of the M1 Mac Multiprocessors + +https://www.reddit.com/r/Python/comments/qog8x3/if_you_are_using_apples_m1_macs_compiling_numpy/ diff --git a/coil_generator-12.ipynb b/coil_generator-12.ipynb index 20cc99e..8eb098c 100644 --- a/coil_generator-12.ipynb +++ b/coil_generator-12.ipynb @@ -106,7 +106,10 @@ "SCREW_HOLE_RADIUS = STATOR_RADIUS\n", "\n", "# Coil params\n", + "# for custom shape\n", "TURNS = 16\n", + "# for spiral\n", + "# TURNS = 17\n", "COIL_CENTER_RADIUS = 19.95\n", "COIL_VIA_RADIUS = 20.95" ] @@ -120,7 +123,12 @@ "# where to put the input pads\n", "INPUT_PAD_RADIUS = STATOR_RADIUS - (PAD_WIDTH / 2 + VIA_DIAM + TRACK_SPACING)\n", "\n", - "USE_SPIRAL = False\n", + "USE_SPIRAL = True\n", + "\n", + "if USE_SPIRAL:\n", + " TURNS = 18\n", + " COIL_VIA_RADIUS = 20.5\n", + " COIL_CENTER_RADIUS = 20.5\n", "\n", "LAYERS = 4" ] @@ -238,18 +246,63 @@ "metadata": {}, "outputs": [], "source": [ - "template_f = []\n", - "for i in range(len(template)):\n", - " template_f.append(template[len(template) - i - len(template) // 2])\n", - "template_f = flip_x(template_f)\n", - "points_f = chaikin(\n", - " optimize_points(flip_x(get_points(template_f, TURNS, TRACK_SPACING + TRACK_WIDTH))),\n", - " 2,\n", - ")\n", - "points_b = chaikin(\n", - " optimize_points(get_points(template, TURNS, TRACK_SPACING + TRACK_WIDTH)), 2\n", - ")\n", + "def get_spiral(turns, start_radius, thickness, layer=Layer.FRONT):\n", + " points = []\n", + " # create a starting point in the center\n", + " for angle in np.arange(0, turns * 360, 1):\n", + " radius = start_radius + thickness * angle / 360\n", + " if layer == Layer.BACK:\n", + " x = radius * np.cos(np.deg2rad(angle + 180))\n", + " y = radius * np.sin(np.deg2rad(angle + 180))\n", + " points.append((x, -y))\n", + " else:\n", + " x = radius * np.cos(np.deg2rad(angle))\n", + " y = radius * np.sin(np.deg2rad(angle))\n", + " points.append((x, y))\n", + " return points" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "if not USE_SPIRAL:\n", + " print(\"Not using spiral\")\n", + " template_f = []\n", + " for i in range(len(template)):\n", + " template_f.append(template[len(template) - i - len(template) // 2])\n", + " template_f = flip_x(template_f)\n", + " points_f = chaikin(\n", + " optimize_points(\n", + " flip_x(get_points(template_f, TURNS, TRACK_SPACING + TRACK_WIDTH))\n", + " ),\n", + " 2,\n", + " )\n", + " points_b = chaikin(\n", + " optimize_points(get_points(template, TURNS, TRACK_SPACING + TRACK_WIDTH)), 2\n", + " )\n", + "else:\n", + " print(\"Using spiral\")\n", + " points_f = get_spiral(\n", + " TURNS, VIA_DIAM / 2 + TRACK_SPACING, TRACK_SPACING + TRACK_WIDTH, Layer.FRONT\n", + " )\n", + " points_b = get_spiral(\n", + " TURNS, VIA_DIAM / 2 + TRACK_SPACING, TRACK_SPACING + TRACK_WIDTH, Layer.BACK\n", + " )\n", "\n", + " points_f = [(0, 0)] + points_f\n", + " points_b = [(0, 0)] + points_b\n", + " print(\"Track points\", len(points_f), len(points_b))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ "points_f = [(COIL_VIA_RADIUS - COIL_CENTER_RADIUS, 0)] + points_f\n", "points_b = [(COIL_VIA_RADIUS - COIL_CENTER_RADIUS, 0)] + points_b\n", "\n", @@ -262,6 +315,47 @@ "print(\"Track points\", len(points_f), len(points_b))" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# write the coil out in a format that can be simulated\n", + "# rotate the points by 90 degrees so that the x axis is horizontal\n", + "pf = rotate(points_f, 90)\n", + "pb = rotate(points_b, 90)\n", + "fname = \"simulations/coils/coil_12_custom\"\n", + "if USE_SPIRAL:\n", + " fname = \"simulations/coils/coil_12_spiral\"\n", + "\n", + "# write the coil out in a format that can be simulated\n", + "p_f = rotate(points_f, 90)\n", + "p_b = rotate(points_b, 90)\n", + "\n", + "with open(fname + \".csv\", \"w\") as f:\n", + " for point in p_f[::-1]:\n", + " f.write(f\"{point[0]},{point[1]},0,0.5\\n\")\n", + "\n", + "# two layer board\n", + "with open(fname + \"-2-layer.csv\", \"wt\") as f:\n", + " for point in p_f[::-1]:\n", + " f.write(f\"{point[0]},{point[1]},0,0.5\\n\")\n", + " for point in p_b:\n", + " f.write(f\"{point[0]},{point[1]},0-0.062,0.5\\n\")\n", + "\n", + "# all four layer board\n", + "with open(fname + \"-4-layer.csv\", \"wt\") as f:\n", + " for point in p_f[::-1]:\n", + " f.write(f\"{point[0]},{point[1]},0,0.5\\n\")\n", + " for point in p_b:\n", + " f.write(f\"{point[0]},{point[1]},0-0.011,0.5\\n\")\n", + " for point in p_f[::-1]:\n", + " f.write(f\"{point[0]},{point[1]},0-(0.011+0.04),0.5\\n\")\n", + " for point in p_b:\n", + " f.write(f\"{point[0]},{point[1]},0-(0.011+0.011+0.04),0.5\\n\")" + ] + }, { "cell_type": "code", "execution_count": null, @@ -270,39 +364,33 @@ "source": [ "# write the coil out in a format that can be simulated\n", "# shift the coils aronnd to make connections a bit easier\n", - "COIL_ROTATION = -360 / 12\n", + "pf = rotate(points_f, 90)\n", + "pb = rotate(points_b, 90)\n", + "fname = \"simulations/coils/coil_12_custom\"\n", + "if USE_SPIRAL:\n", + " fname = \"simulations/coils/coil_12_spiral\"\n", "\n", - "for i in range(12):\n", - " angle = i * 360 / 12 + COIL_ROTATION\n", - " p_f = scale(translate(rotate(points_f, angle), COIL_CENTER_RADIUS, angle), 0.1)\n", - " p_b = scale(translate(rotate(points_b, angle), COIL_CENTER_RADIUS, angle), 0.1)\n", + "with open(fname + \".csv\", \"w\") as f:\n", + " for point in pf:\n", + " f.write(f\"{point[0]/10},{point[1]/10},0,0.5\\n\")\n", "\n", - " # we'll have positive current through the A coils, and negative current through the B coils\n", - " current = 0\n", - " if i % 3 == 0:\n", - " current = 0.5\n", - " elif i % 3 == 1:\n", - " current = -0.5\n", + "# two layer board\n", + "with open(fname + \"-2-layer.csv\", \"wt\") as f:\n", + " for point in pf[::-1]:\n", + " f.write(f\"{point[0]/10},{point[1]/10},0,0.5\\n\")\n", + " for point in pb:\n", + " f.write(f\"{point[0]/10},{point[1]/10},0-0.062,0.5\\n\")\n", "\n", - " # two layer board\n", - " with open(f\"simulations/coil-2-{i}.txt\", \"wt\") as f:\n", - " for point in p_f[::-1]:\n", - " f.write(f\"{point[0]},{point[1]},0, {current}\\n\")\n", - " for point in p_b:\n", - " f.write(f\"{point[0]},{point[1]},0-0.62,{current}\\n\")\n", - "\n", - " # four layer board\n", - " with open(f\"simulations/coil-4-top-{i}.txt\", \"wt\") as f:\n", - " for point in p_f[::-1]:\n", - " f.write(f\"{point[0]},{point[1]},0,{current}\\n\")\n", - " for point in p_b:\n", - " f.write(f\"{point[0]},{point[1]},0-0.11,{current}\\n\")\n", - "\n", - " with open(f\"simulations/coil-4-bottom-{i}.txt\", \"wt\") as f:\n", - " for point in p_f[::-1]:\n", - " f.write(f\"{point[0]},{point[1]},0-(0.11+0.4),{current}\\n\")\n", - " for point in p_b:\n", - " f.write(f\"{point[0]},{point[1]},0-(0.11+0.11+0.4),{current}\\n\")" + "# all four layer board\n", + "with open(fname + \"-4-layer.csv\", \"wt\") as f:\n", + " for point in pf[::-1]:\n", + " f.write(f\"{point[0]/10},{point[1]/10},0,0.5\\n\")\n", + " for point in pb:\n", + " f.write(f\"{point[0]/10},{point[1]/10},0-0.011,0.5\\n\")\n", + " for point in pf[::-1]:\n", + " f.write(f\"{point[0]/10},{point[1]/10},0-(0.011+0.04),0.5\\n\")\n", + " for point in pb:\n", + " f.write(f\"{point[0]/10},{point[1]/10},0-(0.011+0.011+0.04),0.5\\n\")" ] }, { @@ -538,35 +626,35 @@ "\n", "outer_cuts = (\n", " draw_arc(-45 + nibble_angle_size / 2, 45 - nibble_angle_size / 2, STATOR_RADIUS, 5)\n", - " # + translate(\n", - " # rotate(draw_arc(1, 179, SCREW_HOLE_DRILL_DIAM / 2, 5)[::-1], 135),\n", - " # STATOR_RADIUS,\n", - " # 45,\n", - " # )\n", + " + translate(\n", + " rotate(draw_arc(5, 175, SCREW_HOLE_DRILL_DIAM / 2, 5)[::-1], 135),\n", + " STATOR_RADIUS,\n", + " 45,\n", + " )\n", " + draw_arc(\n", " 45 + nibble_angle_size / 2, 135 - nibble_angle_size / 2, STATOR_RADIUS, 5\n", " )\n", - " # + translate(\n", - " # rotate(draw_arc(1, 179, SCREW_HOLE_DRILL_DIAM / 2, 5), 225)[::-1],\n", - " # STATOR_RADIUS,\n", - " # 135,\n", - " # )\n", + " + translate(\n", + " rotate(draw_arc(5, 175, SCREW_HOLE_DRILL_DIAM / 2, 5), 225)[::-1],\n", + " STATOR_RADIUS,\n", + " 135,\n", + " )\n", " + draw_arc(\n", " 135 + nibble_angle_size / 2, 225 - nibble_angle_size / 2, STATOR_RADIUS, 5\n", " )\n", - " # + translate(\n", - " # rotate(draw_arc(1, 179, SCREW_HOLE_DRILL_DIAM / 2, 5), 315)[::-1],\n", - " # STATOR_RADIUS,\n", - " # 225,\n", - " # )\n", + " + translate(\n", + " rotate(draw_arc(5, 175, SCREW_HOLE_DRILL_DIAM / 2, 5), 315)[::-1],\n", + " STATOR_RADIUS,\n", + " 225,\n", + " )\n", " + draw_arc(\n", " 225 + nibble_angle_size / 2, 315 - nibble_angle_size / 2, STATOR_RADIUS, 5\n", " )\n", - " # + translate(\n", - " # rotate(draw_arc(1, 179, SCREW_HOLE_DRILL_DIAM / 2, 5), 45)[::-1],\n", - " # STATOR_RADIUS,\n", - " # 315,\n", - " # )\n", + " + translate(\n", + " rotate(draw_arc(5, 175, SCREW_HOLE_DRILL_DIAM / 2, 5), 45)[::-1],\n", + " STATOR_RADIUS,\n", + " 315,\n", + " )\n", ")\n", "\n", "edge_cuts = [\n", diff --git a/coil_generator-6.ipynb b/coil_generator-6.ipynb index e6ad2f9..daf459a 100644 --- a/coil_generator-6.ipynb +++ b/coil_generator-6.ipynb @@ -65,6 +65,7 @@ "# PCB Edge size\n", "STATOR_RADIUS = 25\n", "STATOR_HOLE_RADIUS = 5.5\n", + "HOLE_SPACE = 2\n", "\n", "# where to puth the mounting pins\n", "SCREW_HOLE_DRILL_DIAM = 2.3 # 2.3mm drill for a 2mm screw\n", @@ -72,16 +73,41 @@ "\n", "HOLE_SPACING = 0.25\n", "\n", - "# where to put the input pads\n", - "INPUT_PAD_RADIUS = STATOR_RADIUS - (PAD_WIDTH / 2 + VIA_DIAM + TRACK_SPACING)\n", - "\n", "# Coil params\n", "TURNS = 26\n", "COIL_CENTER_RADIUS = 15\n", - "COIL_VIA_RADIUS = 16\n", + "COIL_VIA_RADIUS = 15.3" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Large 30 mm version\n", "\n", - "# where to place the pins\n", - "# CONNECTION_PINS_RADIUS = 16.5\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 = 31\n", + "# COIL_CENTER_RADIUS = 19\n", + "# COIL_VIA_RADIUS = 19.3" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# where to put the input pads\n", + "INPUT_PAD_RADIUS = STATOR_RADIUS - (PAD_WIDTH / 2 + VIA_DIAM + TRACK_SPACING)\n", "\n", "USE_SPIRAL = False\n", "\n", @@ -218,8 +244,8 @@ " optimize_points(get_points(template, TURNS, TRACK_SPACING + TRACK_WIDTH)), 2\n", " )\n", "\n", - " points_f = [(0, 0)] + points_f\n", - " points_b = [(0, 0)] + points_b\n", + " points_f = [(COIL_VIA_RADIUS - COIL_CENTER_RADIUS, 0)] + points_f\n", + " points_b = [(COIL_VIA_RADIUS - COIL_CENTER_RADIUS, 0)] + points_b\n", "\n", " df = pd.DataFrame(points_f, columns=[\"x\", \"y\"])\n", " ax = df.plot.line(x=\"x\", y=\"y\", color=\"blue\")\n", @@ -237,27 +263,7 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "# write the coil out in a format that can be simulated\n", - "\n", - "# two layer board\n", - "with open(\"simulations/coil-2.txt\", \"wt\") as f:\n", - " for point in points_f[::-1]:\n", - " f.write(f\"{point[0]},{point[1]},0,0.5\\n\")\n", - " for point in points_b:\n", - " f.write(f\"{point[0]},{point[1]},0-0.62,0.5\\n\")\n", - "\n", - "# four layer board\n", - "with open(\"simulations/coil-4.txt\", \"wt\") as f:\n", - " for point in points_f[::-1]:\n", - " f.write(f\"{point[0]},{point[1]},0,0.5\\n\")\n", - " for point in points_b:\n", - " f.write(f\"{point[0]},{point[1]},0-0.11,0.5\\n\")\n", - " for point in points_f[::-1]:\n", - " f.write(f\"{point[0]},{point[1]},0-(0.11+0.4),0.5\\n\")\n", - " for point in points_b:\n", - " f.write(f\"{point[0]},{point[1]},0-(0.11+0.11+0.4),0.5\\n\")" - ] + "source": [] }, { "cell_type": "markdown", @@ -309,6 +315,50 @@ " print(\"Using template\")" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Write out coils for simulation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# write the coil out in a format that can be simulated\n", + "# rotate the points by 90 degrees so that the x axis is horizontal\n", + "pf = rotate(points_f, 90)\n", + "pb = rotate(points_b, 90)\n", + "fname = \"simulations/coils/coil_6_custom\"\n", + "if USE_SPIRAL:\n", + " fname = \"simulations/coils/coil_6_spiral\"\n", + "\n", + "with open(fname + \".csv\", \"w\") as f:\n", + " for point in pf:\n", + " f.write(f\"{point[0]/10},{point[1]/10},0,0.5\\n\")\n", + "\n", + "# two layer board\n", + "with open(fname + \"-2-layer.csv\", \"wt\") as f:\n", + " for point in pf[::-1]:\n", + " f.write(f\"{point[0]/10},{point[1]/10},0,0.5\\n\")\n", + " for point in pb:\n", + " f.write(f\"{point[0]/10},{point[1]/10},0-0.062,0.5\\n\")\n", + "\n", + "# all four layer board\n", + "with open(fname + \"-4-layer.csv\", \"wt\") as f:\n", + " for point in pf[::-1]:\n", + " f.write(f\"{point[0]/10},{point[1]/10},0,0.5\\n\")\n", + " for point in pb:\n", + " f.write(f\"{point[0]/10},{point[1]/10},0-0.011,0.5\\n\")\n", + " for point in pf[::-1]:\n", + " f.write(f\"{point[0]/10},{point[1]/10},0-(0.011+0.04),0.5\\n\")\n", + " for point in pb:\n", + " f.write(f\"{point[0]/10},{point[1]/10},0-(0.011+0.011+0.04),0.5\\n\")" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -403,12 +453,12 @@ "tracks_b.append(coil_C_opp_b)\n", "\n", "# connect the front and back coils together\n", - "vias.append(create_via(get_arc_point(angle_A, COIL_CENTER_RADIUS)))\n", - "vias.append(create_via(get_arc_point(angle_B, COIL_CENTER_RADIUS)))\n", - "vias.append(create_via(get_arc_point(angle_C, COIL_CENTER_RADIUS)))\n", - "vias.append(create_via(get_arc_point(angle_A_opp, COIL_CENTER_RADIUS)))\n", - "vias.append(create_via(get_arc_point(angle_B_opp, COIL_CENTER_RADIUS)))\n", - "vias.append(create_via(get_arc_point(angle_C_opp, COIL_CENTER_RADIUS)))\n", + "vias.append(create_via(get_arc_point(angle_A, COIL_VIA_RADIUS)))\n", + "vias.append(create_via(get_arc_point(angle_B, COIL_VIA_RADIUS)))\n", + "vias.append(create_via(get_arc_point(angle_C, COIL_VIA_RADIUS)))\n", + "vias.append(create_via(get_arc_point(angle_A_opp, COIL_VIA_RADIUS)))\n", + "vias.append(create_via(get_arc_point(angle_B_opp, COIL_VIA_RADIUS)))\n", + "vias.append(create_via(get_arc_point(angle_C_opp, COIL_VIA_RADIUS)))\n", "\n", "# connect the front copper opposite coils together\n", "common_connection_radius = SCREW_HOLE_RADIUS - (\n", @@ -426,7 +476,7 @@ "# vias.append(create_via(get_arc_point(angle_C_opp, common_connection_radius)))\n", "\n", "# wires for connecting to opposite coils\n", - "connection_radius1 = STATOR_HOLE_RADIUS + (2 * TRACK_SPACING)\n", + "connection_radius1 = STATOR_HOLE_RADIUS + (2 * TRACK_SPACING + HOLE_SPACING)\n", "connection_radius2 = connection_radius1 + (2 * TRACK_SPACING + VIA_DIAM / 2)\n", "\n", "# draw a 45 degree line from each coil at connection radius 1\n", @@ -556,7 +606,7 @@ "outer_cuts = (\n", " draw_arc(-45 + nibble_angle_size / 2, 45 - nibble_angle_size / 2, STATOR_RADIUS, 5)\n", " + translate(\n", - " rotate(draw_arc(1, 179, SCREW_HOLE_DRILL_DIAM / 2, 5)[::-1], 135),\n", + " rotate(draw_arc(5, 175, SCREW_HOLE_DRILL_DIAM / 2, 5)[::-1], 135),\n", " STATOR_RADIUS,\n", " 45,\n", " )\n", @@ -564,7 +614,7 @@ " 45 + nibble_angle_size / 2, 135 - nibble_angle_size / 2, STATOR_RADIUS, 5\n", " )\n", " + translate(\n", - " rotate(draw_arc(1, 179, SCREW_HOLE_DRILL_DIAM / 2, 5), 225)[::-1],\n", + " rotate(draw_arc(5, 175, SCREW_HOLE_DRILL_DIAM / 2, 5), 225)[::-1],\n", " STATOR_RADIUS,\n", " 135,\n", " )\n", @@ -572,7 +622,7 @@ " 135 + nibble_angle_size / 2, 225 - nibble_angle_size / 2, STATOR_RADIUS, 5\n", " )\n", " + translate(\n", - " rotate(draw_arc(1, 179, SCREW_HOLE_DRILL_DIAM / 2, 5), 315)[::-1],\n", + " rotate(draw_arc(5, 175, SCREW_HOLE_DRILL_DIAM / 2, 5), 315)[::-1],\n", " STATOR_RADIUS,\n", " 225,\n", " )\n", @@ -580,7 +630,7 @@ " 225 + nibble_angle_size / 2, 315 - nibble_angle_size / 2, STATOR_RADIUS, 5\n", " )\n", " + translate(\n", - " rotate(draw_arc(1, 179, SCREW_HOLE_DRILL_DIAM / 2, 5), 45)[::-1],\n", + " rotate(draw_arc(5, 175, SCREW_HOLE_DRILL_DIAM / 2, 5), 45)[::-1],\n", " STATOR_RADIUS,\n", " 315,\n", " )\n", @@ -593,7 +643,7 @@ "\n", "# dump out the json version\n", "json_result = dump_json(\n", - " filename=\"coils_6.json\",\n", + " filename=f\"coils_6_{STATOR_RADIUS}.json\",\n", " track_width=TRACK_WIDTH,\n", " pin_diam=PIN_DIAM,\n", " pin_drill=PIN_DRILL,\n", diff --git a/coil_generator-radial.ipynb b/coil_generator-radial.ipynb new file mode 100644 index 0000000..ed599ea --- /dev/null +++ b/coil_generator-radial.ipynb @@ -0,0 +1,579 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import pandas as pd\n", + "import numpy as np\n", + "\n", + "import matplotlib.pyplot as plt\n", + "from skspatial.objects import LineSegment, Line, Vector\n", + "\n", + "# some helper functions\n", + "from helpers import (\n", + " get_arc_point,\n", + " draw_arc,\n", + " rotate,\n", + " translate,\n", + " flip_y,\n", + " flip_x,\n", + " optimize_points,\n", + " chaikin,\n", + " rotate_point,\n", + " scale,\n", + ")\n", + "from pcb_json import (\n", + " dump_json,\n", + " plot_json,\n", + " create_via,\n", + " create_pad,\n", + " create_pin,\n", + " create_track,\n", + " create_silk,\n", + " create_silk,\n", + " create_mounting_hole,\n", + ")\n", + "\n", + "from enum import Enum\n", + "\n", + "Layer = Enum(\"Layer\", \"FRONT BACK\")" + ] + }, + { + "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" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Standard 25 mm version\n", + "\n", + "# PCB Edge size\n", + "STATOR_RADIUS = 25\n", + "STATOR_HOLE_RADIUS = 5.5\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 = 12\n", + "COIL_CENTER_RADIUS = 16\n", + "COIL_VIA_RADIUS = 17" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Large 30 mm version\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" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# where to put the input pads\n", + "INPUT_PAD_RADIUS = STATOR_RADIUS - (PAD_WIDTH / 2 + VIA_DIAM + TRACK_SPACING)\n", + "\n", + "USE_SPIRAL = False\n", + "\n", + "LAYERS = 4" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Radial coil generator" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "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", + " points.append(get_arc_point(angle, inner_radius))\n", + " points.append(get_arc_point(angle, outer_radius))\n", + " points.append(\n", + " (\n", + " get_arc_point(angle, outer_radius)[0],\n", + " get_arc_point(angle, outer_radius)[1] + spacing,\n", + " )\n", + " )\n", + " points.append(get_arc_point(angle + spacing_angle, inner_radius))\n", + " return points" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "points_f = get_points(\n", + " TRACK_SPACING + TRACK_WIDTH,\n", + " (COIL_CENTER_RADIUS - 10),\n", + " COIL_CENTER_RADIUS + 10,\n", + " -15,\n", + " 15,\n", + ")\n", + "\n", + "df = pd.DataFrame(points_f, columns=[\"x\", \"y\"])\n", + "ax = df.plot.line(x=\"x\", y=\"y\", color=\"blue\")\n", + "ax.axis(\"equal\")\n", + "\n", + "# calculat the total length of the track to compute the resistance\n", + "total_length_front = 0\n", + "for i in range(len(points_f) - 1):\n", + " total_length_front += np.linalg.norm(\n", + " np.array(points_f[i + 1]) - np.array(points_f[i])\n", + " )\n", + "print(\"Total length front\", total_length_front)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "spacing = TRACK_SPACING + TRACK_WIDTH\n", + "inner_radius = COIL_CENTER_RADIUS - 10\n", + "outer_radius = COIL_CENTER_RADIUS + 10\n", + "start_angle = -15\n", + "end_angle = 15\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", + "# now calculate the points be iterating from start_angle to end_angle with the spacing_angle\n", + "i = 0\n", + "count = (end_angle - start_angle) // spacing_angle\n", + "print(spacing_angle, count)\n", + "# for angle in np.arange(start_angle, end_angle, spacing_angle):\n", + "for angle in [start_angle, end_angle]:\n", + " i = i + 1\n", + " current = 1\n", + " if i == 2:\n", + " current = -1\n", + " with open(f\"simulations/radial_{i}.csv\", \"w\") as f:\n", + " points = [\n", + " get_arc_point(angle, inner_radius),\n", + " get_arc_point(angle, outer_radius),\n", + " ]\n", + " points = translate(rotate(points, 90), -22.5, 90)\n", + " for point in points:\n", + " f.write(f\"{point[0]}/10,{point[1]}/10,0,{current}\\n\")\n", + "\n", + "\n", + "# # write the coil out in a format that can be simulated\n", + "p_f = translate(rotate(points_f, 90), -22.5, 90)\n", + "p_b = translate(rotate(points_f, 90), -22.5, 90)\n", + "\n", + "# df = pd.DataFrame(p_f, columns=[\"x\", \"y\"])\n", + "# ax = df.plot.line(x=\"x\", y=\"y\", color=\"blue\")\n", + "# ax.axis(\"equal\")\n", + "\n", + "with open(\"simulations/coils/coil_rad.csv\", \"w\") as f:\n", + " for point in p_f:\n", + " f.write(f\"{point[0]},{point[1]},0,0.5\\n\")\n", + "\n", + "# two layer board\n", + "with open(\"simulations/coils/coil-rad-2-layer.csv\", \"wt\") as f:\n", + " for point in p_f:\n", + " f.write(f\"{point[0]},{point[1]},0,0.5\\n\")\n", + " for point in p_b:\n", + " f.write(f\"{point[0]},{point[1]},{0-0.062},0.5\\n\")\n", + "\n", + "# all four layer board\n", + "with open(\"simulations/coils/coil-rad-4-layer.csv\", \"wt\") as f:\n", + " for point in p_f:\n", + " f.write(f\"{point[0]},{point[1]},0,0.5\\n\")\n", + " for point in p_b:\n", + " f.write(f\"{point[0]},{point[1]},{0-0.011},0.5\\n\")\n", + " for point in p_f:\n", + " f.write(f\"{point[0]},{point[1]},{0-(0.011+0.04)},0.5\\n\")\n", + " for point in p_b:\n", + " f.write(f\"{point[0]},{point[1]},{0-(0.011+0.011+0.04)},0.5\\n\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Generate PCB Layout" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# calculat the total length of the track to compute the resistance\n", + "total_length_front = 0\n", + "for i in range(len(points_f) - 1):\n", + " total_length_front += np.linalg.norm(\n", + " np.array(points_f[i + 1]) - np.array(points_f[i])\n", + " )\n", + "print(\"Total length front\", total_length_front)\n", + "\n", + "total_length_back = 0\n", + "for i in range(len(points_b) - 1):\n", + " total_length_back += np.linalg.norm(\n", + " np.array(points_b[i + 1]) - np.array(points_b[i])\n", + " )\n", + "print(\"Total length back\", total_length_back)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "vias = []\n", + "tracks_f = []\n", + "tracks_b = []\n", + "pads = []\n", + "pins = []\n", + "mounting_holes = []\n", + "silk = []\n", + "\n", + "\n", + "# shift the coils aronnd to make connections a bit easier\n", + "COIL_ROTATION = -360 / 12\n", + "\n", + "coil_angles = []\n", + "for i in range(12):\n", + " angle = i * 360 / 12 + COIL_ROTATION\n", + " coil_angles.append(angle)\n", + "\n", + "# the main coils\n", + "coil_labels = [\"A\", \"B\", \"C\"]\n", + "coils_f = []\n", + "coils_b = []\n", + "for i in range(12):\n", + " angle = coil_angles[i]\n", + " if (i // 3) % 2 == 0:\n", + " coil_A_f = translate(rotate(points_f, angle), COIL_CENTER_RADIUS, angle)\n", + " coil_A_b = translate(rotate(points_b, angle), COIL_CENTER_RADIUS, angle)\n", + " else:\n", + " # slightly nudge the coils so that they don't overlap when flipped\n", + " coil_A_f = translate(rotate(flip_y(points_f), angle), COIL_CENTER_RADIUS, angle)\n", + " coil_A_b = translate(rotate(flip_y(points_b), angle), COIL_CENTER_RADIUS, angle)\n", + " # keep track of the coils\n", + " coils_f.append(coil_A_f)\n", + " coils_b.append(coil_A_b)\n", + "\n", + " tracks_f.append(coil_A_f)\n", + " tracks_b.append(coil_A_b)\n", + " vias.append(create_via(get_arc_point(angle, COIL_VIA_RADIUS)))\n", + " silk.append(\n", + " create_silk(get_arc_point(angle, COIL_CENTER_RADIUS), coil_labels[i % 3])\n", + " )\n", + "\n", + "# raidus for connecting the bottoms of the coils together\n", + "connection_radius1 = STATOR_HOLE_RADIUS + 3 * TRACK_SPACING\n", + "\n", + "# create tracks to link the A coils around the center\n", + "connection_via_radius_A = connection_radius1 + 3 * TRACK_SPACING + VIA_DIAM / 2\n", + "coil_A1_A2_inner = (\n", + " [get_arc_point(coil_angles[0], connection_via_radius_A)]\n", + " + draw_arc(COIL_ROTATION, coil_angles[3], connection_radius1)\n", + " + [get_arc_point(coil_angles[3], connection_via_radius_A)]\n", + ")\n", + "tracks_f.append(coil_A1_A2_inner)\n", + "coil_A3_A4_inner = (\n", + " [get_arc_point(coil_angles[6], connection_via_radius_A)]\n", + " + draw_arc(coil_angles[6], coil_angles[9], connection_radius1)\n", + " + [get_arc_point(coil_angles[9], connection_via_radius_A)]\n", + ")\n", + "tracks_f.append(coil_A3_A4_inner)\n", + "# connect up the bottoms of the A coils\n", + "coils_b[0].append(coil_A1_A2_inner[0])\n", + "coils_b[3].append(coil_A1_A2_inner[-1])\n", + "coils_b[6].append(coil_A3_A4_inner[0])\n", + "coils_b[9].append(coil_A3_A4_inner[-1])\n", + "# add the vias to stitch them together\n", + "vias.append(create_via(coil_A1_A2_inner[0]))\n", + "vias.append(create_via(coil_A1_A2_inner[-1]))\n", + "vias.append(create_via(coil_A3_A4_inner[0]))\n", + "vias.append(create_via(coil_A3_A4_inner[-1]))\n", + "\n", + "# create tracks to link the B coils around the center - this can all be done on the bottom layer\n", + "coil_B1_B2_inner = draw_arc(coil_angles[1], coil_angles[4], connection_radius1)\n", + "tracks_b.append(coil_B1_B2_inner)\n", + "coil_B3_B4_inner = draw_arc(coil_angles[7], coil_angles[10], connection_radius1)\n", + "tracks_b.append(coil_B3_B4_inner)\n", + "# connect up the bottoms of the A coils\n", + "coils_b[1].append(coil_B1_B2_inner[0])\n", + "coils_b[4].append(coil_B1_B2_inner[-1])\n", + "coils_b[7].append(coil_B3_B4_inner[0])\n", + "coils_b[10].append(coil_B3_B4_inner[-1])\n", + "\n", + "# create tracks to link the C coils around the center\n", + "connection_via_radius_C = connection_via_radius_A + 3 * TRACK_SPACING + VIA_DIAM / 2\n", + "coil_C1_C2_inner = draw_arc(coil_angles[2], coil_angles[5], connection_via_radius_C)\n", + "tracks_f.append(coil_C1_C2_inner)\n", + "coil_C3_C4_inner = draw_arc(coil_angles[8], coil_angles[11], connection_via_radius_C)\n", + "tracks_f.append(coil_C3_C4_inner)\n", + "# connect up the bottoms of the B coils\n", + "coils_b[2].append(coil_C1_C2_inner[0])\n", + "coils_b[5].append(coil_C1_C2_inner[-1])\n", + "coils_b[8].append(coil_C3_C4_inner[0])\n", + "coils_b[11].append(coil_C3_C4_inner[-1])\n", + "# add the vias to stitch them together\n", + "vias.append(create_via(coil_C1_C2_inner[0]))\n", + "vias.append(create_via(coil_C1_C2_inner[-1]))\n", + "vias.append(create_via(coil_C3_C4_inner[0]))\n", + "vias.append(create_via(coil_C3_C4_inner[-1]))\n", + "\n", + "# connect the last three coils together\n", + "common_connection_radius = SCREW_HOLE_RADIUS - (SCREW_HOLE_DRILL_DIAM / 2 + 0.5)\n", + "tracks_f.append(draw_arc(coil_angles[9], coil_angles[11], common_connection_radius))\n", + "coils_f[9].append(get_arc_point(coil_angles[9], common_connection_radius))\n", + "coils_f[10].append(get_arc_point(coil_angles[10], common_connection_radius))\n", + "coils_f[11].append(get_arc_point(coil_angles[11], common_connection_radius))\n", + "\n", + "# connect the outer A coils together\n", + "outer_connection_radius_A = SCREW_HOLE_RADIUS - (SCREW_HOLE_DRILL_DIAM / 2 + 0.5)\n", + "tracks_f.append(draw_arc(coil_angles[3], coil_angles[6], outer_connection_radius_A))\n", + "coils_f[3].append(get_arc_point(coil_angles[3], outer_connection_radius_A))\n", + "coils_f[6].append(get_arc_point(coil_angles[6], outer_connection_radius_A))\n", + "\n", + "# connect the outer B coils together\n", + "outer_connection_radius_B = outer_connection_radius_A - TRACK_SPACING - VIA_DIAM / 2\n", + "tracks_b.append(\n", + " [get_arc_point(coil_angles[4], outer_connection_radius_B)]\n", + " + draw_arc(coil_angles[4], coil_angles[7], outer_connection_radius_A)\n", + " + [get_arc_point(coil_angles[7], outer_connection_radius_B)]\n", + ")\n", + "coils_f[4].append(get_arc_point(coil_angles[4], outer_connection_radius_B))\n", + "coils_f[7].append(get_arc_point(coil_angles[7], outer_connection_radius_B))\n", + "vias.append(\n", + " create_via(get_arc_point(4 * 360 / 12 + COIL_ROTATION, outer_connection_radius_B))\n", + ")\n", + "vias.append(\n", + " create_via(get_arc_point(7 * 360 / 12 + COIL_ROTATION, outer_connection_radius_B))\n", + ")\n", + "\n", + "# connect the outer C coilds together\n", + "outer_connection_radius_C = outer_connection_radius_B - TRACK_SPACING - VIA_DIAM / 2\n", + "tracks_b.append(\n", + " draw_arc(\n", + " 5 * 360 / 12 + COIL_ROTATION,\n", + " 8 * 360 / 12 + COIL_ROTATION,\n", + " outer_connection_radius_C,\n", + " )\n", + ")\n", + "coils_f[5].append(\n", + " get_arc_point(5 * 360 / 12 + COIL_ROTATION, outer_connection_radius_C)\n", + ")\n", + "coils_f[8].append(\n", + " get_arc_point(8 * 360 / 12 + COIL_ROTATION, outer_connection_radius_C)\n", + ")\n", + "vias.append(\n", + " create_via(get_arc_point(5 * 360 / 12 + COIL_ROTATION, outer_connection_radius_C))\n", + ")\n", + "vias.append(\n", + " create_via(get_arc_point(8 * 360 / 12 + COIL_ROTATION, outer_connection_radius_C))\n", + ")\n", + "\n", + "# create the pads for connecting the inputs to the coils\n", + "silk.append(\n", + " create_silk((INPUT_PAD_RADIUS - PAD_HEIGHT - 2.5, PAD_PITCH), \"C\", \"b\", 2.5, -900)\n", + ")\n", + "silk.append(create_silk((INPUT_PAD_RADIUS - PAD_HEIGHT - 2.5, 0), \"B\", \"b\", 2.5, -900))\n", + "silk.append(\n", + " create_silk((INPUT_PAD_RADIUS - PAD_HEIGHT - 2.5, -PAD_PITCH), \"A\", \"b\", 2.5, -900)\n", + ")\n", + "\n", + "pads.append(create_pad((INPUT_PAD_RADIUS, -PAD_PITCH), PAD_WIDTH, PAD_HEIGHT, \"b\"))\n", + "pads.append(create_pad((INPUT_PAD_RADIUS, 0), PAD_WIDTH, PAD_HEIGHT, \"b\"))\n", + "pads.append(create_pad((INPUT_PAD_RADIUS, PAD_PITCH), PAD_WIDTH, PAD_HEIGHT, \"b\"))\n", + "\n", + "# connect coil A to the top pad\n", + "pad_connection_point_x = INPUT_PAD_RADIUS\n", + "pad_angle = np.rad2deg(np.arcsin(PAD_PITCH / pad_connection_point_x))\n", + "coils_f[0].append(get_arc_point(coil_angles[0], pad_connection_point_x))\n", + "vias.append(create_via(get_arc_point(coil_angles[0], pad_connection_point_x)))\n", + "# connect coil B to the middle pad\n", + "coils_f[1].append((pad_connection_point_x + PAD_WIDTH / 2 + VIA_DIAM / 2, 0))\n", + "vias.append(create_via(((pad_connection_point_x + PAD_WIDTH / 2 + VIA_DIAM / 2, 0))))\n", + "# connect coil C to the bottom pad\n", + "coils_f[2].append(get_arc_point(coil_angles[2], pad_connection_point_x))\n", + "vias.append(create_via(get_arc_point(coil_angles[2], pad_connection_point_x)))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# if we are doing four layers then duplicate the front and back layers on (front and inner1), (inner2 and back)\n", + "tracks_in1 = []\n", + "tracks_in2 = []\n", + "if LAYERS == 4:\n", + " tracks_in1 = tracks_b.copy()\n", + " tracks_in2 = tracks_f.copy()\n", + "\n", + "# these final bits of wiring up to the input pads don't need to be duplicated\n", + "tracks_b.append(\n", + " [(pad_connection_point_x + PAD_WIDTH / 2, 0), (pad_connection_point_x, 0)]\n", + ")\n", + "tracks_b.append(draw_arc(coil_angles[0], -pad_angle, pad_connection_point_x, 1))\n", + "tracks_b.append(draw_arc(coil_angles[2], pad_angle, pad_connection_point_x, 1))\n", + "\n", + "nibble_angle_size = 360 * SCREW_HOLE_DRILL_DIAM / (2 * np.pi * STATOR_RADIUS)\n", + "\n", + "outer_cuts = (\n", + " draw_arc(-45 + nibble_angle_size / 2, 45 - nibble_angle_size / 2, STATOR_RADIUS, 5)\n", + " + translate(\n", + " rotate(draw_arc(5, 175, SCREW_HOLE_DRILL_DIAM / 2, 5)[::-1], 135),\n", + " STATOR_RADIUS,\n", + " 45,\n", + " )\n", + " + draw_arc(\n", + " 45 + nibble_angle_size / 2, 135 - nibble_angle_size / 2, STATOR_RADIUS, 5\n", + " )\n", + " + translate(\n", + " rotate(draw_arc(5, 175, SCREW_HOLE_DRILL_DIAM / 2, 5), 225)[::-1],\n", + " STATOR_RADIUS,\n", + " 135,\n", + " )\n", + " + draw_arc(\n", + " 135 + nibble_angle_size / 2, 225 - nibble_angle_size / 2, STATOR_RADIUS, 5\n", + " )\n", + " + translate(\n", + " rotate(draw_arc(5, 175, SCREW_HOLE_DRILL_DIAM / 2, 5), 315)[::-1],\n", + " STATOR_RADIUS,\n", + " 225,\n", + " )\n", + " + draw_arc(\n", + " 225 + nibble_angle_size / 2, 315 - nibble_angle_size / 2, STATOR_RADIUS, 5\n", + " )\n", + " + translate(\n", + " rotate(draw_arc(5, 175, SCREW_HOLE_DRILL_DIAM / 2, 5), 45)[::-1],\n", + " STATOR_RADIUS,\n", + " 315,\n", + " )\n", + ")\n", + "\n", + "edge_cuts = [\n", + " outer_cuts,\n", + " draw_arc(0, 360, STATOR_HOLE_RADIUS, 1),\n", + "]\n", + "\n", + "# dump out the json version\n", + "json_result = dump_json(\n", + " filename=f\"coils_12_{STATOR_RADIUS}mm.json\",\n", + " track_width=TRACK_WIDTH,\n", + " pin_diam=PIN_DIAM,\n", + " pin_drill=PIN_DRILL,\n", + " via_diam=VIA_DIAM,\n", + " via_drill=VIA_DRILL,\n", + " vias=vias,\n", + " pins=pins,\n", + " pads=pads,\n", + " silk=silk,\n", + " tracks_f=tracks_f,\n", + " tracks_in1=tracks_in1,\n", + " tracks_in2=tracks_in2,\n", + " tracks_b=tracks_b,\n", + " mounting_holes=mounting_holes,\n", + " edge_cuts=edge_cuts,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# plot the json\n", + "plot_json(json_result)" + ] + } + ], + "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/coil_plugin.py b/coil_plugin.py index 30f182b..e1d099e 100644 --- a/coil_plugin.py +++ b/coil_plugin.py @@ -188,5 +188,34 @@ class CoilPlugin(pcbnew.ActionPlugin): ec.SetPolyPoints(v) board.Add(ec) + # put it on the solder mask as well - who knows why... + for edge_cut in coil_data["edgeCuts"]: + ec = pcbnew.PCB_SHAPE(board) + ec.SetShape(pcbnew.SHAPE_T_POLY) + ec.SetFilled(False) + ec.SetLayer(pcbnew.F_Mask) + ec.SetWidth(int(0.1 * pcbnew.IU_PER_MM)) + v = pcbnew.wxPoint_Vector() + for point in edge_cut: + x = point["x"] + CENTER_X + y = point["y"] + CENTER_Y + v.append(pcbnew.wxPointMM(x, y)) + ec.SetPolyPoints(v) + board.Add(ec) + + for edge_cut in coil_data["edgeCuts"]: + ec = pcbnew.PCB_SHAPE(board) + ec.SetShape(pcbnew.SHAPE_T_POLY) + ec.SetFilled(False) + ec.SetLayer(pcbnew.B_Mask) + ec.SetWidth(int(0.1 * pcbnew.IU_PER_MM)) + v = pcbnew.wxPoint_Vector() + for point in edge_cut: + x = point["x"] + CENTER_X + y = point["y"] + CENTER_Y + v.append(pcbnew.wxPointMM(x, y)) + ec.SetPolyPoints(v) + board.Add(ec) + CoilPlugin().register() # Instantiate and register to Pcbnew]) diff --git a/requirements.txt b/requirements.txt index a0de6da..e693d95 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,4 +3,10 @@ jupyter numpy pandas scikit-spatial -pre-commit \ No newline at end of file +pre-commit +tqdm +plotly +magpylib +pyvista +mayavi +PySide \ No newline at end of file diff --git a/simulations/.DS_Store b/simulations/.DS_Store new file mode 100644 index 0000000..13eb6b2 Binary files /dev/null and b/simulations/.DS_Store differ diff --git a/simulations/biot_savart_v4_3.py b/simulations/biot_savart_v4_3.py index 4c542a5..6c92bae 100644 --- a/simulations/biot_savart_v4_3.py +++ b/simulations/biot_savart_v4_3.py @@ -1,4 +1,4 @@ -''' +""" Biot-Savart Magnetic Field Calculator v4.3 Mingde Yin Ryan Zazo @@ -6,7 +6,7 @@ Ryan Zazo All lengths are in cm, B-field is in G from - https://github.com/vuthalab/biot-savart -''' +""" import numpy as np import matplotlib.pyplot as plt @@ -14,17 +14,18 @@ from mpl_toolkits.mplot3d import Axes3D import matplotlib.cm as cm import matplotlib.ticker as ticker -''' +""" Feature Wishlist: improve plot_coil with different colors for different values of current accelerate integrator to use meshgrids directly instead of a layer of for loop get parse_coil to use vectorized function instead of for loop -''' +""" + def parse_coil(filename): - ''' + """ Parses 4 column CSV into x,y,z,I slices for coil. Each (x,y,z,I) entry defines a vertex on the coil. @@ -35,157 +36,219 @@ def parse_coil(filename): - There are 2 amps of current running between points 1 and 2 - There are 3 amps of current running between points 2 and 3 - The last bit of current is functionally useless. - ''' - with open(filename, "r") as f: return np.array([[eval(i) for i in line.split(",")] for line in f.read().splitlines()]).T + """ + with open(filename, "r") as f: + return np.array( + [[eval(i) for i in line.split(",")] for line in f.read().splitlines()] + ).T + def slice_coil(coil, steplength): - ''' + """ Slices a coil into pieces of size steplength. If the coil is already sliced into pieces smaller than that, this does nothing. - ''' + """ + def interpolate_points(p1, p2, parts): - ''' + """ Produces a series of linearly spaced points between two given points in R3+I Linearly interpolates X,Y,Z; but keeps I the SAME i.e. (0, 2, 1, 3), (3, 4, 2, 5), parts = 2: (0, 2, 1, 3), (1.5, 3, 1.5, 3), (3, 4, 2, 5) - ''' - return np.column_stack((np.linspace(p1[0], p2[0], parts+1), np.linspace(p1[1], p2[1], parts+1), - np.linspace(p1[2], p2[2], parts+1), p1[3] * np.ones((parts+1)))) + """ + return np.column_stack( + ( + np.linspace(p1[0], p2[0], parts + 1), + np.linspace(p1[1], p2[1], parts + 1), + np.linspace(p1[2], p2[2], parts + 1), + p1[3] * np.ones((parts + 1)), + ) + ) - newcoil = np.zeros((1, 4)) # fill column with dummy data, we will remove this later. + newcoil = np.zeros( + (1, 4) + ) # fill column with dummy data, we will remove this later. - segment_starts = coil[:,:-1] - segment_ends = coil[:,1:] + segment_starts = coil[:, :-1] + segment_ends = coil[:, 1:] # determine start and end of each segment - segments = segment_ends-segment_starts + segments = segment_ends - segment_starts segment_lengths = np.apply_along_axis(np.linalg.norm, 0, segments) # create segments; determine start and end of each segment, as well as segment lengths # chop up into smaller bits (elements) - stepnumbers = (segment_lengths/steplength).astype(int) + stepnumbers = (segment_lengths / steplength).astype(int) # determine how many steps we must chop each segment into for i in range(segments.shape[1]): - newrows = interpolate_points(segment_starts[:,i], segment_ends[:,i], stepnumbers[i]) + newrows = interpolate_points( + segment_starts[:, i], segment_ends[:, i], stepnumbers[i] + ) # set of new interpolated points to feed in newcoil = np.vstack((newcoil, newrows)) - if newcoil.shape[0] %2 != 0: newcoil = np.vstack((newcoil, newcoil[-1,:])) + if newcoil.shape[0] % 2 != 0: + newcoil = np.vstack((newcoil, newcoil[-1, :])) ## Force the coil to have an even number of segments, for Richardson Extrapolation to work - return newcoil[1:,:].T # return non-dummy columns + return newcoil[1:, :].T # return non-dummy columns + def calculate_field(coil, x, y, z): - ''' + """ Calculates magnetic field vector as a result of some position and current x, y, z, I [In the same coordinate system as the coil] Coil: Input Coil Positions, already sub-divided into small pieces using slice_coil x, y, z: position in cm - + Output B-field is a 3-D vector in units of G - ''' - FACTOR = 0.1 # = mu_0 / 4pi when lengths are in cm, and B-field is in G + """ + FACTOR = 0.1 # = mu_0 / 4pi when lengths are in cm, and B-field is in G def bs_integrate(start, end): - ''' + """ Produces tiny segment of magnetic field vector (dB) using the midpoint approximation over some interval TODO for future optimization: Get this to work with meshgrids directly - ''' - dl = (end-start).T - mid = (start+end)/2 - position = np.array((x-mid[0], y-mid[1], z-mid[2])).T + """ + dl = (end - start).T + mid = (start + end) / 2 + position = np.array((x - mid[0], y - mid[1], z - mid[2])).T # relative position vector - mag = np.sqrt((x-mid[0])**2 + (y-mid[1])**2 + (z-mid[2])**2) + mag = np.sqrt((x - mid[0]) ** 2 + (y - mid[1]) ** 2 + (z - mid[2]) ** 2) # magnitude of the relative position vector - return start[3] * np.cross(dl[:3], position) / np.array((mag ** 3, mag ** 3, mag ** 3)).T + return ( + start[3] + * np.cross(dl[:3], position) + / np.array((mag**3, mag**3, mag**3)).T + ) # Apply the Biot-Savart Law to get the differential magnetic field # current flowing in this segment is represented by start[3] B = 0 # midpoint integration with 1 layer of Richardson Extrapolation - starts, mids, ends = coil[:,:-1:2], coil[:,1::2], coil[:,2::2] + starts, mids, ends = coil[:, :-1:2], coil[:, 1::2], coil[:, 2::2] - for start, mid, end in np.nditer([starts, mids, ends], flags=['external_loop'], order='F'): + for start, mid, end in np.nditer( + [starts, mids, ends], flags=["external_loop"], order="F" + ): # use numpy fast indexing - fullpart = bs_integrate(start, end) # stage 1 richardson - halfpart = bs_integrate(start, mid) + bs_integrate(mid, end) # stage 2 richardson + fullpart = bs_integrate(start, end) # stage 1 richardson + halfpart = bs_integrate(start, mid) + bs_integrate( + mid, end + ) # stage 2 richardson + + B += ( + 4 / 3 * halfpart - 1 / 3 * fullpart + ) # richardson extrapolated midpoint rule + + return ( + B * FACTOR + ) # return SUM of all components as 3 (x,y,z) meshgrids for (Bx, By, Bz) component when evaluated using produce_target_volume - B += 4/3 * halfpart - 1/3 * fullpart # richardson extrapolated midpoint rule - - return B * FACTOR # return SUM of all components as 3 (x,y,z) meshgrids for (Bx, By, Bz) component when evaluated using produce_target_volume def produce_target_volume(coil, box_size, start_point, vol_resolution): - ''' - Generates a set of field vector values for each tuple (x, y, z) in the box. -​ - Coil: Input Coil Positions in format specified above, already sub-divided into small pieces - box_size: (x, y, z) dimensions of the box in cm - start_point: (x, y, z) = (0, 0, 0) = bottom left corner position of the box - vol_resolution: Spatial resolution (in cm) - ''' - x = np.linspace(start_point[0], box_size[0] + start_point[0],int(box_size[0]/vol_resolution)+1) - y = np.linspace(start_point[1], box_size[1] + start_point[1],int(box_size[1]/vol_resolution)+1) - z = np.linspace(start_point[2], box_size[2] + start_point[2],int(box_size[2]/vol_resolution)+1) + """ + Generates a set of field vector values for each tuple (x, y, z) in the box. + ​ + Coil: Input Coil Positions in format specified above, already sub-divided into small pieces + box_size: (x, y, z) dimensions of the box in cm + start_point: (x, y, z) = (0, 0, 0) = bottom left corner position of the box + vol_resolution: Spatial resolution (in cm) + """ + x = np.linspace( + start_point[0], + box_size[0] + start_point[0], + int(box_size[0] / vol_resolution) + 1, + ) + y = np.linspace( + start_point[1], + box_size[1] + start_point[1], + int(box_size[1] / vol_resolution) + 1, + ) + z = np.linspace( + start_point[2], + box_size[2] + start_point[2], + int(box_size[2] / vol_resolution) + 1, + ) # Generate points at regular spacing, incl. end points - - Z, Y, X = np.meshgrid(z, y, x, indexing='ij') + + Z, Y, X = np.meshgrid(z, y, x, indexing="ij") # NOTE: Requires axes to be flipped in order for meshgrid to have the correct dimensional order - return calculate_field(coil, X,Y,Z) + return calculate_field(coil, X, Y, Z) + def get_field_vector(targetVolume, position, start_point, volume_resolution): - ''' + """ Returns the B vector [Bx, By, Bz] components in a generated Target Volume at a given position tuple (x, y, z) in a coordinate system start_point: (x, y, z) = (0, 0, 0) = bottom left corner position of the box volume_resolution: Division of volumetric meshgrid (generate a point every volume_resolution cm) - ''' - relativePosition = ((np.array(position) - np.array(start_point)) / volume_resolution).astype(int) + """ + relativePosition = ( + (np.array(position) - np.array(start_point)) / volume_resolution + ).astype(int) # adjust to the meshgrid's system - if (relativePosition < 0).any(): return ("ERROR: Out of bounds! (negative indices)") + if (relativePosition < 0).any(): + return "ERROR: Out of bounds! (negative indices)" - try: return targetVolume[relativePosition[0], relativePosition[1], relativePosition[2], :] - except: return ("ERROR: Out of bounds!") + try: + return targetVolume[ + relativePosition[0], relativePosition[1], relativePosition[2], : + ] + except: + return "ERROR: Out of bounds!" # basic error checking to see if you actually got a correct input/output -''' + +""" - If you are indexing a targetvolume meshgrid on your own, remember to account for the offset (starting point), and spatial resolution - You will need an index like -''' +""" -def write_target_volume(input_filename,output_filename, box_size, start_point, - coil_resolution=1, volume_resolution=1): - ''' + +def write_target_volume( + input_filename, + output_filename, + box_size, + start_point, + coil_resolution=1, + volume_resolution=1, +): + """ Takes a coil specified in input_filename, generates a target volume, and saves the generated target volume to output_filename. box_size: (x, y, z) dimensions of the box in cm start_point: (x, y, z) = (0, 0, 0) = bottom left corner position of the box AKA the offset coil_resolution: How long each coil subsegment should be volume_resolution: Division of volumetric meshgrid (generate a point every volume_resolution cm) - ''' - coil = parse_coil(input_filename) + """ + coil = parse_coil(input_filename) chopped = slice_coil(coil, coil_resolution) - targetVolume = produce_target_volume(chopped, box_size, start_point, volume_resolution) + targetVolume = produce_target_volume( + chopped, box_size, start_point, volume_resolution + ) - with open(output_filename, "wb") as f: np.save(f, targetVolume) + with open(output_filename, "wb") as f: + np.save(f, targetVolume) # stored in standard numpy pickle form + def read_target_volume(filename): - ''' + """ Takes the name of a saved target volume and loads the B vector meshgrid. Returns None if not found. - ''' + """ targetVolume = None try: @@ -195,12 +258,22 @@ def read_target_volume(filename): except: pass + ## plotting routines -def plot_fields(Bfields,box_size,start_point,vol_resolution,which_plane='z',level=0,num_contours=50): - ''' - Plots the set of Bfields in the given region, at the specified resolutions. - + +def plot_fields( + Bfields, + box_size, + start_point, + vol_resolution, + which_plane="z", + level=0, + num_contours=50, +): + """ + Plots the set of Bfields in the given region, at the specified resolutions. + Bfields: A 4D array of the Bfield. box_size: (x, y, z) dimensions of the box in cm start_point: (x, y, z) = (0, 0, 0) = bottom left corner position of the box AKA the offset @@ -208,269 +281,367 @@ def plot_fields(Bfields,box_size,start_point,vol_resolution,which_plane='z',leve which_plane: Plane to plot on, can be "x", "y" or "z" level : The "height" of the plane. For instance the Z = 5 plane would have a level of 5 num_contours: THe amount of contours on the contour plot. - - ''' + + """ # filled contour plot of Bx, By, and Bz on a chosen slice plane - X = np.linspace(start_point[0], box_size[0] + start_point[0],int(box_size[0]/vol_resolution)+1) - Y = np.linspace(start_point[1], box_size[1] + start_point[1],int(box_size[1]/vol_resolution)+1) - Z = np.linspace(start_point[2], box_size[2] + start_point[2],int(box_size[2]/vol_resolution)+1) + X = np.linspace( + start_point[0], + box_size[0] + start_point[0], + int(box_size[0] / vol_resolution) + 1, + ) + Y = np.linspace( + start_point[1], + box_size[1] + start_point[1], + int(box_size[1] / vol_resolution) + 1, + ) + Z = np.linspace( + start_point[2], + box_size[2] + start_point[2], + int(box_size[2] / vol_resolution) + 1, + ) print(Z) - if which_plane=='x': + if which_plane == "x": converted_level = np.where(X >= level) - B_sliced = [Bfields[converted_level[0][0],:,:,i].T for i in range(3)] - x_label,y_label = "y","z" - x_array,y_array = Y,Z - elif which_plane=='y': + B_sliced = [Bfields[converted_level[0][0], :, :, i].T for i in range(3)] + x_label, y_label = "y", "z" + x_array, y_array = Y, Z + elif which_plane == "y": converted_level = np.where(Y >= level) - B_sliced = [Bfields[:,converted_level[0][0],:,i].T for i in range(3)] - x_label,y_label = "x","z" - x_array,y_array = X,Z + B_sliced = [Bfields[:, converted_level[0][0], :, i].T for i in range(3)] + x_label, y_label = "x", "z" + x_array, y_array = X, Z else: converted_level = np.where(Z >= level) print(converted_level[0][0]) - B_sliced = [Bfields[:,:,converted_level[0][0],i].T for i in range(3)] - x_label,y_label = "x","y" - x_array,y_array = X,Y - - Bmin,Bmax = np.amin(B_sliced),np.amax(B_sliced) - - component_labels = ['x','y','z'] - fig,axes = plt.subplots(nrows=4,ncols=1,figsize=(10,40)) + B_sliced = [Bfields[:, :, converted_level[0][0], i].T for i in range(3)] + x_label, y_label = "x", "y" + x_array, y_array = X, Y + + Bmin, Bmax = np.amin(B_sliced), np.amax(B_sliced) + + component_labels = ["x", "y", "z"] + fig, axes = plt.subplots(nrows=4, ncols=1, figsize=(10, 40)) axes[0].set_ylabel(y_label + " (cm)") for i in range(3): - contours = axes[i].contourf(x_array,y_array,B_sliced[i], - vmin=Bmin,vmax=Bmax, - cmap=cm.magma,levels=num_contours) + contours = axes[i].contourf( + x_array, + y_array, + B_sliced[i], + vmin=Bmin, + vmax=Bmax, + cmap=cm.magma, + levels=num_contours, + ) axes[i].set_xlabel(x_label + " (cm)") - axes[i].set_title("$\mathcal{B}$"+"$_{}$".format(component_labels[i])) - + axes[i].set_title("$\mathcal{B}$" + "$_{}$".format(component_labels[i])) + axes[3].set_aspect(100) - fig.colorbar(contours,cax=axes[3],extend='both') - + fig.colorbar(contours, cax=axes[3], extend="both") + plt.tight_layout() plt.show() - + return plt + + def plot_coil(*input_filenames): - ''' + """ Plots one or more coils in space. - - input_filenames: Name of the files containing the coils. + + input_filenames: Name of the files containing the coils. Should be formatted appropriately. - ''' + """ fig = plt.figure() tick_spacing = 2 - ax = fig.add_subplot(111, projection='3d') + ax = fig.add_subplot(111, projection="3d") ax.set_xlabel("$x$ (cm)") ax.set_ylabel("$y$ (cm)") ax.set_zlabel("$z$ (cm)") for input_filename in input_filenames: coil_points = np.array(parse_coil(input_filename)) - - ax.plot3D(coil_points[0,:],coil_points[1,:],coil_points[2,:],lw=2) - for axis in [ax.xaxis,ax.yaxis,ax.zaxis]: + + ax.plot3D( + coil_points[0, :], coil_points[1, :], coil_points[2, :], lw=2, color="red" + ) + for axis in [ax.xaxis, ax.yaxis, ax.zaxis]: axis.set_major_locator(ticker.MultipleLocator(tick_spacing)) plt.tight_layout() + # set the size of the plot + fig.set_size_inches(10, 10) plt.show() - - -def create_B_x_rectangle(name,p0=[-21.59,-38.1,-21.59,1],L = 76.20,W= 43.18): - ''' + # save the plot + fig.savefig("coil.png", dpi=300) + + +def plot_coil2(coil_points): + """ + Plots one or more coils in space. + + input_filenames: Name of the files containing the coils. + Should be formatted appropriately. + """ + fig = plt.figure() + tick_spacing = 2 + ax = fig.add_subplot(111, projection="3d") + ax.set_xlabel("$x$ (cm)") + ax.set_ylabel("$y$ (cm)") + ax.set_zlabel("$z$ (cm)") + + coil_points = np.array(coil_points) + ax.plot3D( + coil_points[0, :], coil_points[1, :], coil_points[2, :], lw=2, color="red" + ) + + for axis in [ax.xaxis, ax.yaxis, ax.zaxis]: + axis.set_major_locator(ticker.MultipleLocator(tick_spacing)) + plt.tight_layout() + # set the size of the plot + fig.set_size_inches(5, 5) + plt.show() + + +def create_B_x_rectangle(name, p0=[-21.59, -38.1, -21.59, 1], L=76.20, W=43.18): + """ Creates a rectangle of the Y-Z plane that produces a B_x field. - + name: filename to output to. Should be a .txt file. p0: [x0,y0,z0,Current] Starting point of the rectangle. L: Length (on Z) W: Width (on y) - ''' - f = open(name,"w") - - p1 = [p0[0],p0[1]+W,p0[2],p0[3]] - p2 = [p0[0],p0[1]+W,p0[2]+L ,p0[3]] - p3 = [p0[0],p0[1],p0[2]+L,p0[3]] + """ + f = open(name, "w") + p1 = [p0[0], p0[1] + W, p0[2], p0[3]] + p2 = [p0[0], p0[1] + W, p0[2] + L, p0[3]] + p3 = [p0[0], p0[1], p0[2] + L, p0[3]] line = str(p0) - line = line[1:len(line)-1] + "\n" + line = line[1 : len(line) - 1] + "\n" f.write(line) line = str(p1) - line = line[1:len(line)-1] + "\n" + line = line[1 : len(line) - 1] + "\n" f.write(line) line = str(p2) - line = line[1:len(line)-1] + "\n" + line = line[1 : len(line) - 1] + "\n" f.write(line) line = str(p3) - line = line[1:len(line)-1] + "\n" + line = line[1 : len(line) - 1] + "\n" f.write(line) line = str(p0) - line = line[1:len(line)-1] + "\n" + line = line[1 : len(line) - 1] + "\n" f.write(line) f.close() -def create_B_y_rectangle(name,p0=[-21.59,-38.1,-21.59,1], L = 76.20, D= 43.18): +def create_B_y_rectangle(name, p0=[-21.59, -38.1, -21.59, 1], L=76.20, D=43.18): - ''' + """ Creates a rectangle of the X-Z plane that produces a B_y field. - + name: filename to output to. Should be a .txt file. p0: [x0,y0,z0,Current] Starting point of the rectangle. L: Length (on Z) D: Depth (on X) - ''' - f = open(name,"w") + """ + f = open(name, "w") - p1 = [p0[0], p0[1] , p0[2]+L, p0[3]] - p2 = [p0[0] + D , p0[1] , p0[2]+L, p0[3]] - p3 = [p0[0] + D , p0[1], p0[2], p0[3]] + p1 = [p0[0], p0[1], p0[2] + L, p0[3]] + p2 = [p0[0] + D, p0[1], p0[2] + L, p0[3]] + p3 = [p0[0] + D, p0[1], p0[2], p0[3]] line = str(p0) - line = line[1:len(line)-1] + "\n" + line = line[1 : len(line) - 1] + "\n" f.write(line) line = str(p1) - line = line[1:len(line)-1] + "\n" + line = line[1 : len(line) - 1] + "\n" f.write(line) line = str(p2) - line = line[1:len(line)-1] + "\n" + line = line[1 : len(line) - 1] + "\n" f.write(line) line = str(p3) - line = line[1:len(line)-1] + "\n" + line = line[1 : len(line) - 1] + "\n" f.write(line) line = str(p0) - line = line[1:len(line)-1] + "\n" + line = line[1 : len(line) - 1] + "\n" f.write(line) f.close() -def create_B_z_rectangle(name,p0=[-26.67,-26.67,-26.67,1], H= 53.340, DD = 53.340): - ''' + +def create_B_z_rectangle(name, p0=[-26.67, -26.67, -26.67, 1], H=53.340, DD=53.340): + """ Creates a rectangle of the X-Y plane that produces a B_z field. - + name: filename to output to. Should be a .txt file. p0: [x0,y0,z0,Current] Starting point of the rectangle. H: Height (on Y) DD: Depth (on X) - ''' - f = open(name,"w") + """ + f = open(name, "w") p1 = [p0[0] + DD, p0[1], p0[2], p0[3]] - p2 = [p0[0] + DD, p0[1]+H, p0[2], p0[3]] - p3 = [p0[0], p0[1]+H, p0[2], p0[3]] + p2 = [p0[0] + DD, p0[1] + H, p0[2], p0[3]] + p3 = [p0[0], p0[1] + H, p0[2], p0[3]] line = str(p0) - line = line[1:len(line)-1] + "\n" + line = line[1 : len(line) - 1] + "\n" f.write(line) line = str(p1) - line = line[1:len(line)-1] + "\n" + line = line[1 : len(line) - 1] + "\n" f.write(line) line = str(p2) - line = line[1:len(line)-1] + "\n" + line = line[1 : len(line) - 1] + "\n" f.write(line) line = str(p3) - line = line[1:len(line)-1] + "\n" + line = line[1 : len(line) - 1] + "\n" f.write(line) line = str(p0) - line = line[1:len(line)-1] + "\n" + line = line[1 : len(line) - 1] + "\n" f.write(line) f.close() - + + def helmholtz_coils(fname1, fname2, numSegments, radius, spacing, current): - ''' + """ Creates a pair of Helmholtz Coils that are parallel to the X-Y plane. - + fname1: Name of the file where the first coil will be saved. fname2: Name of the file where the second coil will be saved. numSegments: Number of segments per coil radius: Radius of the coils spacing: Spacing between the coils. The first coil will be located at -spacing/2 and the 2nd coil will be located at spacing/2 on the Z plane current: The current that goest through each coil. - ''' - f = open(fname1,"w") + """ + f = open(fname1, "w") line = "" - for i in range(0,numSegments,1): - line = str(np.cos(2*np.pi*(i)/(numSegments-1))*radius) + "," + str(np.sin(2*np.pi*(i)/(numSegments-1))*radius) + "," + str(-spacing/2.0) + "," + str(current) + "\n" + for i in range(0, numSegments, 1): + line = ( + str(np.cos(2 * np.pi * (i) / (numSegments - 1)) * radius) + + "," + + str(np.sin(2 * np.pi * (i) / (numSegments - 1)) * radius) + + "," + + str(-spacing / 2.0) + + "," + + str(current) + + "\n" + ) f.write(line) f.close() - f = open(fname2,"w") + f = open(fname2, "w") line = "" - for i in range(0,numSegments,1): - line = str(np.cos(2*np.pi*(i)/(numSegments-1))*radius) + "," + str(np.sin(2*np.pi*(i)/(numSegments-1))*radius) + "," + str(spacing/2.0) + "," + str(current) + "\n" + for i in range(0, numSegments, 1): + line = ( + str(np.cos(2 * np.pi * (i) / (numSegments - 1)) * radius) + + "," + + str(np.sin(2 * np.pi * (i) / (numSegments - 1)) * radius) + + "," + + str(spacing / 2.0) + + "," + + str(current) + + "\n" + ) f.write(line) f.close() - - + + def create_Bx_circle(fname, numSegments, radius, spacing, current, center): - ''' + """ Creates a coil on the Y-Z plane that produces a B_x field. - + fname: Name of the file where the first coil will be saved. numSegments: Number of segments per coil radius: Radius of the coil spacing: Spacing between the coil and the Y-Z plane current: The current that goest through the coil. center: (y,z) The center of the coil on the Y-Z plane - ''' - f = open(fname,"w") + """ + f = open(fname, "w") line = "" - for i in range(0,numSegments,1): - line = str(spacing) + "," + str(np.cos(2*np.pi*(i)/(numSegments-1))*radius + center[0]) + "," + str(np.sin(2*np.pi*(i)/(numSegments-1))*radius + center[1]) + "," + str(current) + "\n" + for i in range(0, numSegments, 1): + line = ( + str(spacing) + + "," + + str(np.cos(2 * np.pi * (i) / (numSegments - 1)) * radius + center[0]) + + "," + + str(np.sin(2 * np.pi * (i) / (numSegments - 1)) * radius + center[1]) + + "," + + str(current) + + "\n" + ) f.write(line) f.close() - - + + def create_By_circle(fname, numSegments, radius, spacing, current, center): - ''' + """ Creates a coil on the X-Z plane that produces a B_y field. - + fname: Name of the file where the first coil will be saved. numSegments: Number of segments per coil radius: Radius of the coil spacing: Spacing between the coil and the X-Z plane current: The current that goest through the coil. center: (x,z) The center of the coil on the X-Z plane - ''' - f = open(fname,"w") + """ + f = open(fname, "w") line = "" - for i in range(0,numSegments,1): - line = str(np.cos(2*np.pi*(i)/(numSegments-1))*radius + center[0]) + "," + str(spacing) + "," + str(np.sin(2*np.pi*(i)/(numSegments-1))*radius + center[1]) + "," + str(current) + "\n" + for i in range(0, numSegments, 1): + line = ( + str(np.cos(2 * np.pi * (i) / (numSegments - 1)) * radius + center[0]) + + "," + + str(spacing) + + "," + + str(np.sin(2 * np.pi * (i) / (numSegments - 1)) * radius + center[1]) + + "," + + str(current) + + "\n" + ) f.write(line) f.close() - - + + def create_Bz_circle(fname, numSegments, radius, spacing, current, center): - ''' + """ Creates a coil on the X-Y plane that produces a B_z field. - + fname: Name of the file where the first coil will be saved. numSegments: Number of segments per coil radius: Radius of the coil spacing: Spacing between the coil and the X-Y plane current: The current that goest through the coil. center: (x,y) The center of the coil on the X-Y plane - ''' - f = open(fname,"w") + """ + f = open(fname, "w") line = "" - for i in range(0,numSegments,1): - line = str(np.cos(2*np.pi*(i)/(numSegments-1))*radius + center[0]) + "," + str(np.sin(2*np.pi*(i)/(numSegments-1))*radius + center[1]) + "," + str(spacing) + "," + str(current) + "\n" + for i in range(0, numSegments, 1): + line = ( + str(np.cos(2 * np.pi * (i) / (numSegments - 1)) * radius + center[0]) + + "," + + str(np.sin(2 * np.pi * (i) / (numSegments - 1)) * radius + center[1]) + + "," + + str(spacing) + + "," + + str(current) + + "\n" + ) f.write(line) f.close() - diff --git a/simulations/coil_6_spiral_coil.svg b/simulations/coil_6_spiral_coil.svg new file mode 100644 index 0000000..5df132c --- /dev/null +++ b/simulations/coil_6_spiral_coil.svg @@ -0,0 +1 @@ + \ No newline at end of file diff --git a/simulations/compare_coil_fields.ipynb b/simulations/compare_coil_fields.ipynb new file mode 100644 index 0000000..8a095df --- /dev/null +++ b/simulations/compare_coil_fields.ipynb @@ -0,0 +1,359 @@ +{ + "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": [ + "# set the size larger\n", + "plt.rcParams[\"figure.figsize\"] = (10, 10)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# COIL_NAME = \"coil_6_spiral\"\n", + "# COIL_NAME = \"coil_6_custom\"\n", + "# COIL_NAME = \"coil_12_spiral\"\n", + "COIL_NAME = \"coil_12_custom\"\n", + "COIL = \"coils/\" + COIL_NAME + \".csv\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "bs.plot_coil(COIL)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "bs.write_target_volume(\n", + " COIL,\n", + " COIL_NAME + \"volume.vol\",\n", + " (4, 4, 4),\n", + " (-2, -2, -2),\n", + " 0.1,\n", + " 0.1,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# read in the results\n", + "total_volume = bs.read_target_volume(COIL_NAME + \"volume.vol\")\n", + "\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)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import plotly.graph_objects as go\n", + "\n", + "X, Y, Z = np.mgrid[-2:2:41j, -2:2:41j, -2:2:41j]\n", + "# values = np.sqrt(total_volume[:,:,:,0]*total_volume[:,:,:,0] +\n", + "# total_volume[:,:,:,1]*total_volume[:,:,:,1] +\n", + "# total_volume[:,:,:,2]*total_volume[:,:,:,2])\n", + "values = total_volume[:, :, :, 2]\n", + "fig = go.Figure(\n", + " data=go.Volume(\n", + " x=X.flatten(),\n", + " y=Y.flatten(),\n", + " z=Z.flatten(),\n", + " value=np.abs(values.flatten()),\n", + " isomin=0.1,\n", + " isomax=0.5,\n", + " opacity=0.1, # needs to be small to see through all surfaces\n", + " surface_count=17, # needs to be a large number for good volume rendering\n", + " )\n", + ")\n", + "\n", + "coil = bs.parse_coil(\"coils/\" + COIL_NAME + \".csv\")\n", + "coil = coil[0:3, :]\n", + "# add the coil as a line plot\n", + "fig.add_trace(\n", + " go.Scatter3d(\n", + " x=coil[0, :],\n", + " y=coil[1, :],\n", + " z=coil[2, :],\n", + " mode=\"lines\",\n", + " line=dict(color=\"red\", width=5),\n", + " )\n", + ")\n", + "# make the plotly graph bigger\n", + "fig.update_layout(\n", + " autosize=False,\n", + " width=1400,\n", + " height=1000,\n", + ")\n", + "fig.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "slice = [total_volume[:, 21, :, 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, x_contour.shape[0], 1)\n", + "y = np.arange(0, z_contour.shape[1], 1)\n", + "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", + "\n", + "# color = 2 * np.log(np.hypot(x_contour, z_contour))\n", + "color = np.log(mag)\n", + "stream_plot = plt.streamplot(\n", + " x,\n", + " y,\n", + " x_contour,\n", + " z_contour,\n", + " color=color,\n", + " linewidth=1,\n", + " cmap=plt.cm.Reds,\n", + " density=2,\n", + " arrowstyle=\"->\",\n", + " arrowsize=1.5,\n", + ")\n", + "# set the axis to be equal\n", + "plt.axis(\"equal\")\n", + "\n", + "# set the plt size to 5x5\n", + "plt.rcParams[\"figure.figsize\"] = (10, 10)\n", + "\n", + "# save the figure\n", + "plt.savefig(COIL_NAME + \"_field_lines.png\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.cm as cm\n", + "\n", + "# get the total left and right field in the X direction\n", + "slice = [total_volume[:, :, 21, i].T for i in range(3)]\n", + "x_contour = slice[0]\n", + "y_contour = slice[1]\n", + "z_contour = slice[2]\n", + "\n", + "# vmin = np.min([np.min(x_contour), np.min(z_contour), np.min(y_contour)])\n", + "# vmax = np.max([np.max(x_contour), np.max(z_contour), np.max(y_contour)])\n", + "\n", + "vmin = np.min(x_contour)\n", + "vmax = np.max(x_contour)\n", + "\n", + "# -7.970943527635176 7.775784181783158\n", + "vmin = -7.970943527635176\n", + "vmax = 7.775784181783158\n", + "\n", + "print(vmin, vmax)\n", + "\n", + "left_x = np.sum(x_contour[x_contour < 0])\n", + "right_x = np.sum(x_contour[x_contour > 0])\n", + "\n", + "print(\"Left X: \", left_x)\n", + "print(\"Right X: \", right_x)\n", + "\n", + "plt.rcParams[\"figure.figsize\"] = (5, 5)\n", + "\n", + "# make the plot larger\n", + "plt.rcParams[\"figure.figsize\"] = (10, 10)\n", + "\n", + "# plt.imshow(x_contour, cmap='hot', interpolation='nearest')\n", + "plt.contourf(\n", + " x_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", + "\n", + "print(vmin, vmax)\n", + "\n", + "# 6 coil custom\n", + "# Left X: -565.6790142743889\n", + "# Right X: 544.1011578804321\n", + "\n", + "# 6 coil spiral\n", + "# Left X: -512.1107132193335\n", + "# Right X: 534.1946450523811\n", + "\n", + "# 12 coil spiral\n", + "# Left X: -218.54696343677543\n", + "# Right X: 233.84835661493196\n", + "\n", + "# 12 coil custom\n", + "# Left X: -361.95894587646944\n", + "# Right X: 347.13174239269904" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# plot the magitude of the slice\n", + "plt.imshow(np.log(mag), cmap=\"hot\", interpolation=\"nearest\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "slice = [total_volume[:, :, 25, i].T for i in range(3)]\n", + "x_contour = slice[0]\n", + "y_contour = slice[1]\n", + "z_contour = slice[2]\n", + "\n", + "mag = np.sqrt(x_contour**2 + y_contour**2 + z_contour**2)\n", + "\n", + "print(np.max(mag), np.min(mag))\n", + "\n", + "# plot a histogram of the magnitude\n", + "# plt.hist(mag.flatten(), bins=10)\n", + "\n", + "plt.hist(z_contour.flatten(), bins=10)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# show the z contour as a 3d surface plot\n", + "fig = plt.figure()\n", + "\n", + "# change the figure size to 16.9\n", + "ax = fig.add_subplot(111, projection=\"3d\")\n", + "\n", + "xG_ = xG[10:30, 10:30]\n", + "yG_ = yG[10:30, 10:30]\n", + "z_contour_ = z_contour[10:30, 10:30]\n", + "\n", + "print(xG.shape, yG.shape, z_contour.shape)\n", + "\n", + "# change the axis limits to 10 to 30\n", + "ax.set_xlim(10, 30)\n", + "ax.set_ylim(10, 30)\n", + "surf = ax.plot_surface(\n", + " xG_, yG_, z_contour_, cmap=cm.plasma, linewidth=0, antialiased=False, alpha=0.5\n", + ")\n", + "# add the coil to the plot as a set of lines - the X and Y coordinates of the coil are in the first two elements of the coil array\n", + "coil = bs.parse_coil(\"coils/\" + COIL_NAME + \".csv\")\n", + "coil = coil[0:2, :]\n", + "# need to scale the coil by 10 and shift it by 20, 20\n", + "coil = coil * 10 + 20\n", + "# draw the coil in yellow\n", + "ax.plot(coil[0, :], coil[1, :], 0, color=\"black\")\n", + "\n", + "# hide the axes\n", + "ax.set_axis_off()\n", + "\n", + "# save to png\n", + "plt.savefig(COIL_NAME + \"_field_surface_z.png\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# write out the coil as SVG\n", + "coil = bs.parse_coil(\"coils/\" + COIL_NAME + \".csv\")\n", + "coil = coil[0:2, :]\n", + "with open(COIL_NAME + \"_coil.svg\", \"w\") as f:\n", + " # write the header\n", + " f.write('')\n", + " # the X and Y points of the coil the X and Y coordinates of the coil are in the first two elements of the coil array\n", + " # write out the coil as a set of lines\n", + " print(len(coil[0]) - 1)\n", + " # write out the coil as a polyline\n", + " f.write('')\n", + " # write the footer\n", + " f.write(\"\")" + ] + } + ], + "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/helpers.py b/simulations/helpers.py new file mode 100644 index 0000000..234576b --- /dev/null +++ b/simulations/helpers.py @@ -0,0 +1,111 @@ +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 new file mode 100644 index 0000000..c10f20e --- /dev/null +++ b/simulations/magnet_field.ipynb @@ -0,0 +1,607 @@ +{ + "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 new file mode 100644 index 0000000..a41c36c --- /dev/null +++ b/simulations/magnet_field_single_coil.ipynb @@ -0,0 +1,525 @@ +{ + "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/simulate_coil.ipynb b/simulations/simulate_coil.ipynb index d11dc2d..31ffbe3 100644 --- a/simulations/simulate_coil.ipynb +++ b/simulations/simulate_coil.ipynb @@ -7,7 +7,8 @@ "outputs": [], "source": [ "import numpy as np\n", - "import biot_savart_v4_3 as bs" + "import biot_savart_v4_3 as bs\n", + "import matplotlib.pyplot as plt" ] }, { @@ -16,7 +17,7 @@ "metadata": {}, "outputs": [], "source": [ - "bs.plot_coil(\"coil-4-top-1.txt\")" + "bs.plot_coil(\"coil-4-top-1.csv\")" ] }, { @@ -25,7 +26,7 @@ "metadata": {}, "outputs": [], "source": [ - "bs.plot_coil(\"coil-4-bottom-1.txt\")" + "bs.plot_coil(\"coil-4-bottom-1.csv\")" ] }, { @@ -36,7 +37,7 @@ "source": [ "for i in range(0, 12):\n", " bs.write_target_volume(\n", - " f\"coil-4-top-{i}.txt\",\n", + " f\"coil-4-top-{i}.csv\",\n", " f\"targetvol-top-{i}.vol\",\n", " (6, 6, 1),\n", " (-3, -3, -0.5),\n", @@ -44,7 +45,7 @@ " 0.1,\n", " )\n", " bs.write_target_volume(\n", - " f\"coil-4-bottom-{i}.txt\",\n", + " f\"coil-4-bottom-{i}.csv\",\n", " f\"targetvol-bottom-{i}.vol\",\n", " (6, 6, 1),\n", " (-3, -3, -0.5),\n", @@ -74,6 +75,62 @@ "# 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": {