Dispersion-Chromaticity.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![alt text](images/uspas.png)\n",
    "# Fundamentals of Accelerator Physics and Technology \n",
    "### (with Simulations and Measurements Lab)\n",
    "# Computer Lab: Dispersion and Chromaticity in Simulated Beam Transport\n",
    "##### Authors: K. Ruisard, N. Evans, N. Neveu and L. Dovlatyan\n",
    "\n",
    "## We will be simulating beam transport in simple beamlines. To do this, we will be running the Elegant code in the [Sirepo](https://beta.sirepo.com) environment. Questions to be turned in for credit are in **bold** and numbered.\n",
    "### Python Notes: \n",
    "- You wont need much python specific code for this lab, but here are some tips:\n",
    "- Press shift+enter to execute a cell, or use the play button at the top of the window\n",
    "- Make sure you execute cells in order, or re-execute cells if you change something at the top of the notebook.\n",
    "- You can also execute the whole notebook by using 'Run all cells' under the 'Run' tab.\n",
    "\n",
    "</br>\n",
    "Also helpful: Shift+right click brings up OS/browser right-click menu, can copy image or save.\n",
    "\n",
    "----------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 1. Dispersion and Chromaticity\n",
    "\n",
    "### A) Dispersion in FODO lattice\n",
    "\n",
    "In other labs, we simulate beamlines without dipoles. Dipoles create dispersion ($\\eta$ or D, depending who you talk to), which describes the evolution of transverse size for particles/beams that are off-momentum. Off-momentum particles experience a different force in the dipole than the design orbit, as they have different magnetic rigidity $B\\rho=\\frac{p}{q}$. Dispersion is a function of s and can have a periodic (matched) solution, similar to the $\\beta$ and $\\alpha$ functions.\n",
    "\n",
    "The initial conditions for this simulation:\n",
    "\n",
    "| Species  | Election  |   \n",
    "|---|---| \n",
    "| Energy  | 1 GeV  |   \n",
    "|  X emittance | $\\epsilon_x = 6$ mm-mrad  |  \n",
    "|  Y emittance |  $\\epsilon_y = 6$ mm-mrad |  \n",
    "|  Quadrupole geometric strength | $K = 0.6$ m$^{-2}$|\n",
    "|  FODO cell length | $L = 5$ m|\n",
    "\n",
    "Open simulation “FODObeamline-with-dipole” \n",
    "> https://www.sirepo.com/elegant#/source/SMHkFUFD\n",
    "\n",
    "Initially, this simulation lattice is a FODO cell that is identical to the cell used in the FODO-transport lab. \n",
    "\n",
    "- Go to Visualization tab at top of the page.\n",
    "- Run simulation to run over 1 cell (make sure beamline \"FODOcell\" is selected). \n",
    "- verify that calculated solution is matched. If not, under Tab: Control -> twiss_output, Change Matched to Yes and re-run simulation. \n",
    "- Note down matched beam size (Sx and Sy in the run_setup.sigma plot) at the location of the quads QF and QD for future reference. \n",
    "\n",
    "\n",
    "Next, add a dipole element to your FODO cell: \n",
    "- Go to the 'Lattice Tab' at the top of the screen\n",
    "- Under 'Beamlines', hover over FODOcell, and click 'Edit'\n",
    "- Delete the 0.5 meter drift \"D3\" and, in its place, add the RBEN element \"DIPO\" \n",
    "- \"DIPO\" is a 20-degree dipole with rectangular geometry (as opposed to sector geometry, which is also an option).\n",
    "- Run simulation over modified FODO cell\n",
    "\n",
    "\n",
    "#### A comment on dipole edge focusing\n",
    "First, scroll to the sigma_output plot to observe the effect of the dipole on the beam moments.\n",
    "\n",
    "Note that the sigma_output plot is generated by tracking a particle distribution through the cell. This distribution does not include a momentum spread, so dispersion should not affect the results of this plot. \n",
    "\n",
    "However, you will observe that Sx looks slightly different from the case without a dipole. There is a small de-focusing kick at the entrance of the dipole, and a small focusing kick at the exit of the dipole. These two kicks cancel, and the periodic solution in the horizontal plane is almost identical to the FODO cell without a dipole. \n",
    "\n",
    "In contrast, the average value of Sy has decreased significantly, and there is a focusing kick at the entrance and exit. The dipole has a net focusing effect in the vertical plane. \n",
    "\n",
    "The reason for these differences is an effect called *edge focusing*, which is caused by fringe fields that extend past the edge of the dipole. An exploration of edge focusing is beyond the scope of this lab, but it is good to be aware that this is another way in which dipoles affect the lattice and beam parameters. Edge focusing is an effect which does not depend on beam energy spread. In contrast, the strength of the effects caused by dispersion and chromaticity do depend on energy spread in the particle distribution - which is mainly determined by the RF system used for acceleration.\n",
    "\n",
    "#### Effect of dispersion on beam size\n",
    "\n",
    "Scroll down and to the twiss_output plot.\n",
    "Observe dispersion $\\eta_x$ (etax). \n",
    "\n",
    "hint: to see dispersion function at an appropriate scale, turn off $\\beta_x$ and $\\beta_y$ curves by un-checking the boxes.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q1) What is minimum dispersion in this lattice? At what location does this occur (focusing quad, defocusing quad or drift)?\\\n",
    "If we compare this lattice to an identical FODO cell without the dipole, would the horizontal beam size at this location be larger or smaller?  \n",
    "What about the vertical beam size?**\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q2) Assuming a 0.1% momentum spread ($\\frac{\\delta p}{p}=0.001$) in the beam, what is the horizontal beam size we expect in the focusing quadrupole QF? How does this compare to our beam size without energy spread? \\\n",
    "What is the vertical beam size in QF with and without energy spread?**\n",
    "\n",
    "Use max eta and the $\\epsilon_x$ from above.\n",
    "\n",
    "you may find this formula helpful: $\\sigma_x^2 = \\epsilon_x \\beta_x+ \\eta_{max}^2 \\left( \\frac{\\Delta p}{p_0} \\right)^2$, where $\\Delta p/p $ is the fractional momentum spread.\n",
    "\n",
    "You can check your calculation against simulation by introducing a momentum spread. Under the \"Control\" tab, edit \"bunched_beam\", Page 3, \"Sigma DP,\" where Sigma DP $ = \\Delta p / p_0$. The result will not be identical but should be within ~0.2 mm. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "|           | $\\delta = 0$ | $\\delta = 0.001$ |  \n",
    "|---        |---           |---               | \n",
    "|$\\sigma_x$ |              |                  |\n",
    "|$\\sigma_y$ |              |                  |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### B) Designing a zero-dispersion insert\n",
    "\n",
    "Open simulation DispersionFree:\n",
    "> https://www.sirepo.com/elegant#/source/dW4biLsM\n",
    "\n",
    "The lattice cell in this simulation has two 18-degree bends and five quadrupoles.  \n",
    "The quadrupole fields are initially set to zero and matching is disabled.  \n",
    "- Run the simulation and observe the evolution of the dispersion $\\eta_x$ through the double bend.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q3) What is $\\eta_x$ and $\\eta_x'$ at the end of the cell?**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Go to 'Lattice' at the top of the page.\n",
    "- “Turn on” the middle quadrupole (Q1) at $k_1=1\\ m^{-2}$. Note the effect this has on the dispersion function.\n",
    "- Adjust the strength of Q1 to find the solution for which dispersion is zero after the two bends.\n",
    "- Hint: you will need to restart the simulation after you change K. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q4) What is the strength you found?** Give your answer with at least two decimal places.\n",
    "\n",
    "\n",
    "$k_1$ ="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Enable matching (Under the Control tab -> twiss_output command -> set “Matched” to “yes”).  \n",
    "- Run simulation; you should see that beamline is unstable, i.e. there are no periodic solutions. \n",
    "Let’s add some additional focusing to correct that.\n",
    "- “Turn on” the flanking quadrupoles by setting Q2: $k_1=1.33\\ m^{-2}$ and Q3: $k_1=-1.59\\ m^{-2}$. \n",
    "- Re-run to find matched solution.\n",
    "\n",
    "Congrats! You’ve made a double-bend achromat. This is useful for light sources, where sections with zero dispersion are used for placing radiation-generating inserts like undulators and wigglers. Note that dispersion is cancelled for only part of the cell - it is not zero everywhere. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q5) What is the maximum dispersion in this cell?**\n",
    "\n",
    "$\\eta_x = $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Momentum acceptance limited by dispersion\n",
    "\n",
    "Because the effect of dispersion depends on the beam momentum spread, it presents a limit to the amount of spread that can be transported without losing particles.\n",
    "\n",
    "**Q6) Assuming the vacuum chamber is a pipe with 2.5 cm radius, what is the largest momentum spread that can be tolerated before particles start to hit the chamber wall ?** \n",
    "\n",
    "You will need to use this equation again: $\\sigma_x^2 = \\epsilon_x \\beta_x+ \\eta^2 \\left( \\frac{\\Delta p}{p_0} \\right)^2$ \n",
    "\n",
    "**Q7) If momentum spread exceeds this limit, where in the cell will beam loss occur?**\n",
    "\n",
    "You can check your answers to Q6) and Q7) against simulation by introducing a momentum spread into the beam. To do this, follow the same instructions as under Q2). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### C) Chromaticity in a ring\n",
    "\n",
    "Repeat your cell 10 times to create a ring. \n",
    "- Under the Lattice tab, make a new beamline. \n",
    "- Fill it with DBA elements.\n",
    "- Rerun the simulation under the Visualization tab.\n",
    "- Make sure you select your new beamline from the drop down menu.\n",
    "\n",
    "\n",
    "Although sections of this lattice are dispersion-free, this ring still has chromaticity. In other words, although off-momentum particles do not increase the beam size in the dispersion-free sections, they still feel a different focusing strength than the *synchronous* or on-energy particle and, as a results, have a slightly different *tune*. Recall that tune is the frequency of particle oscillation around the design particle.  \n",
    "\n",
    "First we will calculate the tune of the on-energy particles, then calculate the spread of tune that results from momentum spread. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "**Q8) Record x and y tunes (to 3 significant figures):**\n",
    "\n",
    "The tunes are \"nux\" and \"nuy\" in the \"Output Parameters\" list\n",
    "\n",
    "$\\nu_x = $   \n",
    "$\\nu_y = $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q9) In our FODO cell, we apply equal focusing strengths in both planes and as a result, $\\nu_x=\\nu_y$. Looking at the beta functions for one DBA cell, explain why you would not necessarily expect the tunes to be equal.**\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q10) For a 0.1% momentum spread in the beam, what is the spread of tunes due to chromaticity?**\n",
    "\n",
    "\n",
    "In the Output Parameters list, look for dnux/dp and dnuy/dp which are the chromaticity values: \n",
    "\n",
    "$C_x = $\n",
    "\n",
    "$C_y = $\n",
    "\n",
    "Use $\\Delta \\nu = C \\frac{\\Delta p}{p_0}$  to calculate tune spread:\n",
    "\n",
    "$\\Delta \\nu_x = $\n",
    "\n",
    "$\\Delta \\nu_y = $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### Momentum acceptance limit due to chromaticity\n",
    "\n",
    "When operating a ring, we have to carefully choose the tune to avoid resonant conditions, which may cause instability in the particle trajectories. Small errors in our lattice can resonate with the beam motion when there is a rational relationship between the tunes. \n",
    "\n",
    "A rational relationship is one that satisfies $m \\nu_x + n \\nu_y = p$ for any integers $m$,$n$,$p$. The order of the resonance is $|m|+|n|$. Generally speaking, driving terms for higher-order resonances are weaker. \n",
    "\n",
    "Which resonance lines and orders of resonance need to be avoided are unique to each accelerator facility. Here we will assume that the ring we just designed needs to be operated away from all resonances at 3rd order and below. We will use the plot below to visualize how the chromatic tune spread might affect our operation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "## This cell includes the function that makes the plot below. You don't need to look at it, just execute it once.\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.lines import Line2D\n",
    "import numpy as np\n",
    "\n",
    "def tunediagram(order=range(1,4),integer=[0,0],lines=[1,1,1,1],colors='ordered',linestyle='-',fig=plt.gcf(),legend=True):\n",
    "    '''\n",
    "    plot resonance diagram up to specified order\n",
    "    mx + ny = p\n",
    "    x = (p-ny)/m\n",
    "    x = 1 where y = (p-m)/n\n",
    "    EXAMPLE:\n",
    "       tunediagram(order=[1,2,3])\n",
    "       tunediagram([1,2,3],integer=[6,8],lines=[0,0,1,1])\n",
    "       tunediagram([1,2,3],integer=[6,8],lines=[0,0,1,1], colors='black3', linestyle='--')\n",
    "    INPUT:\n",
    "      1. order - array of tune orders to plot. e.g. up to 3rd order, [1,2,3]\n",
    "      2. integers - integer part of the tune to plot in x,y, default is [0,0]. e.g. plot from 6-7 and 9-10, integer=[6,9]\n",
    "      2. lines - a boolean of which resonance lines to plot. e.g. [vertical,horizontal,sum,diff]. e.g. plot only vert/horz lines, lines = [1,1,0,0]\n",
    "      4. colors - option to plot lines in different colors. default is ordered. color options only go up to 10th order\n",
    "        ordered - each order resonance is a different color\n",
    "        black - all lines are black\n",
    "        blackX - X here is a number. all resonances X+ will be in black. e.g. black3 means plot resonances 1-2 in color and 3,4,5,... in black\n",
    "      6. linestyle - linestyle option from matplotlib\n",
    "      7. fig - pass in a handle to a figure, otherwise things will be plotted in the current figure.\n",
    "    Written by Levon Dovlatyan\n",
    "    University of Maryland, Department of Physics\n",
    "    Oct 2018\n",
    "    '''\n",
    "    # define some variables\n",
    "    pval = 40 # increase for more/higher order lines\n",
    "    p = np.linspace(0,pval,pval+1)\n",
    "    \n",
    "    qxmin,qymin = integer[0],integer[1]\n",
    "    qxmax,qymax = qxmin+1,qymin+1\n",
    "    \n",
    "    # define different colors, up to 10th order\n",
    "    color = ['C0','C1','C2','C3','C4','C5','C6','C7','C8','C9']\n",
    "    if colors == 'black':\n",
    "        color = ['k']*10\n",
    "    elif colors[0:-1] == 'black':\n",
    "        idx = int(colors[-1])\n",
    "        color = color[0:idx-1] + (['k']*10)[idx-1:]\n",
    "    \n",
    "    # adjust plot limits\n",
    "    plt.xlim((qxmin-0.01, qxmax+0.01))\n",
    "    plt.ylim((qymin-0.01, qymax+0.01))\n",
    "    \n",
    "    # hold legend objects\n",
    "    legend_elements = []\n",
    "    \n",
    "    # Plotting formula\n",
    "    # we plot resonances in reverse order\n",
    "    for i in order[::-1]:\n",
    "        # -- legend\n",
    "        legend_elements.append(Line2D([0], [0], color=color[i-1], linestyle=linestyle, label='order %i'%i))\n",
    "        \n",
    "        # -- resonance line\n",
    "        m = np.linspace(-i,i,2*i+1)\n",
    "        n1 = (i-np.abs(m))\n",
    "        n2 = -1*n1\n",
    "        for j in range(0,m.size,1):\n",
    "            # check to see equation is divided by 0 or not\n",
    "            # ver & hor res lines\n",
    "            if ((n1[j] == 0 and lines[1]) or (m[j] == 0 and lines[0])):\n",
    "                # vertical lines\n",
    "                if n1[j] == 0 and lines[1]:\n",
    "                    plt.vlines(p/m[j],qymin,qymax,color=color[i-1],linestyle=linestyle);\n",
    "                # horizontal lines\n",
    "                if m[j] == 0 and lines[0]:\n",
    "                    plt.hlines(p/n1[j],qxmin,qxmax,color=color[i-1],linestyle=linestyle);\n",
    "                    plt.hlines(p/n2[j],qxmin,qxmax,color=color[i-1],linestyle=linestyle);\n",
    "            # sum and dif res lines\n",
    "            elif not(n1[j] == 0) and not(m[j] == 0):\n",
    "                # resonance sum lines\n",
    "                if lines[2]:\n",
    "                    if np.sign(m[j]) > 0:\n",
    "                        plt.plot([[qxmin]*p.size,[qxmax]*p.size],[p/n2[j] - np.array(m[j]*qxmin/n2[j]), p/n2[j] - np.array(m[j]*qxmax/n2[j])],color=color[i-1],linestyle=linestyle);\n",
    "                    else:\n",
    "                        plt.plot([[qxmin]*p.size,[qxmax]*p.size],[p/n1[j] - np.array(m[j]*qxmin/n1[j]), p/n1[j] - np.array(m[j]*qxmax/n1[j])],color=color[i-1],linestyle=linestyle);\n",
    "                # resonance dif lines\n",
    "                if lines[3]:\n",
    "                    if np.sign(m[j]) > 0:\n",
    "                        plt.plot([[qxmin]*p.size,[qxmax]*p.size],[p/n1[j] - np.array(m[j]*qxmin/n1[j]), p/n1[j] - np.array(m[j]*qxmax/n1[j])],color=color[i-1],linestyle=linestyle);\n",
    "                    else:\n",
    "                        plt.plot([[qxmin]*p.size,[qxmax]*p.size],[p/n2[j] - np.array(m[j]*qxmin/n2[j]), p/n2[j] - np.array(m[j]*qxmax/n2[j])],color=color[i-1],linestyle=linestyle);\n",
    "\n",
    "    if legend:\n",
    "        legend_elements.reverse()\n",
    "        plt.legend(handles=legend_elements, loc='upper left',bbox_to_anchor=(1,1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# This cell creates a resonance diagram with the tune footprint for our ring indicated.\n",
    "# for this cell to work, make sure you execute the previous cell.\n",
    "\n",
    "##########################\n",
    "## Change these parameters\n",
    "##########################\n",
    "nux = 0\n",
    "nuy = 0\n",
    "sigma_dp = 0\n",
    "Cx = 0\n",
    "Cy = 0\n",
    "resonance_order = 2\n",
    "\n",
    "##########################\n",
    "## don't change below\n",
    "##########################\n",
    "delta_nux = Cx*sigma_dp\n",
    "delta_nuy = Cy*sigma_dp\n",
    "\n",
    "\n",
    "#from tunediagram import tunediagram\n",
    "import matplotlib.pyplot as plt\n",
    "plt.figure(figsize=[6,6])\n",
    "\n",
    "# -- plot resonance lines\n",
    "tunediagram(order=range(1,resonance_order+1),integer=[int(nux),int(nuy)])\n",
    "plt.title('Resonance diagram')\n",
    "plt.ylabel('vertical tune'); plt.xlabel('horizontal tune');\n",
    "\n",
    "# -- plot tune footprint over resonance lines\n",
    "if (nux!=0) & (nuy!=0):\n",
    "    plt.plot([nux],[nuy],'ko',label='ring tune')\n",
    "    plt.plot([nux+delta_nux,nux-delta_nux],[nuy+delta_nuy,nuy-delta_nuy],'k.')\n",
    "    plt.plot([nux-delta_nux,nux+delta_nux],[nuy-delta_nuy,nuy+delta_nuy],'k-');\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q11) Assuming we can tolerate resonances of order 4 and above, but cannot operate on resonances of order 3 and below, how big of a momentum spread can we tolerate in this ring?** Using the plot above, set `resonance_order=3` and plug in the values for tunes and chromaticity that you recorded above. Adjust the momentum spread `sigma_dp` until your tune spread (the black line) crosses a resonance line. Give your answer to the nearest 0.1%\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q12) Comparing your answers to Q6 and Q11, is the momentum acceptance of this ring limited by chromaticity or by dispersion?**"
   ]
  }
 ],
 "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
}