Radia_Ex03.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# RADIA Example 3: Simple U46 Hybrid Undulator\n",
    "\n",
    "This example creates and solves a simple U46 Hybrid Undulator made with\n",
    "rectangular blocks.\n",
    "\n",
    "To obtain an acceptable precision for the magnetic field, it is essential to\n",
    "increase the segmentation, which, of course, requires more memory and CPU time.\n",
    "The initial segmentation in this example has been set to `nsp = [2,2,5]` for\n",
    "the pole and `nsm = [1,3,1]` for the magnet. These settings yield an absolute\n",
    "error of order 1% in the peak field of the undulator. The user can vary both\n",
    "`nsp` and `nsm` (in the section **_Create the undulator magnet_**) and assess the consequences\n",
    "for memory usage, CPU time, and convergence of the peak field value. As a rule,\n",
    "the iron will require finer segmentation because its magnetization is more\n",
    "sensitive to the external field.\n",
    "\n",
    "One can also—for undulators such as this one—compute field integrals and forces\n",
    "on the magnets and poles. In addition, one can build a wedge-pole undulator by\n",
    "using polyhedrons or extruded polygons instead of rectangular blocks. To learn\n",
    "more, read the [Radia documentation](\n",
    "  https://www.esrf.eu/Accelerators/Groups/InsertionDevices/Software/Radia/Documentation\n",
    "  \"RADIA documentation at ESRF\"),\n",
    "experiment with the code, and develop experience.\n",
    "\n",
    "## _Import Radia and other packages_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import scipy.constants as sc\n",
    "import matplotlib.pyplot as plt\n",
    "import time as tm\n",
    "\n",
    "import radia as rad\n",
    "import ipywidgets\n",
    "from jupyter_rs_radia import radia_viewer\n",
    "\n",
    "# import an illustration of this magnet\n",
    "from IPython.display import Image\n",
    "Image(filename=('./U46HybridUnd.png'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## _Define function to build the undulator_\n",
    "\n",
    "We here define a function that constructs the desired undulator geometry in Radia.\n",
    "It constructs in Radia an array of rectangular permanent magnets and iron blocks,\n",
    "each with given dimensions, magnetization, segmentation, material properties,\n",
    "and colors (for plotting). It then constructs their mirror symmetric counterparts.\n",
    "The function is generic and can be used to build almost any hybrid undulator—though\n",
    "without optimizing the design of terminations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def undulator(lp, mp, np, cp,  lm, mm, nm, cm,  gap, gapOffset, numPer):\n",
    "    \"\"\"\n",
    "    create hybrid undulator magnet\n",
    "    arguments:\n",
    "      lp = [lpx, lpy, lpz] = dimensions of the iron poles / mm\n",
    "      mp = magnetic properties of the iron poles (M-H curve)\n",
    "      np = segmentation of the iron poles\n",
    "      cp = [r,g,b] = color for the iron poles\n",
    "      lm = [lmx, lmy, lmz] = dimensions of the magnet blocks / mm\n",
    "      mm = magnetic properties of the magnet blocks (remanent magnetization)\n",
    "      nm = segmentation of the magnet blocks\n",
    "      cm = [r,g,b] = color for the magnet blocks\n",
    "      gap = undulator gap / mm\n",
    "      gapOffset = vertical offset / mm of the magnet blocks w/rt the poles\n",
    "      numPer = number of full periods of the undulator magnetic field\n",
    "    return: Radia representations of\n",
    "      undulator group, poles, permanent magnets\n",
    "    \"\"\"\n",
    "    zer = [0,0,0]\n",
    "\n",
    "    # full magnet will be assembled into this Radia group\n",
    "    grp = rad.ObjCnt([])\n",
    "\n",
    "    # principal poles and magnet blocks in octant(+,+,–)\n",
    "    # -- half pole\n",
    "    y = lp[1]/4\n",
    "    pole = rad.ObjFullMag([lp[0]/4,    y,    -lp[2]/2 - gap/2],\n",
    "                          [lp[0]/2, lp[1]/2,  lp[2]          ], zer, np, grp, mp, cp)\n",
    "    y += lp[1]/4\n",
    "    # -- magnet and pole pairs\n",
    "    mDir = -1\n",
    "    for i in range(0, numPer):\n",
    "        initM = [0, mDir, 0]; mDir *= -1\n",
    "        y += lm[1]/2\n",
    "        magnet = rad.ObjFullMag([lm[0]/4,   y,   -lm[2]/2 - gap/2 - gapOffset],\n",
    "                                [lm[0]/2, lm[1],  lm[2]                      ], initM, nm, grp, mm, cm)\n",
    "        y += (lm[1] + lp[1])/2\n",
    "        pole = rad.ObjFullMag([lp[0]/4,   y,   -lp[2]/2 - gap/2],\n",
    "                              [lp[0]/2, lp[1],  lp[2]          ], zer, np, grp, mp, cp)\n",
    "        y += lp[1]/2\n",
    "    # -- end magnet block\n",
    "    initM = [0, mDir, 0]\n",
    "    y += lm[1]/4;\n",
    "    magnet = rad.ObjFullMag([lm[0]/4,   y,     -lm[2]/2 - gap/2 - gapOffset],\n",
    "                            [lm[0]/2, lm[1]/2,  lm[2]                      ], initM, nm, grp, mm, cm)\n",
    "\n",
    "    # use mirror symmetry to define the full undulator\n",
    "    rad.TrfZerPerp(grp, zer, [1,0,0]) # reflect in the (y,z) plane\n",
    "    rad.TrfZerPara(grp, zer, [0,0,1]) # reflect in the (x,y) plane\n",
    "    rad.TrfZerPerp(grp, zer, [0,1,0]) # reflect in the (z,x) plane\n",
    "    \n",
    "    return grp, pole, magnet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this version supplies a work-around for the viz bug\n",
    "def undulatorX(lp, mp, np, cp,  lm, mm, nm, cm,  gap, gapOffset, numPer):\n",
    "    \"\"\"\n",
    "    create hybrid undulator magnet\n",
    "    arguments:\n",
    "      lp = [lpx, lpy, lpz] = dimensions of the iron poles / mm\n",
    "      mp = magnetic properties of the iron poles (M-H curve)\n",
    "      np = segmentation of the iron poles\n",
    "      cp = [r,g,b] = color for the iron poles\n",
    "      lm = [lmx, lmy, lmz] = dimensions of the magnet blocks / mm\n",
    "      mm = magnetic properties of the magnet blocks (remanent magnetization)\n",
    "      nm = segmentation of the magnet blocks\n",
    "      cm = [r,g,b] = color for the magnet blocks\n",
    "      gap = undulator gap / mm\n",
    "      gapOffset = vertical offset / mm of the magnet blocks w/rt the poles\n",
    "      numPer = number of full periods of the undulator magnetic field\n",
    "    return: Radia representations of\n",
    "      undulator group, poles, permanent magnets\n",
    "    \"\"\"\n",
    "    zer = [0,0,0]\n",
    "\n",
    "    # full magnet will be assembled in the following Radia groups\n",
    "    grp_m = rad.ObjCnt([])            # iron magnet\n",
    "    grp_p = rad.ObjCnt([])            # permanent magnet poles\n",
    "    grp = rad.ObjCnt([grp_m, grp_p])  # full magnet\n",
    "\n",
    "    # principal poles and magnet blocks in octant(+,+,–)\n",
    "    # -- half pole\n",
    "    y = lp[1]/4\n",
    "    pole = rad.ObjFullMag([lp[0]/4,    y,    -lp[2]/2 - gap/2],\n",
    "                          [lp[0]/2, lp[1]/2,  lp[2]          ], zer, np, grp_p, mp, cp)\n",
    "    y += lp[1]/4\n",
    "    # -- magnet and pole pairs\n",
    "    mDir = -1\n",
    "    for i in range(0, numPer):\n",
    "        initM = [0, mDir, 0]; mDir *= -1\n",
    "        y += lm[1]/2\n",
    "        magnet = rad.ObjFullMag([lm[0]/4,   y,   -lm[2]/2 - gap/2 - gapOffset],\n",
    "                                [lm[0]/2, lm[1],  lm[2]                      ], initM, nm, grp_m, mm, cm)\n",
    "        y += (lm[1] + lp[1])/2\n",
    "        pole = rad.ObjFullMag([lp[0]/4,   y,   -lp[2]/2 - gap/2],\n",
    "                              [lp[0]/2, lp[1],  lp[2]          ], zer, np, grp_p, mp, cp)\n",
    "        y += lp[1]/2\n",
    "    # -- end magnet block\n",
    "    initM = [0, mDir, 0]\n",
    "    y += lm[1]/4;\n",
    "    magnet = rad.ObjFullMag([lm[0]/4,   y,     -lm[2]/2 - gap/2 - gapOffset],\n",
    "                            [lm[0]/2, lm[1]/2,  lm[2]                      ], initM, nm, grp_m, mm, cm)\n",
    "\n",
    "    # use mirror symmetry to define the full undulator\n",
    "    rad.TrfZerPerp(grp_m, zer, [1,0,0]) # reflect in the (y,z) plane\n",
    "    rad.TrfZerPara(grp_m, zer, [0,0,1]) # reflect in the (x,y) plane\n",
    "    rad.TrfZerPerp(grp_m, zer, [0,1,0]) # reflect in the (z,x) plane\n",
    "    \n",
    "    # use mirror symmetry to define the full undulator\n",
    "    rad.TrfZerPerp(grp_p, zer, [1,0,0]) # reflect in the (y,z) plane\n",
    "    rad.TrfZerPara(grp_p, zer, [0,0,1]) # reflect in the (x,y) plane\n",
    "    rad.TrfZerPerp(grp_p, zer, [0,1,0]) # reflect in the (z,x) plane\n",
    "    \n",
    "    return grp, pole, magnet"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## _Define functions for the magnetic materials and fields_\n",
    "\n",
    "The first of the following functions returns a Radia representation of the magnetization curves _M_ (in Tesla) vs _H_ (in Amp/m) for both the pole tips and the magnets. The material for the poles is defined by calling the function `rad.MatSatIsoTab` with the data table `ma`, which contains the paired values of _H_ and _M_, both expressed in Tesla. (The conversion of _H_ in Amp/m to Tesla is accomplished simply by multiplying _H_ by _μ<sub>o</sub>_, i.e. `4*pi/10^7`.) The material for the permanent magnets is defined by calling `rad.MatStd`, with arguments _material type_ and _remanent magnetization_. We will later plot the resulting _M_(_H_) curves for verification.\n",
    "\n",
    "The next function extracts magnetic field information from a given magnetic material defined by Radia. This allows us to plot the _M_(_H_) curves tocompare with our input data. And the third function calculates (for a _solved_ Radia magnet) the vertical field along the Y (longitudinal) axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# magnetic constant (vacuum magnetic permeability)\n",
    "mu0 = 4 * np.pi / 1e7\n",
    "\n",
    "def materials(H, M,  smat, rm):\n",
    "    \"\"\"\n",
    "    define magnetic materials for the undulator poles and magnets\n",
    "    arguments:\n",
    "      H    = list of magnetic field values / (Amp/m)\n",
    "      M    = corresponding magnetization values / T\n",
    "      smat = material type string\n",
    "      rm   = remanent magnetization / T\n",
    "    return: Radia representations of ...\n",
    "      pole-tip material, magnet material\n",
    "    \"\"\"\n",
    "    # -- magnetic property of poles\n",
    "    #convH = 4. * np.pi / 1.0e7\n",
    "    ma = [[mu0 * H[i], M[i]] for i in range(len(H))]\n",
    "    mp = rad.MatSatIsoTab(ma)\n",
    "    # -- permanent magnet material\n",
    "    mm = rad.MatStd(smat, rm)\n",
    "\n",
    "    return mp, mm\n",
    "\n",
    "\n",
    "def getMagnMaterCompMvsH(meshH, ind, cmpnH,  cmpnM):\n",
    "    \"\"\"\n",
    "    extract magnetization M_comp vs field strength H_comp from material data\n",
    "    arguments:\n",
    "      meshH = how to mesh the H axis: [hmin, hmax, nh]\n",
    "      ind   = magnet material index\n",
    "      cmpnH = component of H ('x', 'y', or 'z')\n",
    "      cmpnM = component of M ('x', 'y', or 'z')\n",
    "    return:\n",
    "      given component of magnetization along the given mesh in H\n",
    "    \"\"\"\n",
    "    hMin = meshH[0]; hMax = meshH[1]; nh = meshH[2]\n",
    "    hStep = (hMax - hMin) / (nh - 1)\n",
    "    m = [ 0 ] * nh\n",
    "    sCmpnM = 'm' + cmpnM\n",
    "    h = hMin\n",
    "    H = [0,0,0]\n",
    "    for i in range(nh):\n",
    "        if(  cmpnH == 'x'): H[0] = h\n",
    "        elif(cmpnH == 'y'): H[1] = h\n",
    "        elif(cmpnH == 'z'): H[2] = h\n",
    "        m[i] = rad.MatMvsH(ind, sCmpnM, H)\n",
    "        h += hStep\n",
    "    return m\n",
    "\n",
    "\n",
    "def calcField(grp, per, numPer):\n",
    "    \"\"\"\n",
    "    compute vertical magnetic field vs longitudinal position (Y)\n",
    "    arguments:\n",
    "      grp    = magnet (Radia group)\n",
    "      per    = period / mm\n",
    "      numPer = number of full magnetic periods\n",
    "    return:\n",
    "      Bz(y), meshY\n",
    "    \"\"\"\n",
    "    xc = 0 # transverse X coordinate\n",
    "    zc = 0 # transverse Z coordinate\n",
    "    yMax = per * (numPer + 1) / 2.\n",
    "    yMin = -yMax\n",
    "    ny = 501\n",
    "    yStep = (yMax - yMin) / ny\n",
    "    y = np.arange(yMin, yMin + ny * yStep, yStep)\n",
    "    \n",
    "    points = [[xc, yi, zc] for yi in y]\n",
    "    BzVsY = rad.Fld(grp, 'bz', points)\n",
    "\n",
    "    return BzVsY, [yMin, yMax, ny]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## _Create the undulator magnet_\n",
    "\n",
    "First enter and plot the data describing the magnetization curve for the iron pole tips."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ~iron type Va Permendur\n",
    "ironH = [  0.8,   1.5,     2.2,    3.6,    5.0,     6.8,     9.8,    18.0,\n",
    "          28.0,   37.5,   42.0,   55.0,   71.5,    80.0,    85.0,    88.0,\n",
    "          92.0,  100.0,  120.0,  150.0,  200.0,   300.0,   400.0,   600.0,\n",
    "         800.0, 1000.0, 2000.0, 4000.0, 6000.0, 10000.0, 25000.0, 40000.0]\n",
    "ironM = [0.000998995, 0.00199812, 0.00299724, 0.00499548, 0.00699372, 0.00999145, 0.0149877, 0.0299774,\n",
    "         0.0499648,   0.0799529,  0.0999472,  0.199931,   0.49991,    0.799899,   0.999893, 1.09989,\n",
    "         1.19988,     1.29987,    1.41985,    1.49981,    1.59975,    1.72962,    1.7995,   1.89925,\n",
    "         1.96899,     1.99874,    2.09749,    2.19497,    2.24246,    2.27743,    2.28958,  2.28973]\n",
    "\n",
    "# full M-H curve\n",
    "plt.figure(figsize=(7,4.3))\n",
    "plt.plot([mu0 * h for h in ironH], ironM)\n",
    "plt.title(r'Magneteization Curve $M(H)$')\n",
    "plt.xlabel(r'$\\mu_o H$ / T')\n",
    "plt.ylabel(r'$M$ / T')\n",
    "plt.grid()\n",
    "plt.show()\n",
    "\n",
    "# leading edge of M-H curve\n",
    "plt.figure(figsize=(7,4.3))\n",
    "plt.plot([mu0 * h for h in ironH], ironM)\n",
    "plt.title(r'Magneteization Curve $M(H)$')\n",
    "plt.xlabel(r'$\\mu_o H$ / T')\n",
    "plt.ylabel(r'$M$ / T')\n",
    "plt.xlim([0,0.001])\n",
    "plt.ylim([0,2.0])\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then set the magnet parameters and materials, and “build the magnet”."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set parameters for this undulator\n",
    "# -- general parameters\n",
    "numPer =  2    # number of full magnetic periods\n",
    "per = 46       # period / mm\n",
    "gap = 20       # gap height / mm\n",
    "gapOffset = 1  # \n",
    "\n",
    "# -- parameters for the iron poles\n",
    "lp =  [45,5,25] # dimensions / mm\n",
    "nsp = [2,2,5]   # pole-tip segmentation (try varying this setting)\n",
    "cp =  [1,0,1]   # assign color\n",
    "ll = per/2 - lp[1]\n",
    "\n",
    "# -- parameters for the magnet blocks\n",
    "lm  = [65,ll,45] # dimensions / mm\n",
    "nsm = [1,3,1]    # magnet-block segmentation (try varying this)\n",
    "cm  = [0,1,1]    # assign color\n",
    "\n",
    "# -- magnetic materials\n",
    "# pole tips: ~iron type Va Permendur\n",
    "# permanent magnets: NdFeB with 1.2 Tesla remanent magnetization\n",
    "mp, mm = materials(ironH, ironM, 'NdFeB', 1.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# then build the undulator\n",
    "und, pole, magnet = undulator(lp, mp, nsp, cp, lm, mm, nsm, cm, gap, gapOffset, numPer)\n",
    "\n",
    "# set up the radia viewer and display the magnet\n",
    "rv = radia_viewer.RadiaViewer()\n",
    "rv.add_geometry('Simple Dipole Steerer', und)\n",
    "rv.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# then build the undulator (work-around version)\n",
    "und, pole, magnet = undulator(lp, mp, nsp, cp, lm, mm, nsm, cm, gap, gapOffset, numPer)\n",
    "undX, pole, magnet = undulatorX(lp, mp, nsp, cp, lm, mm, nsm, cm, gap, gapOffset, numPer)\n",
    "\n",
    "# set up the radia viewer and display the magnet\n",
    "rvX = radia_viewer.RadiaViewer()\n",
    "rvX.add_geometry('Iron Poles',   pole)\n",
    "rvX.add_geometry('Permanent Magnets',   magnet)\n",
    "rvX.add_geometry('Simple Dipole Steerer',   und)\n",
    "rvX.add_geometry('Simple Dipole Steerer X', undX)\n",
    "rvX.display()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here plot the magnetization curve (extracted from Radia) for the iron pole tips. These should appear identical to the corresponding plots above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# full M-H curve\n",
    "meshH= [0.00, mu0 * ironH[-1], 201]\n",
    "Mpole = getMagnMaterCompMvsH(meshH, pole, 'y', 'y')\n",
    "\n",
    "plt.figure(figsize=(7,4.3))\n",
    "plt.plot(np.linspace(meshH[0],meshH[1],meshH[2]), Mpole)\n",
    "plt.title(r'Magneteization Curve $M(H)$')\n",
    "plt.xlabel(r'$\\mu_oH$ / T')\n",
    "plt.ylabel(r'$M$ / T')\n",
    "plt.grid()\n",
    "plt.show()\n",
    "\n",
    "# leading edge of the M-H curve\n",
    "meshH= [0.000, 0.001, 101]\n",
    "Mpole = getMagnMaterCompMvsH(meshH, pole, 'y', 'y')\n",
    "\n",
    "plt.figure(figsize=(7,4.3))\n",
    "plt.plot(np.linspace(meshH[0],meshH[1],meshH[2]), Mpole)\n",
    "plt.title(r'Magneteization Curve $M(H)$')\n",
    "plt.xlabel(r'$\\mu_oH$ / T')\n",
    "plt.ylabel(r'$M$ / T')\n",
    "plt.xlim([0,0.001])\n",
    "plt.ylim([0,2.0])\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## _Solve for the undulator magnetization_\n",
    "\n",
    "Here we ask Radia to solve for the undulator magnetization—_i.e._ for the magnetization within each sub-block that constitutes the undulator—and also report some statistics concerning that solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# solve for the magnetization\n",
    "prec = 0.0003  # precision for this computation\n",
    "maxIter = 1000 # maximum allowed iterations\n",
    "t0 = tm.time()\n",
    "res = rad.Solve(und, prec, maxIter)\n",
    "t1 = tm.time()\n",
    "\n",
    "print(\"Solved for magnetization in {0:6f} s\".format(t1 - t0))\n",
    "print(\"Relaxation results\")\n",
    "print(\"  number of iterations: {0:5d}\".format(int(res[3])))\n",
    "print(\"  average stability of magnetization at last iteration: {0:.4e} T\".format(res[0]))\n",
    "print(\"  maximum absolute magnetization at last iteration: {0:.5f} T\".format(res[1]))\n",
    "print(\"  maximum H vector at last iteration: {0:.5f} T\".format(res[2]))\n",
    "print(\"Peak magnetic field: {0:.8f} T\".format(rad.Fld(und, 'bz', [0,0,0])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now ask Radia to compute the magnetic field along the undulator’s longitudinal axis. And then plot the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute magnetic field values\n",
    "t2 = tm.time()\n",
    "BzVsY, MeshY = calcField(und, per, numPer)\n",
    "t3 = tm.time()\n",
    "print(\"Computed requested field vaues in {0:6f} s\\n\".format(t3 - t2))\n",
    "\n",
    "# and plot\n",
    "y = np.linspace(MeshY[0], MeshY[1], MeshY[2])\n",
    "plt.figure(figsize=(7,4.3))\n",
    "plt.plot(y, BzVsY, linewidth = 2.0)\n",
    "plt.title(r\"Vertical Magnetic Field, $B_z$, along undulator's $y$ axis\")\n",
    "plt.xlabel('y / mm')\n",
    "plt.ylabel('$B_z$ / T')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## _Check the Material Properties_\n",
    "\n",
    "Here we extract from Radia, and plot, the magnetization curve—magnetization _vs._ field strength—of the pole tip material."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract data for the pole tip material\n",
    "meshH_pole = [-0.0012, 0.0012, 201]\n",
    "M_pole = getMagnMaterCompMvsH(meshH_pole, pole, 'z', 'z')\n",
    "\n",
    "# and plot that data\n",
    "h = np.linspace(1e3*meshH_pole[0], 1e3*meshH_pole[1], meshH_pole[2])\n",
    "plt.figure(figsize=(7,4.3))\n",
    "plt.plot(h, M_pole, linewidth = 2.0)\n",
    "plt.title(r'Magnetization Curve $M(H)$ for pole material')\n",
    "plt.xlabel('magnetic field strength $\\mu_0 H$ / mT')\n",
    "plt.ylabel('magnetization $M$ / T')\n",
    "plt.grid()\n",
    "plt.show()\n",
    "\n",
    "# do this again over a wider range of H values\n",
    "meshH_pole = [-0.01, 0.01, 201]\n",
    "M_pole = getMagnMaterCompMvsH(meshH_pole, pole, 'x', 'x')\n",
    "\n",
    "h = np.linspace(1e3*meshH_pole[0], 1e3*meshH_pole[1], meshH_pole[2])\n",
    "plt.figure(figsize=(7,4.3))\n",
    "plt.plot(h, M_pole, linewidth = 2.0)\n",
    "plt.title(r'Magnetization Curve $M(H)$ for pole material')\n",
    "plt.xlabel('magnetic field strength $\\mu_0 H$ / mT')\n",
    "plt.ylabel('magnetization $M$ / T')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In a similar manner, we extract from Radia, and plot, magnetization curves for the permanent magnet material. In this case, the material has a significant anisotropy, as evidenced by the diffferent results parallel and perpendicular to the so-called “easy” axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract data for the permanent magnet material\n",
    "meshH_mag = [-1, 1, 201]\n",
    "Mpara_mag = getMagnMaterCompMvsH(meshH_mag, magnet, 'y', 'y') # parallel to easY axis\n",
    "Mperp_mag = getMagnMaterCompMvsH(meshH_mag, magnet, 'x', 'x') # perpendicular to easY axis\n",
    "\n",
    "# plot parallel result\n",
    "h = np.linspace(meshH_mag[0], meshH_mag[1], meshH_mag[2])\n",
    "plt.figure(figsize=(7,4.3))\n",
    "plt.plot(h, Mpara_mag, linewidth = 2.0)\n",
    "plt.title(r'Magnetization $M(H)$ parallel to \"easy\" axis')\n",
    "plt.xlabel('magnetic field strength $\\mu_0 H$ / T')\n",
    "plt.ylabel('magnetization $M$ / T')\n",
    "plt.grid()\n",
    "plt.show()\n",
    "\n",
    "# plot perpendicular result\n",
    "plt.figure(figsize=(7,4.3))\n",
    "plt.plot(h, Mperp_mag, linewidth = 2.0)\n",
    "plt.title(r'Magnetization $M(H)$ perpendicular to \"easy\" axis')\n",
    "plt.xlabel('magnetic field strength $\\mu_0 H$ / T')\n",
    "plt.ylabel('magnetization $M$ / T')\n",
    "plt.grid()\n",
    "plt.show()\n",
    "\n",
    "# plot both results\n",
    "plt.figure(figsize=(7,4.3))\n",
    "plt.plot(h, Mpara_mag, linewidth = 2.0)\n",
    "plt.plot(h, Mperp_mag, linewidth = 2.0)\n",
    "plt.title(r'$M(H)$ parallel (blue) and perpendicular (orange) to \"easy\" axis')\n",
    "plt.xlabel('magnetic field strength $\\mu_0 H$ / T')\n",
    "plt.ylabel('magnetization $M$ / T')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## _Optimizing Parameters_\n",
    "\n",
    "In this section, we examine the peak field of the undulator as a function\n",
    "of the horizontal width of the iron pole tips. This field value is needed\n",
    "in calculating the fundamental undulator parameter K, given by\n",
    "$$K = \\frac{e B_0 \\lambda_u}{2\\pi m_e c}.$$\n",
    "The function `maxField` defined here computes the maximum field for a given\n",
    "pole width. The other required arguments allow one to vary the segmentation\n",
    "of both pole-tips and magnet blocks. In the computations that follow, you\n",
    "will see that the choice of segmentation can have a significant impact on\n",
    "the accuracy of values computed for $B_0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def maxField(pw, np, nm, prec = 0.0005, maxIter = 1000):\n",
    "    \"\"\"\n",
    "    compute maximum vertical field for given undulator design\n",
    "    NB: Except for the argument values here, all magnet\n",
    "        and simulation parameters are as defined above\n",
    "    arguments:\n",
    "      pw      = pole width / mm\n",
    "      np      = segmentation of the iron poles, [nx, ny, nz]\n",
    "      nm      = segmentation of the magnet blocks\n",
    "      prec    = precision goal for this computation\n",
    "      maxIter = maximum allowed iterations\n",
    "    return: maximum vertical field / T\n",
    "    \"\"\"\n",
    "    rad.UtiDelAll()\n",
    "    lp0 = lp[0]\n",
    "    lp[0] = pw\n",
    "    mp, mm = materials(ironH, ironM, 'NdFeB', 1.2)\n",
    "    und, pole, magnet = undulator(lp, mp, np, cp, lm, mm, nm, cm, gap, gapOffset, numPer)\n",
    "    res = rad.Solve(und, prec, maxIter)\n",
    "    lp[0] = lp0\n",
    "    return rad.Fld(und, 'bz', [0,0,0])\n",
    "\n",
    "def undulatorK(b0, wl):\n",
    "    \"\"\"\n",
    "    compute undulator K value\n",
    "    arguments:\n",
    "      b0 = peak field / T\n",
    "      wl = undulator period / m\n",
    "    return:\n",
    "      K = (e B_0 \\lambda_u) / (2\\pi m_e c)\n",
    "    \"\"\"\n",
    "    return sc.e * b0 * wl / (2 * sc.pi * sc.m_e * sc.c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "undulatorK(0.491812, 46.e-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Setting both segmentation arrays to `[1,1,1]` results in a quick, but\n",
    "inaccurate, computation. Still, this choice gives us a general range to\n",
    "search for an optimum value of the magnet width."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pw = np.linspace(4, 60, (14 + 1))\n",
    "Bmax = [maxField(w, [1,1,1], [1,1,1], prec = 0.0005) for w in pw]\n",
    "\n",
    "# plot Bmax vs Pole width\n",
    "plt.figure(figsize=(7,4.3))\n",
    "plt.plot(pw, Bmax, linewidth = 2.0)\n",
    "plt.title(r'maximum $B$ field versus pole width')\n",
    "plt.xlabel('pole width $w$ / mm')\n",
    "plt.ylabel('$B_\\text{max}$ / T')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now increase both segmentation arrays, and then compute $B_0$ in a narrower\n",
    "range of pole-tip widths."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pw = np.linspace(24, 48, (16 + 1))\n",
    "Bmax = [maxField(w, [2,2,5], [1,3,1], prec = 0.0003) for w in pw]\n",
    "\n",
    "# plot Bmax vs Pole width\n",
    "plt.figure(figsize=(7,4.3))\n",
    "plt.plot(pw, Bmax, linewidth = 2.0)\n",
    "plt.title(r'maximum $B$ field versus pole width')\n",
    "plt.xlabel('pole width $w$ / mm')\n",
    "plt.ylabel('$B_\\text{max}$ / T')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can experiment with further increases to the segmentation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pw = np.linspace(30, 40, (20 + 1))\n",
    "Bmax = [maxField(w, [4,4,10], [2,6,2], prec = 0.0001) for w in pw]\n",
    "\n",
    "# plot Bmax vs Pole width\n",
    "plt.figure(figsize=(7,4.3))\n",
    "plt.plot(pw, Bmax, linewidth = 2.0)\n",
    "plt.title(r'maximum $B$ field versus pole width')\n",
    "plt.xlabel('pole width $w$ / mm')\n",
    "plt.ylabel('$B_\\text{max}$ / T')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pw = np.linspace(36, 40, (16 + 1))\n",
    "Bmax = [maxField(w, [4,4,10], [2,6,2], prec = 0.0001) for w in pw]\n",
    "\n",
    "# plot Bmax vs Pole width\n",
    "plt.figure(figsize=(7,4.3))\n",
    "plt.plot(pw, Bmax, linewidth = 2.0)\n",
    "plt.title(r'maximum $B$ field versus pole width')\n",
    "plt.xlabel('pole width $w$ / mm')\n",
    "plt.ylabel('$B_\\text{max}$ / T')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pw = np.linspace(36, 40, (16 + 1))\n",
    "Bmax = [maxField(w, [4,4,11], [2,6,2], prec = 0.0001) for w in pw]\n",
    "\n",
    "# plot Bmax vs Pole width\n",
    "plt.figure(figsize=(7,4.3))\n",
    "plt.plot(pw, Bmax, linewidth = 2.0)\n",
    "plt.title(r'maximum $B$ field versus pole width')\n",
    "plt.xlabel('pole width $w$ / mm')\n",
    "plt.ylabel('$B_\\text{max}$ / T')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pw = np.linspace(36, 40, (16 + 1))\n",
    "Bmax = [maxField(w, [4,4,11], [2,6,4], prec = 0.0001) for w in pw]\n",
    "\n",
    "# plot Bmax vs Pole width\n",
    "plt.figure(figsize=(7,4.3))\n",
    "plt.plot(pw, Bmax, linewidth = 2.0)\n",
    "plt.title(r'maximum $B$ field versus pole width')\n",
    "plt.xlabel('pole width $w$ / mm')\n",
    "plt.ylabel('$B_\\text{max}$ / T')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pw = np.linspace(36, 39, (12 + 1))\n",
    "Bmax = [maxField(w, [6,6,11], [4,4,4], prec = 0.0001) for w in pw]\n",
    "\n",
    "# plot Bmax vs Pole width\n",
    "plt.figure(figsize=(7,4.3))\n",
    "plt.plot(pw, Bmax, linewidth = 2.0)\n",
    "plt.title(r'maximum $B$ field versus pole width')\n",
    "plt.xlabel('pole width $w$ / mm')\n",
    "plt.ylabel('$B_\\text{max}$ / T')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "py3"
  },
  "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.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}