Skip to content

Commit

Permalink
Make variable names more readble.
Browse files Browse the repository at this point in the history
  • Loading branch information
@dtabell
dtabell committed Aug 16, 2023
1 parent 8ac129c commit 4280a2d
Showing 1 changed file with 102 additions and 75 deletions.
177 changes: 102 additions & 75 deletions Radia_Ex06.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down Expand Up @@ -423,6 +431,7 @@
{
"cell_type": "markdown",
"metadata": {
"jp-MarkdownHeadingCollapsed": true,
"tags": []
},
"source": [
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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",
Expand All @@ -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",
Expand All @@ -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"
]
},
{
Expand All @@ -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)"
]
},
{
Expand All @@ -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",
Expand Down Expand Up @@ -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": {},
Expand All @@ -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\")"
]
},
{
Expand Down

0 comments on commit 4280a2d

Please sign in to comment.