diff --git a/Radia_Ex06.ipynb b/Radia_Ex06.ipynb index a353d53..48e9997 100644 --- a/Radia_Ex06.ipynb +++ b/Radia_Ex06.ipynb @@ -316,14 +316,22 @@ " return s/np/(ro**q) " ] }, + { + "cell_type": "markdown", + "metadata": { + "jp-MarkdownHeadingCollapsed": true, + "tags": [] + }, + "source": [ + "## _Define a general function to build a quadrupole focusing magnet_" + ] + }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ - "## _Define a general function to build a quadrupole focusing magnet_\n", - "\n", "Here we define a function that creates the quadrupole focusing magnet\n", "of this example. The various arguments (detailed in the function’s\n", "docstring) specify the geometry, material properties, current, and \n", @@ -423,6 +431,7 @@ { "cell_type": "markdown", "metadata": { + "jp-MarkdownHeadingCollapsed": true, "tags": [] }, "source": [ @@ -557,11 +566,20 @@ { "cell_type": "markdown", "metadata": { + "jp-MarkdownHeadingCollapsed": true, + "tags": [] + }, + "source": [ + "## _Plots of the magnetic field and field integrals_" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "jp-MarkdownHeadingCollapsed": true, "tags": [] }, "source": [ - "## _Plots of the magnetic field and field integrals_\n", - "\n", "Here we show plots of magnetic field in the gap and corresponding\n", "field integrals. The field values are obtained by calling `Fld` on\n", "a list of points. One may also use `FldLst`.\n", @@ -698,78 +716,92 @@ " y_hat = [0, 1, 0]\n", " z_hat = [0, 0, 1]\n", "\n", - " # discretize hyperbola that defines the pole tip\n", + " # default segmentation parameters\n", + " nx = 2\n", + " ny = 2\n", + " n1 = [nx,ny,3]\n", + " n2 = [nx,ny,3]\n", + " np3 = 2\n", + " nr3 = ny\n", + " n4 = [nx,3,ny]\n", + " np5 = ceil(np3/2)\n", + " nr5 = ny\n", + "\n", + " # default colors\n", + " ironcolor = [0.0, 0.5, 0.8]\n", + " coilcolor = [1.0, 0.2, 0.0]\n", + "\n", + " # specify the particular iron to use\n", + " ironmat = rad.MatSatIsoFrm([2000, 2], [0.1, 2], [0.1, 2])\n", + "\n", + " # discretize path that defines the pole tip\n", " # :: z^2 - (hyp * y)^2 = z0^2, w/ asymptotes z = ±hyp * y\n", - " tan_np = tan(pi / n_poles)\n", + " tan_np = tan(π / n_poles)\n", " hyp = 1 / tan_np # slope of hyperbola's asymptote\n", " z0 = gap / 2\n", " ym = width / 2\n", - " a_max = asinh(hyp * ym / z0)\n", - " zm = z0 * cosh(a_max)\n", + " zm = hypot(hyp * ym, z0)\n", " # sanity check: pole tip includes all of pole face?\n", " assert zm <= z0 + height, \\\n", " \"Pole tip height too short to accommodate entire curved pole face.\"\n", - " da = a_max / n_curve\n", - " na = n_curve + 3 # go two points beyond so we don't have to extend the array\n", - " qq = [ [(z0 / hyp) * sinh(ia * da), z0 * cosh(ia * da)] for ia in range(na) ]\n", + " # construct hyperbolic path\n", + " dy = ym / n_curve\n", + " ny = n_curve + 3 # go two points beyond so we don't have to extend the array\n", + " Γ_tip = [ [iy * dy, hypot(hyp * iy * dy, z0)] for iy in range(ny) ]\n", " # and\n", " # modify last two points so as to outline lower portion of the (half) pole tip\n", " ht_lower = poletip_frac * height\n", " # sanity check: lower fraction of pole tip includes all of pole face?\n", " assert zm <= z0 + ht_lower, \\\n", " \"Lower fraction of pole tip cannot accommodate entire pole face.\"\n", - " qq[n_curve + 1] = [ym, z0 + ht_lower]\n", - " qq[n_curve + 2] = [ 0, z0 + ht_lower]\n", + " Γ_tip[n_curve + 1] = [ym, z0 + ht_lower]\n", + " Γ_tip[n_curve + 2] = [ 0, z0 + ht_lower]\n", "\n", " # create and segment the lower portion of the (half) pole tip\n", - " g1 = rad.ObjThckPgn(thick / 4, thick / 2, qq)\n", - " rad.ObjDivMag(g1, n1)\n", + " g_tip = rad.ObjThckPgn(thick / 4, thick / 2, Γ_tip)\n", + " rad.ObjDivMag(g_tip, n1)\n", "\n", " # create and segment the upper portion of the (half) pole tip\n", " ht_upper = height - ht_lower\n", - " g2 = rad.ObjRecMag([thick / 4, width / 4, z0 + height - ht_upper / 2],\n", + " g_pole = rad.ObjRecMag([thick / 4, width / 4, z0 + height - ht_upper / 2],\n", " [thick / 2, width / 2, ht_upper])\n", - " rad.ObjDivMag(g2, n2)\n", + " rad.ObjDivMag(g_pole, n2)\n", "\n", " # combine the lower and upper portions of the (half) pole tip\n", - " gg = rad.ObjCnt([g1, g2])\n", + " g_pt = rad.ObjCnt([g_tip, g_pole])\n", " # and\n", " # cut chamfer, then retain desired metal\n", - " gp = rad.ObjCutMag(gg, [thick / 2 - chamfer, 0, z0], [1, 0, -1])[0]\n", + " g_poletip = rad.ObjCutMag(g_pt, [thick / 2 - chamfer, 0, z0], [1, 0, -1])[0]\n", "\n", " # create and segment \"corner\" above (half) pole tip\n", " depth = yoketip_frac * width\n", - " g3 = rad.ObjRecMag([thick / 4, width / 4, z0 + height + depth / 2],\n", + " g_top = rad.ObjRecMag([thick / 4, width / 4, z0 + height + depth / 2],\n", " [thick / 2, width / 2, depth])\n", " cy = [[[0, width / 2, gap / 2 + height], [1, 0, 0]], [0, 0, gap / 2 + height], 2 * depth / width]\n", - " rad.ObjDivMag(g3, [nr3, np3, nx], 'cyl', cy)\n", + " rad.ObjDivMag(g_top, [nr3, np3, nx], 'cyl', cy)\n", "\n", " # create and segment horizontal yoke segment to corner\n", " length = tan_np * (gap / 2 + height) - width / 2\n", - " g4 = rad.ObjRecMag([thick / 4, width / 2 + length / 2, gap / 2 + height + depth / 2],\n", + " g_bar = rad.ObjRecMag([thick / 4, width / 2 + length / 2, gap / 2 + height + depth / 2],\n", " [thick / 2,length, depth])\n", - " rad.ObjDivMag(g4, n4)\n", + " rad.ObjDivMag(g_bar, n4)\n", "\n", " # outline the corner\n", " yc = width / 2 + length\n", " zc = gap / 2 + height\n", - " q_corner = [[yc, zc], [yc, zc + depth], [yc + depth * tan_np, zc + depth]]\n", + " Γ_corner = [[yc, zc], [yc, zc + depth], [yc + depth * tan_np, zc + depth]]\n", " # and\n", " # create and segment yoke corner\n", - " g5 = rad.ObjThckPgn(thick / 4, thick / 2, q_corner)\n", + " g_corner = rad.ObjThckPgn(thick / 4, thick / 2, Γ_corner)\n", " cy = [[[0, yc, zc], [1, 0, 0]], [0, yc, zc + depth], 1]\n", - " rad.ObjDivMag(g5, [nr5, np5, nx], 'cyl', cy)\n", + " rad.ObjDivMag(g_corner, [nr5, np5, nx], 'cyl', cy)\n", "\n", - " # create container for all the metal\n", - " g = rad.ObjCnt([gp, g3, g4, g5])\n", - " gd = rad.ObjCnt([g])\n", + " # create container for the (half) pole tip plus attached crossbar\n", + " g_yoke = rad.ObjCnt([g_poletip, g_top, g_bar, g_corner])\n", + " # specify the iron\n", + " rad.MatApl(g_yoke, ironmat)\n", " # and set color for iron\n", - " rad.ObjDrwAtr(g, ironcolor)\n", - "\n", - " # sanity check: coil_1 does not cross sector diagonal\n", - " # assert y_coil < y1\n", - " # sanity check: coil_2 does not cross sector diagonal\n", - " # assert y_coil < y2\n", + " rad.ObjDrwAtr(g_yoke, ironcolor)\n", "\n", " # create coil1\n", " ht_coil = height - tip_coil_sep\n", @@ -781,6 +813,8 @@ " coil1 = rad.ObjRaceTrk([0, 0, z0 + height - ht_coil / 2], [r_min, r_max],\n", " [thick, width - 2 * r_min],\n", " ht_coil, 3, curr_density)\n", + " # and set color for coil1\n", + " rad.ObjDrwAtr(coil1, coilcolor)\n", "\n", " # create coil2\n", " ht_coil = (height - tip_coil_sep) / 2\n", @@ -789,33 +823,31 @@ " coil2 = rad.ObjRaceTrk([0, 0, z0 + height - ht_coil / 2], [r_max, r2],\n", " [thick, width - 2 * r_min],\n", " ht_coil, 3, curr_density)\n", - "\n", - " # set color for coils\n", - " rad.ObjDrwAtr(coil1, coilcolor)\n", + " # and set color for coil2\n", " rad.ObjDrwAtr(coil2, coilcolor)\n", "\n", - " # apply symmetries to create full pole tip plus attached yoke\n", - " ctr = [0, 0, 0]\n", - " x_hat = [1, 0, 0]\n", - " # reflection in y-z plane, with zero field perpendicular to the plane\n", - " rad.TrfZerPerp(gd, ctr, x_hat)\n", - " # reflection in z-x plane, with zero field perpendicular to the plane\n", - " rad.TrfZerPerp(gd, ctr, y_hat)\n", + " # apply symmetries to create full pole tip plus attached crossbar\n", + " # :: reflection in y-z plane, with zero field perpendicular to the plane\n", + " rad.TrfZerPerp(g_yoke, ctr, x_hat)\n", + " # :: reflection in z-x plane, with zero field perpendicular to the plane\n", + " rad.TrfZerPerp(g_yoke, ctr, y_hat)\n", "\n", - " # create container for full pole tip with attached yoke and coils\n", - " t = rad.ObjCnt([gd, coil1, coil2])\n", + " # create container for full magnet: here iron yoke plus coils in one sector\n", + " g_magnet = rad.ObjCnt([g_yoke, coil1, coil2])\n", "\n", - " # reflection across diagonal plane, with zero field parallel to the plane\n", - " rad.TrfZerPara(t, ctr, [0, cos(pi / n_poles), sin(pi / n_poles)])\n", + " # :: reflection across diagonal plane, with zero field parallel to the plane\n", + " rad.TrfZerPara(g_magnet, ctr, [0, cos(pi / n_poles), sin(pi / n_poles)])\n", " # ==>> at this point we have a matched pair of pole tips\n", " # they subtend an angle 2 * (2π / n_poles) = 4π / n_poles\n", "\n", - " # apply symmetries to create full multipole electromagnet\n", - " rad.TrfMlt(t, rad.TrfRot(ctr, x_hat, 4 * pi / n_poles), int(n_poles / 2))\n", - " rad.MatApl(g, ironmat)\n", - " rad.TrfOrnt(t, rad.TrfRot(ctr, x_hat, pi / n_poles))\n", + " # apply rotation symmetries to create full multipole electromagnet\n", + " rad.TrfMlt(g_magnet, rad.TrfRot(ctr, x_hat, 4 * pi / n_poles), int(n_poles / 2))\n", "\n", - " return t" + " # ensure upright orientation of this multipole\n", + " if n_poles % 4 == 0:\n", + " rad.TrfOrnt(g_magnet, rad.TrfRot(ctr, x_hat, π / n_poles))\n", + "\n", + " return g_magnet" ] }, { @@ -832,22 +864,7 @@ "height = 50. # height of pole tip / mm\n", "chamfer = 8. # size of chamfer\n", "tip_coil_sep = 10.\n", - "curr_density = -3. # current density / (A / mm^2)\n", - "\n", - "# segmentation parameters\n", - "nx = 2\n", - "ny = 2\n", - "n1 = [nx,ny,3]\n", - "n2 = [nx,ny,3]\n", - "np3 = 2\n", - "nr3 = ny\n", - "n4 = [nx,3,ny]\n", - "np5 = ceil(np3/2)\n", - "nr5 = ny\n", - "\n", - "# colors\n", - "ironcolor = [0.0, 1.0, 1.0]\n", - "coilcolor = [1.0, 0.0, 0.0]" + "curr_density = -3. # current density / (A / mm^2)" ] }, { @@ -857,7 +874,7 @@ "outputs": [], "source": [ "rad.UtiDelAll()\n", - "ironmat = rad.MatSatIsoFrm([2000,2],[0.1,2],[0.1,2])\n", + "# ironmat = rad.MatSatIsoFrm([2000,2], [0.1,2], [0.1,2])\n", "t0 = tm.time()\n", "quad = create_multipole(n_poles, thick, width, gap, height, chamfer, tip_coil_sep, curr_density)\n", "t1 = tm.time()\n", @@ -988,6 +1005,16 @@ "uti_plot_show()" ] }, + { + "cell_type": "raw", + "metadata": {}, + "source": [ + "numpy.savetxt(nb_dir + \"BzVy1-new.txt\", BzVy1)\n", + "numpy.savetxt(nb_dir + \"BzVy2-new.txt\", BzVy2)\n", + "numpy.savetxt(nb_dir + \"BzVx1-new.txt\", BzVx1)\n", + "numpy.savetxt(nb_dir + \"BzVx2-new.txt\", BzVx2)" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -1001,10 +1028,10 @@ "metadata": {}, "outputs": [], "source": [ - "BzVy1_sv = numpy.loadtxt(nb_dir + \"BzVy1.txt\")\n", - "BzVy2_sv = numpy.loadtxt(nb_dir + \"BzVy2.txt\")\n", - "BzVx1_sv = numpy.loadtxt(nb_dir + \"BzVx1.txt\")\n", - "BzVx2_sv = numpy.loadtxt(nb_dir + \"BzVx2.txt\")" + "BzVy1_sv = numpy.loadtxt(nb_dir + \"BzVy1-new.txt\")\n", + "BzVy2_sv = numpy.loadtxt(nb_dir + \"BzVy2-new.txt\")\n", + "BzVx1_sv = numpy.loadtxt(nb_dir + \"BzVx1-new.txt\")\n", + "BzVx2_sv = numpy.loadtxt(nb_dir + \"BzVx2-new.txt\")" ] }, {