From a6e3ae3d332926a0f97bebc3447e84ae1c0d4fff Mon Sep 17 00:00:00 2001 From: Cristian Date: Tue, 1 Jun 2021 21:41:50 +0100 Subject: [PATCH] Draft version --- .../IBS_mod/ELENA_PyHEAD_IBS_Example.ipynb | 671 ++++++++++++++++++ .../IBS_mod/ELENA_elements.def | 182 +++++ .../IBS_mod/ELENA_ring_mod.seq | 197 +++++ simulation_notebooks/IBS_mod/PyHEAD_IBS.py | 627 ++++++++++++++++ .../IBS_mod/Quadrupole_settings | 373 ++++++++++ .../IBS_mod/elena_simplified_origin.madx | 122 ++++ 6 files changed, 2172 insertions(+) create mode 100644 simulation_notebooks/IBS_mod/ELENA_PyHEAD_IBS_Example.ipynb create mode 100644 simulation_notebooks/IBS_mod/ELENA_elements.def create mode 100644 simulation_notebooks/IBS_mod/ELENA_ring_mod.seq create mode 100755 simulation_notebooks/IBS_mod/PyHEAD_IBS.py create mode 100644 simulation_notebooks/IBS_mod/Quadrupole_settings create mode 100644 simulation_notebooks/IBS_mod/elena_simplified_origin.madx diff --git a/simulation_notebooks/IBS_mod/ELENA_PyHEAD_IBS_Example.ipynb b/simulation_notebooks/IBS_mod/ELENA_PyHEAD_IBS_Example.ipynb new file mode 100644 index 0000000..43058aa --- /dev/null +++ b/simulation_notebooks/IBS_mod/ELENA_PyHEAD_IBS_Example.ipynb @@ -0,0 +1,671 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# general imports" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "from __future__ import division, print_function\n", + "\n", + "import numpy as np\n", + "np.random.seed(42)\n", + "\n", + "from scipy.constants import m_p, c, e\n", + "\n", + "import matplotlib.pyplot as plt\n", + "%matplotlib inline\n", + "\n", + "from IPython import display" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "# sets the PyHEADTAIL directory etc.\n", + "try:\n", + " from settings import *\n", + "except:\n", + " pass" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# PyHEADTAIL imports" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "PyHEADTAIL v1.14.1\n", + "\n", + "\n" + ] + } + ], + "source": [ + "from PyHEADTAIL.particles.generators import ParticleGenerator, RF_bucket_distribution, gaussian2D,gaussian2D_asymmetrical\n", + "from PyHEADTAIL.trackers.longitudinal_tracking import RFSystems\n", + "from PyHEADTAIL.trackers.rf_bucket import RFBucket\n", + "from PyHEADTAIL.trackers.transverse_tracking import TransverseMap\n", + "from PyHEADTAIL.monitors.monitors import BunchMonitor, SliceMonitor, ParticleMonitor\n", + "from PyHEADTAIL.trackers.detuners import Chromaticity, AmplitudeDetuning\n", + "import sys\n", + "from PyHEADTAIL.particles import generators,particles\n", + "from PyHEADTAIL.particles.generators import RFBucketMatcher\n", + "from PyHEADTAIL.cobra_functions.stats import emittance" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Setting up the machine and beam parameters" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "# General\n", + "\n", + "nturns = 100 #200 #20000 #2**16\n", + "npart =5000\n", + "A=1\n", + "Q=-1\n", + "\n", + "Ekin_per_nucleon = 100000 #0.2e9 # in eV\n", + "intensity =4.5e6 #0.625e11\n", + "nmass=0.9383\n", + "mass_eV = A * nmass * 1e9\n", + "mass_MeV = A * nmass * 1e3\n", + "mass = mass_eV * e / c**2 # in kg\n", + "charge = Q * e # in Coul\n", + "\n", + "Ekin = Ekin_per_nucleon * A\n", + "p0c = np.sqrt(Ekin**2 + 2*Ekin*mass/e * c**2) # in eV\n", + "\n", + "Etot = np.sqrt(p0c**2 + (mass/e)**2 * c**4) * 1e-9 # in GeV\n", + "p0 = p0c / c * e # in SI units\n", + "gamma = np.sqrt(1 + (p0 / (mass * c))**2)\n", + "beta = np.sqrt(1 - gamma**-2)\n", + "\n", + "# Machine\n", + "C = 30.4056\n", + "R = C/(2*np.pi)\n", + "\n", + "alpha = 0.25835\n", + "eta = alpha - 1/gamma**2\n", + "V = [100]\n", + "h = [1]\n", + "phi = [np.pi]\n", + "frequency1=0.143938 #MhZ\n", + "\n", + "# Acceleration\n", + "T0 = C/(beta*c)\n", + "normalisation = 1/C * e/p0 * T0\n", + "\n", + "epsx_rms_fin = 2.5e-6 #35e-6 / 4 # geometrical emittances\n", + "epsy_rms_fin = 2.5e-6 #15e-6 / 4\n", + "\n", + "limit_n_rms_x = 2 #3.5#2\n", + "limit_n_rms_y = 2#3.5#2\n", + "limit_n_rms_z = 2 #3.5#3.4\n", + "\n", + "sig_z = 0.3284 #58 / 4. # in m\n", + "sig_dp =1e-4 #0.5e-3e\n", + "epsx_gauss = epsx_rms_fin #* 1.778 RMS adjustment not needed if limit_n_rms are >>2!\n", + "epsy_gauss = epsy_rms_fin #* 1.82 RMS adjustment not needed if limit_n_rms are >>2!\n", + "\n", + "epsn_x = epsx_gauss * beta * gamma\n", + "epsn_y = epsy_gauss * beta * gamma\n", + "\n", + "beta_z = sig_z / sig_dp\n", + "\n", + "epsn_z = sig_z * sig_dp * 4 * np.pi * p0/e\n", + "\n", + "p_increment = sig_dp * p0\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Dictionary for IBS calculations" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "beam_dict={'Z':Q,\n", + " 'mass':mass_MeV,\n", + " 'KE':Ekin/1e+6,\n", + " 'emit_x':epsx_gauss,\n", + " 'emit_y':epsx_gauss,\n", + " 'dp_p':sig_dp,\n", + " 'N_ptcl':intensity,\n", + " 'sigma_s0':sig_z\n", + " }" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Preparing Twiss files" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + " ++++++++++++++++++++++++++++++++++++++++++++\n", + " + MAD-X 5.06.01 (64 bit, Linux) +\n", + " + Support: mad@cern.ch, http://cern.ch/mad +\n", + " + Release date: 2020.09.01 +\n", + " + Execution date: 2021.06.01 20:19:22 +\n", + " ++++++++++++++++++++++++++++++++++++++++++++\n", + "MAD-X summary: \n", + " Qx= 2.454 \n", + " Qy= 1.416 \n", + " dQx=-3.1574408714289386\n" + ] + } + ], + "source": [ + "from cpymad.madx import Madx\n", + "qx = 2.454 #18.95#84 #tune_range_qx[0]\n", + "qy = 1.416 #18.94#73 #tune_range_qy[0]\n", + "install_apertures = False\n", + "apply_errors = False\n", + "Thin=False\n", + "pystl_device = \"opencl:0.0\"\n", + "e_seed = 24\n", + "\n", + "qqx, qqy = int(np.round((qx%1) * 100)), int(np.round((qy%1) * 100))\n", + "filename_errors = \"errors_{qqx}_{qqy}_{eseed:d}\".format(\n", + " qqx=qqx, qqy=qqy, eseed=e_seed)\n", + "\n", + "madx=Madx()\n", + "madx.options.echo = False\n", + "madx.options.warn = False\n", + "madx = Madx(stdout=False)\n", + "madx.call('elena_simplified_origin.madx')\n", + "\n", + "\n", + "madx.command.beam(particle='ion', mass=A*nmass, charge=Q, energy=Etot) # /Q for RF to have proton beam!\n", + "\n", + "title=\"Thick elements\"\n", + "\n", + "madx.use(sequence='ELENA') \n", + " \n", + "if Thin==False:\n", + " \n", + " madx.input('''match, sequence=ELENA; \n", + " global, sequence=ELENA, q1={qx}, q2={qy};\n", + " vary, name=kq1, step=0.00001;\n", + " vary, name=kq2, step=0.00001;\n", + " lmdif, calls=500, tolerance=1.0e-10;\n", + " endmatch;\n", + " '''.format(qx=qx, qy=qy)\n", + " )\n", + " \n", + " twiss = madx.twiss();\n", + "\n", + "madx_seq = madx.sequence.ELENA\n", + "print(f\"MAD-X summary: \\n Qx= {twiss.summary.q1} \\n Qy= {twiss.summary.q2} \\n dQx={twiss.summary.dq1*beta}\")\n", + "qx_act=twiss.summary.q1\n", + "qy_act=twiss.summary.q2\n", + "dq1=twiss.summary.dq1*beta\n", + "dq2=twiss.summary.dq2*beta" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# IBS functions" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "from PyHEAD_IBS import IBS_Martini_rates,IBS_BM_rates,IBS_kick" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Rates X 0.22132545210160895 # Rates Y -0.31346560046982186 # Rates Z 51.2198147405692 #\n" + ] + } + ], + "source": [ + "IBS_Martini=IBS_Martini_rates(twiss,beam_dict,True,60,60,100,12.5,1)\n", + "print(f'Rates X {IBS_Martini.rsx} # Rates Y {IBS_Martini.rsy} # Rates Z {IBS_Martini.rss} #')" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Rates X 0.22069657967379055 # Rates Y -0.3062356075066439 # Rates Z 51.69187245553312 #\n" + ] + } + ], + "source": [ + "IBS_BM=IBS_BM_rates(twiss,beam_dict,True,60,60,100,12.5,1)\n", + "print(f'Rates X {IBS_BM.rsx} # Rates Y {IBS_BM.rsy} # Rates Z {IBS_BM.rss} #')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "chroma = Chromaticity(Qp_x=[dq1], Qp_y=[dq2]) # note the lists!" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Mean optics from MAD-X IBS calculations" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "betx = 4.20436E+00 bety = 3.48406E+00 Dx = 1.21510E+00 Dy = 0.00000E+00\n", + "alfx = 4.01651E-16 alfy = 4.81982E-16 Dpx = -4.38165E-17 Dpy = 0.00000E+00\n", + "1/betx = 5.05807E-01 1/bety = 2.93247E-01\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "n_segments = 1\n", + "C = twiss.summary.length\n", + "R = C / (2.*np.pi)\n", + "s = np.arange(0, n_segments + 1) * C / n_segments\n", + "\n", + "alpha_x_inj = 0.\n", + "alpha_y_inj = 0.\n", + "\n", + "beta_x_inj =4.20436 #R/qx_act\n", + "beta_y_inj =3.48406 #R/qy_act\n", + "\n", + "D_x_inj=1.2151\n", + "D_y_inj=0\n", + "\n", + "alpha_x = alpha_x_inj * np.ones(n_segments)\n", + "beta_x = beta_x_inj * np.ones(n_segments)\n", + "alpha_y = alpha_y_inj * np.ones(n_segments)\n", + "beta_y = beta_y_inj * np.ones(n_segments)\n", + "\n", + "D_x = np.ones(n_segments)*D_x_inj\n", + "D_y = np.ones(n_segments)*D_y_inj\n", + "\n", + "dq1=twiss.summary.dq1*beta\n", + "dq2=twiss.summary.dq2*beta\n", + "\n", + "trans_map_smooth = TransverseMap(\n", + "s, alpha_x, beta_x, D_x, alpha_y, beta_y, D_y, qx_act, qx_act, [chroma])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "longitudinal_map = RFSystems(\n", + " C, h, V, phi, [alpha], gamma, \n", + " p_increment=0, charge=-e, mass=m_p,D_x=D_x_inj, D_y=0)\n", + "\n", + "rfbucket = longitudinal_map.get_bucket(gamma=gamma)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "np.random.seed(3)\n", + "\n", + "pyht_beam = generators.generate_Gaussian6DTwiss(\n", + " npart, 1, charge, mass, C, gamma,\n", + " alpha_x_inj, alpha_y_inj,beta_x_inj, beta_y_inj,\n", + " beta_z, epsn_x, epsn_y, epsn_z,\n", + " dispersion_x=D_x_inj, #D_x_0 if D_x_0 else None,\n", + " dispersion_y=0, #D_y_0 if D_y_0 else None,\n", + " limit_n_rms_x=limit_n_rms_x**2, limit_n_rms_y=limit_n_rms_y**2, \n", + " limit_n_rms_z=limit_n_rms_z**2,\n", + ")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "epsn_x_ratio_sq = np.sqrt(epsn_x / pyht_beam.epsn_x())\n", + "epsn_y_ratio_sq = np.sqrt(epsn_y / pyht_beam.epsn_y())\n", + "dp_ratio=sig_dp / np.std( pyht_beam.dp)\n", + "sigz_ratio=sig_z/np.std( pyht_beam.z)\n", + "\n", + "pyht_beam.x *= epsn_x_ratio_sq\n", + "pyht_beam.xp *= epsn_x_ratio_sq\n", + "pyht_beam.y *= epsn_y_ratio_sq\n", + "pyht_beam.yp *= epsn_y_ratio_sq\n", + "pyht_beam.z *= sigz_ratio\n", + "pyht_beam.dp *= dp_ratio" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import h5py as hp\n", + "def read_all_data(bfile,data):\n", + " bunchdata = hp.File(bfile + '.h5')\n", + "\n", + " # Bunchdata\n", + " bdata = bunchdata['Bunch']\n", + " \n", + " n_turns = len(bdata[data])\n", + " _ = np.empty(n_turns)\n", + " for key in bdata.keys():\n", + " _[:] = bdata[key][:]\n", + " \n", + " bunchdata.close()\n", + "\n", + "def read_n_plot_data(bfile,nturns,data,factor):\n", + " bunchdata = hp.File(bfile + '.h5')\n", + "\n", + " fig = plt.figure(figsize=(16, 16))\n", + " ax1 = fig.add_subplot(311)\n", + "\n", + " ax1.plot(bunchdata['Bunch'][data][:nturns]*factor)\n", + " ax1.set_xlabel('turns', fontsize=20)\n", + " ax1.set_ylabel(f'{data} of bunch', fontsize=20)\n", + " ax1.set_title('BunchMonitor', fontsize=20)\n", + " \n", + " plt.tight_layout()\n", + " \n", + " bunchdata.close()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# Monitors\n", + "bunch_filename = 'bunch_mon'\n", + "particle_filename = 'particle_mon'\n", + "\n", + "bunch_monitor = BunchMonitor(\n", + " filename=bunch_filename, n_steps=nturns, \n", + " parameters_dict={'Q_x': qx_act}, write_buffer_every=2\n", + ")\n", + "particle_monitor = ParticleMonitor(\n", + " filename=particle_filename, stride=10, \n", + " parameters_dict={'Q_x': qx_act}\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "os.remove(bunch_filename + '.h5')\n", + "os.remove(particle_filename + '.h5part')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "map_ = list(trans_map_smooth)+[longitudinal_map]\n", + "dt=twiss.summary.length/(pyht_beam.beta*c)\n", + "rxl=[]\n", + "ryl=[]\n", + "rsl=[]\n", + "time=[]\n", + "\n", + "for i in range(nturns):\n", + " t=i+1\n", + " for m_ in map_:\n", + " m_.track(pyht_beam)\n", + " #Update beam dictionary \n", + " beam_dict={'Z':Q,\n", + " 'mass':mass_MeV,\n", + " 'KE':Ekin/1e+6,\n", + " 'emit_x':pyht_beam.epsn_x()/pyht_beam.betagamma,\n", + " 'emit_y':pyht_beam.epsn_y()/pyht_beam.betagamma,\n", + " 'dp_p':pyht_beam.sigma_dp(),\n", + " 'N_ptcl':intensity,\n", + " 'sigma_s0':pyht_beam.sigma_z() \n", + " } \n", + " #Calculate rates\n", + " IBS=IBS_Martini_rates(twiss,beam_dict,True,60,60,100,12.5,1)\n", + " print(f'Time = {dt*t} Rates X {IBS.rsx} # Rates Y {IBS.rsy} # Rates Z {IBS.rss} #')\n", + " \n", + " growth_rx=IBS.rsx\n", + " growth_ry=IBS.rsy\n", + " growth_rs=IBS.rss\n", + " \n", + " rxl.append(growth_rx)\n", + " ryl.append(growth_ry)\n", + " rsl.append(growth_rs)\n", + " time.append(dt*t)\n", + " \n", + " #Apply IBS kick\n", + " IBSMap=IBS_kick(growth_rx,growth_ry,growth_rs, dt ,Dx=D_x_inj,Dy=0)\n", + " IBSMap.track(pyht_beam)\n", + " \n", + " bunch_monitor.dump(pyht_beam) \n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "with hp.File(bunch_filename + '.h5') as bunchdata:\n", + " print(bunchdata['Bunch']['epsn_x'][0]/pyht_beam.betagamma)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "with hp.File(bunch_filename + '.h5') as bunchdata:\n", + "\n", + " fig = plt.figure(figsize=(8,5))\n", + "\n", + " data='ryl'\n", + " plt.plot(rsl)\n", + " plt.xlabel('turns', fontsize=20)\n", + " plt.ylabel(f'{data} of bunch', fontsize=20)\n", + " plt.title('BunchMonitor', fontsize=20)\n", + " plt.twinx()\n", + " plt.plot((bunchdata['Bunch']['sigma_z'][0:nturns]))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "with hp.File(bunch_filename + '.h5') as bunchdata:\n", + " \n", + " exb=bunchdata['Bunch']['epsn_x'][0:nturns]\n", + " eyb=bunchdata['Bunch']['epsn_y'][0:nturns]\n", + " sigzb=bunchdata['Bunch']['sigma_z'][0:nturns]\n", + " b=list(zip( time,exb,eyb,dpl,sigzb,rxl,ryl,rsl)) \n", + " np.savetxt(\"PyH_everything after.csv\", np.array(b), delimiter=\"\\t\") " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure(figsize=(16, 16))\n", + "ax1 = fig.add_subplot(311)\n", + "\n", + "data='ryl'\n", + "ax1.plot(rxl)\n", + "\n", + "#ax1.plot(bunchdata['Bunch'][data][:nturns]*1)\n", + "ax1.set_xlabel('turns', fontsize=20)\n", + "ax1.set_ylabel(f'{data} of bunch', fontsize=20)\n", + "ax1.set_title('BunchMonitor', fontsize=20)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure(figsize=(16, 16))\n", + "ax1 = fig.add_subplot(311)\n", + "\n", + "data='ryl'\n", + "ax1.plot(ryl)\n", + "#ax1.plot(bunchdata['Bunch'][data][:nturns]*1)\n", + "ax1.set_xlabel('turns', fontsize=20)\n", + "ax1.set_ylabel(f'{data} of bunch', fontsize=20)\n", + "ax1.set_title('BunchMonitor', fontsize=20)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.8.3" + }, + "varInspector": { + "cols": { + "lenName": 16, + "lenType": 16, + "lenVar": 40 + }, + "kernels_config": { + "python": { + "delete_cmd_postfix": "", + "delete_cmd_prefix": "del ", + "library": "var_list.py", + "varRefreshCmd": "print(var_dic_list())" + }, + "r": { + "delete_cmd_postfix": ") ", + "delete_cmd_prefix": "rm(", + "library": "var_list.r", + "varRefreshCmd": "cat(var_dic_list()) " + } + }, + "types_to_exclude": [ + "module", + "function", + "builtin_function_or_method", + "instance", + "_Feature" + ], + "window_display": false + } + }, + "nbformat": 4, + "nbformat_minor": 1 +} diff --git a/simulation_notebooks/IBS_mod/ELENA_elements.def b/simulation_notebooks/IBS_mod/ELENA_elements.def new file mode 100644 index 0000000..e874adb --- /dev/null +++ b/simulation_notebooks/IBS_mod/ELENA_elements.def @@ -0,0 +1,182 @@ +! +! ELENA ring elements. +! +! +! Changelog: +! - 5/06/2018: started with Pavel''s optics. + + + +/****************************************************************************************************** +! Parameters used for building dipoles. +! The bending radius of ELENA magnet is equal 0.927m, and its length is equal to rho*phi +******************************************************************************************************/ + LBM := 0.927*Pi/3; + BME := PI*16.45/180.; + +/**************************************************************************************************** +! The value of bending magnet edge angle has been modified after comparison of tracking particles +! through the manufactured magnet with OPERA program and tracking of particles through the designed +! bending magnet which is described with MADX. +! +! To compare these two magnets, many particles initially have been placed on ellipse well distanced +! from magnet to take into account extended fields. To fit ellipse inclination after tracking with MADX +! (Olav Berrig) to that after tracking with OPERA (Alexander Aloev,visitor from Moscow State University), +! edge angle has been varied, and then corrected from 17 degrees (which was design value) to 16.45 degrees. +! +! In addition to that, tracking showed visible distortion of the ellipse shape, which was successfully +! corrected by introducing sextupolar component on both edge of magnet. It was resulted from the +! nonlinearities od edge field, which was derived analytically by Pavel Belochitskii) +******************************************************************************************************/ + BMH: SBEND, L=LBM, ANGLE=PI/3,E1=BME, E2=BME, FINT=0.424, HGAP=0.076/2,apertype=ellipse,aperture={0.03,0.03}; + +/****************************************************************************************************** +! Sometimes one needs to cut magnet in 2 parts, for example for better plots of Twiss functiions, +! or for better calculations of lifetimes, or better knowledge of average values etc. +! In this case it worth to cut bending magnets in pieces +******************************************************************************************************/ + BMH1: SBEND, L=LBM/2,ANGLE=PI/6,E1=BME, E2=0 , FINT=0.424, HGAP=0.076/2, FINTX=0, apertype=ellipse,aperture={0.03,0.03}; + BMH2: SBEND, L=LBM/2,ANGLE=PI/6,E1=0 , E2=BME, FINT=0, HGAP=0.076/2, FINTX=0.424, apertype=ellipse,aperture={0.03,0.03}; + +/****************************************************************************************************** +! thin quadrupole length can be installed on both sides of ELENA bending magnet to modify edge angle +******************************************************************************************************/ + BME_edge_angle_error: MULTIPOLE, KNL={0,PI*16.45/180*0.001}; + +/****************************************************************************************************** +! A multipole definition of SBEND allowing to add multipole errors +******************************************************************************************************/ + +!OPTION, THIN_FOC=true; + + edgeBMH1: DIPEDGE,H=(PI/3)/(LBM), E1=BME , FINT=0.424, HGAP=0.076/2; + MBMH1: multipole,LRAD=LBM/2, KNL={PI/6},apertype=ellipse,aperture={0.03,0.03}; + MBMH2: multipole,LRAD=LBM/2, KNL={PI/6},apertype=ellipse,aperture={0.03,0.03}; + edgeBMH2: DIPEDGE,H=(PI/3)/(LBM), E1=BME , FINT=0.424, HGAP=0.076/2; + +/****************************************************************************************************** +! quadrupoles are divided in two equal parts +******************************************************************************************************/ + Q1H: QUADRUPOLE, L=0.125, K1:=KQ1,apertype=ellipse,aperture={0.03,0.03}; + Q2H: QUADRUPOLE, L=0.125, K1:=KQ2,apertype=ellipse,aperture={0.03,0.03}; + Q3H: QUADRUPOLE, L=0.125, K1:=KQ3,apertype=ellipse,aperture={0.03,0.03}; + +/****************************************************************************************************** +! skew quadrupoles used for the residual coupling compensation +******************************************************************************************************/ + QSK1: QUADRUPOLE, L=0.15, K1S:=KSQ1,apertype=ellipse,aperture={0.03,0.03}; + QSK2: QUADRUPOLE, L=0.15, K1S:=KSQ2,apertype=ellipse,aperture={0.03,0.03}; + +/***************************************************************************************************** +! These quadrupoles could be used like trims in case asymmetry in ELENA will be needed +! (for example, to put quadrupole trims in electron cooler section) +*****************************************************************************************************/ +! Q1T: QUADRUPOLE, L=0.125, K1:=KQ1 + dKQ1; +! Q2T: QUADRUPOLE, L=0.125, K1:=KQ2 + dKQ2; +! Q3T: QUADRUPOLE, L=0.125, K1:=KQ3 + dKQ3; + +/*************************************************************************************************** +! These quadrupoles will be used to take into account tilt of normal quadrupoles +****************************************************************************************************/ +! TQ1: multipole, KNL:={0,K1TQ1}; +! TQ2: multipole, KNL:={0,K1TQ2}; +! TQ3: multipole, KNL:={0,K1TQ3}; + + +/****************************************************************************************************** +! Electron cooler solenoid is called MSOL, two solenoids SCOMP used for coupling compensation. +! Cooler solenoid may be divided in pieces to take into account effect +! of space charge of antiproton beam (Laslett tune shift) +******************************************************************************************************/ + MSOL: SOLENOID, L=0.650, KS:=KMSOL,apertype=ellipse,aperture={0.03,0.03}; + !!Disable Ecool + MSOL1: SOLENOID, L=0.650-((0.763 + 0.1016)/2), KS:=KMSOL,apertype=ellipse,aperture={0.03,0.03}; + !!MSOL1: drift, L=0.650-((0.763 + 0.1016)/2); + !!Disable Ecool + MSOL2: SOLENOID, L=(0.763 + 0.1016)/2 , KS:=KMSOL,apertype=ellipse,aperture={0.03,0.03}; + !!MSOL2: drift, L=(0.763 + 0.1016)/2; + SCOMP: SOLENOID, L=0.360, KS:=KSCOMP,apertype=ellipse,aperture={0.03,0.03}; + +/****************************************************************************************************** +! kicks from toroid coils of electron cooler taken into account +******************************************************************************************************/ +!!Disable Ecool + +TOR1: KICKER, L=0., HKICK:= Tkick, VKICK:=0,apertype=ellipse,aperture={0.03,0.03}; +TOR2: KICKER, L=0., HKICK:=-Tkick, VKICK:=0,apertype=ellipse,aperture={0.03,0.03}; + +!!TOR1: drift, L=0.0; +!!TOR2: drift, L=0.0; + + +/****************************************************************************************************** +! sextupoles used for chromaticity correction +******************************************************************************************************/ + SF: SEXTUPOLE, L=0.15, K2:=KSF,apertype=ellipse,aperture={0.03,0.03}; + SD: SEXTUPOLE, L=0.15, K2:=KSD,apertype=ellipse,aperture={0.03,0.03}; + + +/****************************************************************************************************** +! dipoles for orbit correction +******************************************************************************************************/ + DHV: KICKER, L=0., HKICK:=0., VKICK:=0.,apertype=ellipse,aperture={0.03,0.03}; + +/****************************************************************************************************** +! Kyoto-style correctors around electron cooler +******************************************************************************************************/ + DHVEC1: KICKER, L=0., HKICK:=0., VKICK:=0.,apertype=ellipse,aperture={0.03,0.03}; + DHVEC2: KICKER, L=0., HKICK:=0., VKICK:=0.,apertype=ellipse,aperture={0.03,0.03}; + + + +/****************************************************************************************************** +! injection kicker +******************************************************************************************************/ + KFI: HKICKER, l=0.432, KICK:=0,apertype=ellipse,aperture={0.03,0.03}; + + +/****************************************************************************************************** +! extraction kickers +! In fact, both are electrostatic elements, and for simulations have to be treated in a proper way +******************************************************************************************************/ + DFH1: HKICKER,L=0.5, KICK:=0.,apertype=ellipse,aperture={0.03,0.03}; +! Increased kicker to LNE50 length from 0.5 to 0.6 as ABT transferline optics files +! davide Sep 2018 + DFH2: HKICKER,L=0.6, KICK:=0.,apertype=ellipse,aperture={0.03,0.03}; + +/****************************************************************************************************** +! RF cavity +******************************************************************************************************/ + RFC: RFCAVITY, L=0.3, VOLT:=RFvoltage,apertype=ellipse,aperture={0.03,0.03}; + + + +/*************************************************************************************************** +! Instrumentation +****************************************************************************************************/ + BPM: MONITOR, L=0,apertype=ellipse,aperture={0.03,0.03}; + USL: INSTRUMENT, L=0.425,apertype=ellipse,aperture={0.03,0.03}; + SCRAP: INSTRUMENT, L=0.550,apertype=ellipse,aperture={0.03,0.03}; + VACP: INSTRUMENT, L=0.,apertype=ellipse,aperture={0.03,0.03}; + !BTV: INSTRUMENT, L=0.392/2,apertype=ellipse,aperture={0.03,0.03}; ! should probably be deleted - davide Jun2018 + BTV: INSTRUMENT, L=0.,apertype=ellipse,aperture={0.03,0.03}; + Qmeter: INSTRUMENT, L=0.350,apertype=ellipse,aperture={0.03,0.03}; + + +/*************************************************************************************************** +! Markers +****************************************************************************************************/ + MSS: MARKER; + EBEAM: MARKER; + MECOOL: MARKER; + KFISTART: MARKER; + SMI: MARKER; + SMIENDM: MARKER; + SMIEND: MARKER; + BTV_MID: MARKER; + M10: MARKER; + M20: MARKER; + MAXBETX: MARKER; + + +RETURN; diff --git a/simulation_notebooks/IBS_mod/ELENA_ring_mod.seq b/simulation_notebooks/IBS_mod/ELENA_ring_mod.seq new file mode 100644 index 0000000..2113c95 --- /dev/null +++ b/simulation_notebooks/IBS_mod/ELENA_ring_mod.seq @@ -0,0 +1,197 @@ +! +! ELENA ring sequence +! +! +! Changelog: +! 5/06/2018: +! - started with Pavel''s optics. +! - cleaned up some commented lines and added some markers +! - redefined randomly some names of the e-cool solenoids +! - TODO: cleanup useless lines. re-arrange markers "badly" used here and there. + + + +TITLE, "ELENA ring"; + +! it seems well compatible with https://edms.cern.ch/ui/file/1252808/AJ/ad_lm___0071-vAJ.pdf! +! note, LBM = 0.927*pi/3 = 0.9708. This 0.97 here and there is an approximation, I guess... +SECTION1: SEQUENCE, REFER=CENTRE, L=4.50 + 0.97; + LNR.MCCAY.0105: DHV, at = 0.3206; + LNR.BPMEB.0110: BPM, at = 0.3206; + LNR.MMIDINJ: MARKER, at = 2.25; ! basically the middle of injection section + LNR.MBMIE.0115: MARKER, at = 0.15035 + 2.27925; ! placed inside of injection septum + LNR.BTVPA.0117: BTV, at = 2.9984 - 0.392/4; + LNR.BTVPA.0118: BTV, at = 2.9984 + 0.392/4; + LNR.MKKFH.0120: KFI, at = 3.5906; + LNR.BPMEB.0125: BPM, at = (0.97-LBM)/2 + 4.179 ; + LNR.MCCAY.0130: DHV, at = (0.97-LBM)/2 + 4.179; + LNR.MBHEK.0135a: BMH1, at = 4.5 - (LBM - 0.97) + LBM/2 - LBM/4; + LNR.MBHEK.0135: MARKER, at = 4.5 - (LBM - 0.97) + LBM/2; +! LNR.MBHEK.0135: BMH, at = 4.5 - (LBM - 0.97) + LBM/2; + LNR.MBHEK.0135b: BMH2, at = 4.5 - (LBM - 0.97) + LBM/2 + LBM/4; +ENDSEQUENCE; + +SECTION2: SEQUENCE, REFER=CENTRE, L=3.8964-(LBM-0.97) + LBM; + LNR.MXNAD.0201: SF, at = (0.97 - LBM)/2 + 0.275; + LNR.MQNLG.0205: Q1H, at = (0.97 - LBM)/2 + 0.625 - 0.125/2; + LNR.MQNLGm.0205: MARKER, at = (0.97 - LBM)/2 + 0.625; + LNR.MQNLG.0206: Q1H, at = (0.97 - LBM)/2 + 0.625 + 0.125/2; + LNR.MQNLG.0210: Q2H, at = (0.97 - LBM)/2 + 1.0982 - 0.125/2; + LNR.BPMEA.0215: BPM, at = (0.97 - LBM)/2 + 1.0982; + LNR.MQNLG.0211: Q2H, at = (0.97 - LBM)/2 + 1.0982 + 0.125/2; + LNR.MXNAD.0220: SD, at = (0.97 - LBM)/2 + 1.4132 + 0.150/2; + LNR.APULA.0225: USL, at = (0.97 - LBM)/2 + 1.88195; + LNR.APULA.0227: USL, at = (0.97 - LBM)/2 + 2.39945; + LNR.MCCAY.0230: DHV, at = (0.97 - LBM)/2 + 3.1884 - 0.209/2; + LNR.MQNLG.0235: Q3H, at = (0.97 - LBM)/2 + 3.4714 - 0.125/2; + LNR.MQNLG.0235m: MARKER, at = (0.97 - LBM)/2 + 3.4714; + LNR.BPMEA.0240: BPM, at = (0.97 - LBM)/2 + 3.4714; + LNR.MQNLG.0236: Q3H, at = (0.97 - LBM)/2 + 3.4714 + 0.125/2; + LNR.MBHEK.0245a: BMH1, at = 3.8964 - (LBM - 0.97) + LBM/2 - LBM/4; + LNR.MBHEK.0245: MARKER, at = 3.8964 - (LBM - 0.97) + LBM/2; +! LNR.MBHEK.0245: BMH, at = 3.8964 - (LBM - 0.97) + LBM/2; + LNR.MBHEK.0245b: BMH2, at = 3.8964 - (LBM - 0.97) + LBM/2 + LBM/4; +ENDSEQUENCE; + +SECTION3: SEQUENCE, REFER=CENTRE, L=3.8964-(LBM-0.97) + LBM; + LNR.MQNLG.0305: Q3H, at = (0.97 - LBM)/2 + 0.425 - 0.125/2; + LNR.MQNLG.0305m: MARKER, at = (0.97 - LBM)/2 + 0.425; + LNR.MQNLG.0306: Q3H, at = (0.97 - LBM)/2 + 0.425 + 0.125/2; +/*** Trim quads for adjustment of beta functions around cooler, not implemented ***/ +! LNR.MQNLG.0305a: Q3T, at = (0.97 - LBM)/2 + 0.425 - 0.125/2; +! LNR.MQNLG.0305b: Q3T, at = (0.97 - LBM)/2 + 0.425 + 0.125/2; + LNE.MAR.5000: MARKER, at = (0.97 - LBM)/2 + 1.115 - 0.3; ! marker where LNE50 optics starts + LNR.ZDFA.0310: DFH2, at = (0.97 - LBM)/2 + 1.115; +! M20, at = (0.97 - LBM)/2 + 2.0086 - 0.375; + LNR.BQKMA.0312: Qmeter, at = (0.97 - LBM)/2 + 2.0086; + LNR.MQNLG.0315: Q2H, at = (0.97 - LBM)/2 + 2.7982 - 0.125/2; + LNR.MQNLG.0315m: MARKER, at = (0.97 - LBM)/2 + 2.7982; + LNR.BPMEA.0320: BPM, at = (0.97 - LBM)/2 + 2.7982; + LNR.MQNLG.0316: Q2H, at = (0.97 - LBM)/2 + 2.7982 + 0.125/2; +/*** Trim quads for adjustment of beta functions around cooler, not implemented ***/ +! LNR.MQNLG.0315a: Q2T, at = (0.97 - LBM)/2 + 2.7982 - 0.125/2; +! LNR.MQNLG.0315b: Q2T, at = (0.97 - LBM)/2 + 2.7982 + 0.125/2; + LNR.MQNLG.0325: Q1H, at = (0.97 - LBM)/2 + 3.2714 - 0.125/2; + LNR.MQNLG.0326: Q1H, at = (0.97 - LBM)/2 + 3.2714 + 0.125/2; +/*** Trim quads for adjustment of beta functions around cooler, not implemented ***/ +! LNR.MQNLG.0325a: Q1T, at = (0.97 - LBM)/2 + 3.2714 - 0.125/2; +! LNR.MQNLG.0325b: Q1T, at = (0.97 - LBM)/2 + 3.2714 + 0.125/2; + LNR.MCCAY.0330: DHV, at = (0.97 - LBM)/2 + 3.6059; + LNR.MBHEK.0335a: BMH1, at = 3.8964 - (LBM - 0.97) + LBM/2 - LBM/4; + LNR.MBHEK.0335: MARKER, at = 3.8964 - (LBM - 0.97) + LBM/2; +! LNR.MBHEK.0335: BMH, at = 3.8964 - (LBM - 0.97) + LBM/2; + LNR.MBHEK.0335b: BMH2, at = 3.8964 - (LBM - 0.97) + LBM/2 + LBM/4; +ENDSEQUENCE; + +SECTION4: SEQUENCE, REFER=CENTRE, L=4.5-(LBM-0.97) + LBM; +! LNR.MCCAY.0405: DHV41, at = (0.97 - LBM)/2 + 0.321; +! LNR.MCCAY.0405: DHV, at = (0.97 - LBM)/2 + 0.321 - 0.209/4; + LNR.MCCAY.0405: DHV, at = (0.97 - LBM)/2 + 2.25 - 1.92125; +! LNR.BPMEB.0407: BPM, at = (0.97 - LBM)/2 + 0.321; + LNR.BPMEB.0407: BPM, at = (0.97 - LBM)/2 + 2.25 - 1.92125; + LNR.MLNAF.0410: SCOMP, at = (0.97-LBM)/2 + 0.725; +! LNR.MCCAY.0420: DHVEC1, at = (0.97 - LBM)/2 + 2.25 - 1.01382; + LNR.MCCAY.0420: DHV, at = (0.97 - LBM)/2 + 2.25 - 1.01382; +!!!!! start E-COOLER !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +/*********** Kyoto-style correctors integrated into electron cooler ***********/ + LNR.MCVEA.0425: DHVEC2, at = (0.97 - LBM)/2 + 2.25 - 0.8698 ; + LNR.ECTOR.0427xx: TOR1, at = (0.97 - LBM)/2 + 2.25 - 0.66; + LNR.ECSOL.0427xx: MSOL1, at = (0.97 - LBM)/2 + 2.25 - (0.65/2+(0.763 + 0.1016)/4); + LNR.BPMEC.0428: BPM, at = (0.97 - LBM)/2 + 2.25 - (0.763 + 0.1016)/2; +!! MSOL, at = (0.97 - LBM)/2 + 2.25 - 0.65/2; +!! EBEAM, at = (0.97 - LBM)/2 + 2.25 - 0.65; +!! MSOL2, at = (0.97 - LBM)/2 + 2.25 - (0.753 + 0.1016)/4; + LNR.ECSOL.0427yy: MSOL2, at = (0.97 - LBM)/2 + 2.25 - (0.763 + 0.1016)/4; + EBEAM, at = (0.97 - LBM)/2 + 2.25; + MECOOL, at = (0.97 - LBM)/2 + 2.25; + EBEAM, at = (0.97 - LBM)/2 + 2.25; +!! MSOL, at =(0.97 - LBM)/2 + 2.25 + 0.325; + LNR.ECSOL.0433yy: MSOL2, at = (0.97 - LBM)/2 + 2.25 + (0.763 + 0.1016)/4; +!! EBEAM, at = (0.97 - LBM)/2 + 2.25 + 0.65; + LNR.BPMEC.0432: BPM, at = (0.97 - LBM)/2 + 2.25 + (0.763 + 0.1016)/2; + LNR.ECSOL.0433xx: MSOL1, at = (0.97 - LBM)/2 + 2.25 + (0.65/2+(0.763 + 0.1016)/4); + LNR.ECTOR.0433xx: TOR2, at = (0.97 - LBM)/2 + 2.25 + 0.66; +/*********** Kyoto-style correctors integrated into electron cooler ***********/ + LNR.MCVEA.0435: DHVEC2, at = (0.97 - LBM)/2 + 2.25 + 0.8698; +!!!!! end E-COOLER !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + LNR.MCCAY.0440: DHVEC1, at = (0.97 - LBM)/2 + 2.25 + 1.01382; + LNR.MLNAF.0450: SCOMP, at = (0.97 - LBM)/2 + 4.392 - 0.617; +! LNR.MCCAY.0460: DHV, at = (0.97 - LBM)/2 + 4.392 - 0.213 - 0.209/4; + LNR.MCCAY.0460: DHV, at = (0.97 - LBM)/2 + 2.25 + 1.92125; +! LNR.BPMEB.0457: BPM, at = (0.97 - LBM)/2 + 4.392 - 0.213; + LNR.BPMEB.0457: BPM, at = (0.97 - LBM)/2 + 2.25 + 1.92125; + LNR.MBHEK.0470a: BMH1, at = 4.5 - (LBM - 0.97) + LBM/2 - LBM/4; + LNR.MBHEK.0470: MARKER, at = 4.5 - (LBM - 0.97) + LBM/2; +! LNR.MBHEK.0470: BMH, at = 4.5 - (LBM - 0.97) + LBM/2; + LNR.MBHEK.0470b: BMH2, at = 4.5 - (LBM - 0.97) + LBM/2 + LBM/4; +ENDSEQUENCE; + +SECTION5: SEQUENCE, REFER=CENTRE, L=3.8964 - (LBM - 0.97) + LBM; + LNR.MXNAD.0505: SF, at = (0.97 - LBM)/2 + 0.200 + 0.150/2; +/*** Trim quads for adjustment of beta functions around cooler, not implemented ***/ +! Q1T, at = (0.97 - LBM)/2 + 0.625 - 0.125/2; +! Q1T, at = (0.97 - LBM)/2 + 0.625 + 0.125/2; + LNR.MQLNG.0510: Q1H, at = (0.97 - LBM)/2 + 0.625 - 0.125/2; + LNR.MQLNG.0510m: MARKER, at = (0.97 - LBM)/2 + 0.625; + LNR.MQLNG.0511: Q1H, at = (0.97 - LBM)/2 + 0.625 + 0.125/2; +/*** Trim quads for adjustment of beta functions around cooler, not implemented ***/ +! Q2T, at = (0.97 - LBM)/2 + 1.0982 - 0.125/2; +! Q2T, at = (0.97 - LBM)/2 + 1.0982 + 0.125/2; + LNR.MQNLG.0515: Q2H, at = (0.97 - LBM)/2 + 1.0982 - 0.125/2; + LNR.MQNLG.0515m: MARKER, at = (0.97 - LBM)/2 + 1.0982; + LNR.BPMEA.0520: BPM, at = (0.97 - LBM)/2 + 1.0982; + LNR.MQNLG.0516: Q2H, at = (0.97 - LBM)/2 + 1.0982 + 0.125/2; + LNR.MXNAD.0525: SD, at = (0.97 - LBM)/2 + 1.4132 + 0.150/2; + LNR.ACWO2.0530: RFC, at = (0.97 - LBM)/2 + 1.7972; +! VACP, at = (0.97-LBM)/2 + 2.0449; + LNR.MCCAY.0535: DHV, at = (0.97 - LBM)/2 + 2.2972; + LNR.MQSAB.0540: QSK1, at = (0.97 - LBM)/2 + 2.5822; +! SCRAP, at = (0.97-LBM)/2 + 2.9114; + LNR.BSSHV.0545: SCRAP, at = (0.97 - LBM)/2 + 2.98015; +/*** Trim quads for adjustment of beta functions around cooler, not implemented ***/ +! LNR.MQNLG.0550a: Q3T, at = (0.97 - LBM)/2 + 3.4714 - 0.125/2; +! LNR.MQNLG.0550b: Q3T, at = (0.97 - LBM)/2 + 3.4714 + 0.125/2; + LNR.MQNLG.0550: Q3H, at = (0.97 - LBM)/2 + 3.4714 - 0.125/2; + LNR.MQNLG.0550m: MARKER, at = (0.97 - LBM)/2 + 3.4714; + LNR.BPMEA.0555: BPM, at = (0.97 - LBM)/2 + 3.4714; + LNR.MQNLG.0551: Q3H, at = (0.97 - LBM)/2 + 3.4714 + 0.125/2; + LNR.MBHEK.0560a: BMH1, at = 3.8964 - (LBM - 0.97) + LBM/2 - LBM/4; + LNR.MBHEK.0560: MARKER, at = 3.8964 - (LBM - 0.97) + LBM/2; +! LNR.MBHEK.0560: BMH, at = 3.8964 - (LBM - 0.97) + LBM/2; + LNR.MBHEK.0560b: BMH2, at = 3.8964 - (LBM - 0.97) + LBM/2 + LBM/4; + ENDSEQUENCE; + +SECTION6: SEQUENCE, REFER=CENTRE, L=3.8964 - (LBM - 0.97) + LBM; + LNR.MQNLG.0605: Q3H, at = (0.97 - LBM)/2 + 0.425 - 0.125/2; + LNR.MQNLG.0605m: MARKER, at = (0.97 - LBM)/2 + 0.425; + LNR.MQNLG.0606: Q3H, at = (0.97 - LBM)/2 + 0.425 + 0.125/2; + LNR.ZDFHL.0610: DFH1, at = (0.97 - LBM)/2 + 1.115; +! M20, at = (0.97 - LBM)/2 + 1.115 + 0.45; + LNR.BQKMA.0612: Qmeter, at = (0.97 - LBM)/2 + 2.0086; + LNR.MCCAY.0615: DHV, at = (0.97 - LBM)/2 + 2.286 + 0.209/2; + LNR.MQLNG.0620: Q2H, at = (0.97 - LBM)/2 + 2.7982 - 0.125/2; + LNR.BPMEA.0625: BPM, at = (0.97 - LBM)/2 + 2.7982; + LNR.MQLNG.0626: Q2H, at = (0.97 - LBM)/2 + 2.7982 + 0.125/2; + LNR.MQLNG.0630: Q1H, at = (0.97 - LBM)/2 + 3.2714 - 0.125/2; + LNR.MQLNG.0631: Q1H, at = (0.97 - LBM)/2 + 3.2714 + 0.125/2; + LNR.MQSAB.0635: QSK2, at = (0.97 - LBM)/2 + 3.591; + LNR.MBHEK.0640a: BMH1, at = 3.8964 - (LBM - 0.97) + LBM/2 - LBM/4; + LNR.MBHEK.0640: MARKER, at = 3.8964 - (LBM - 0.97) + LBM/2; +! LNR.MBHEK.0640: BMH, at = 3.8964 - (LBM - 0.97) + LBM/2; + LNR.MBHEK.0640b: BMH2, at = 3.8964 - (LBM - 0.97) + LBM/2 + LBM/4; +ENDSEQUENCE; + +! Davide - April 2019: +! added begin/end markers +ELENA: SEQUENCE,REFER=ENTRY,L= 2*4.5 - 6*(LBM - 0.97) + 6*LBM + 4*3.8964; + LNR.BEGIN: MARKER, at = 0; + SECTION1, at = 0; + SECTION2, at = 4.5 - (LBM - 0.97) + LBM; + SECTION3, at = 4.5 - 2*(LBM - 0.97) + 2*LBM + 3.8964; + SECTION4, at = 4.5 - 3*(LBM - 0.97) + 3*LBM + 2*3.8964; + SECTION5, at = 2*4.5 - 4*(LBM - 0.97) + 4*LBM + 2*3.8964; + SECTION6, at = 2*4.5 - 5*(LBM - 0.97) + 5*LBM + 3*3.8964; + LNR.END: MARKER, at = 2*4.5 - 6*(LBM - 0.97) + 6*LBM + 4*3.8964; +ENDSEQUENCE; + +RETURN; diff --git a/simulation_notebooks/IBS_mod/PyHEAD_IBS.py b/simulation_notebooks/IBS_mod/PyHEAD_IBS.py new file mode 100755 index 0000000..fd75d7a --- /dev/null +++ b/simulation_notebooks/IBS_mod/PyHEAD_IBS.py @@ -0,0 +1,627 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Apr 28 01:38:50 2021 + +@author: Volodymyr +""" +from math import sqrt, pi,sin,cos +from scipy.constants import e, m_p, c, epsilon_0 +from scipy.constants import physical_constants +import numpy as np +import copy +import random + +class IBS_Martini_rates(object): + + def __init__(self, twiss,beam_dict,bunch,nu_x,nu_y,nu_z,log_c,fact): + + + self.beam_dict=beam_dict + self.twiss=twiss + + self.nu_x=nu_x + self.nv_y=nu_y + self.nu_z=nu_z + + self.log_c=log_c + + + #Ion beam parameters + + self.Z = self.beam_dict['Z'] + self.mass =self.beam_dict['mass'] + self.KE=self.beam_dict['KE'] + + self.k_ke = 8.9875517873681764E9; #Coulomb constant + self.emit_x=self.beam_dict['emit_x'] + self.emit_y=self.beam_dict['emit_y'] + self.dp_p=self.beam_dict['dp_p'] + self.N_ptcl=self.beam_dict['N_ptcl'] + self.sigma_s0=self.beam_dict['sigma_s0'] + + gamma = 1+self.KE/self.mass; + beta = sqrt(gamma*gamma-1)/gamma; + + self.gamma=gamma + self.beta=beta + + self.r_i = self.k_ke*self.Z**2*e*1e-6/self.mass; #Classical ion radius + self.bunched = bunch; + + self.energy_spread = beta*beta*self.dp_p; + p0= gamma*self.mass*1e6*e*beta/c; + + ##Lattice functions + s_v=self.twiss['s'] + betx_v=self.twiss['betx'] + alfx_v=self.twiss['alfx'] + mux_v=self.twiss['mux'] + dx_v=self.twiss['dxst']*fact + dpx_v=self.twiss['dpxst']*fact + bety_v=self.twiss['bety'] + alfy_v=self.twiss['alfy'] + muy_v=self.twiss['muy'] + dy_v=self.twiss['dyst']*fact + dpy_v=self.twiss['dpyst']*fact + + l_element_v=[0] + + n_element=len(betx_v) + + for n in range(n_element-1): + l_element_v.append(s_v[n+1]-s_v[n]) + + circ=twiss.summary.length + + f0 = beta*c/circ; #revolution frequency. + w0= 2*pi*f0; #angular frequency. + + ##Containers for 3D grid integration + + #Containers for 3D grid integration + class TrigonometryStorageUV: + sin_u2_cos_v2=g1=g2_1=g2_2=0 + + def __init__(self,sin_u2_cos_v2,g1,g2_1,g2_2): + self.sin_u2_cos_v2 = sin_u2_cos_v2 + self.g1 = g1 + self.g2_1 = g2_1 + self.g2_2 = g2_2 + + class TrigonometryStorageV: + sin_v=cos_v=0 + + def __init__(self,sin_v,cos_v ): + self.sin_v = sin_v + self.cos_v = cos_v + + class TrigonometryStorageU: + sin_u=sin_u2=cos_u2=g3=0 + storage_uv=[] + + def __init__(self,sin_u,sin_u2,cos_u2,g3,storage_uv): + self.sin_u = sin_u + self.sin_u2 = sin_u2 + self.cos_u2 = cos_u2 + self.g3 = g3 + self.storage_uv = storage_uv + + class OpticalStorage: + a=b2=c2=d2=dtld=k1=k2=k3=0 + + def __init__(self,a,b2,c2,d2,dtld,k1,k2,k3): + self.a = a + self.b2 = b2 + self.c2 = c2 + self.d2 = d2 + self.dtld = dtld + self.k1 = k1 + self.k2 = k2 + self.k3 = k3 + + rx_ibs=ry_ibs=rs_ibs=0 + + ###Solver section### + + sigma_xbet=[None] * n_element + sigma_y=[None] * n_element + alfx2=[None] * n_element + sigma_xp=[None] * n_element + alfy2=[None] * n_element + sigma_yp=[None] * n_element + + #Bunch properties at each element + for i in range(n_element): + sigma_xbet[i]=sqrt(betx_v[i]*self.emit_x); + sigma_y[i]=sqrt(bety_v[i]*self.emit_y); + alfx2[i]=alfx_v[i]**2; + sigma_xp[i]=sqrt((1+alfx2[i])*self.emit_x/betx_v[i]); + alfy2[i]=alfy_v[i]**2; + sigma_yp[i]=sqrt((1+alfy2[i])*self.emit_y/bety_v[i]); + + #abcdk + os_all=[None] * n_element + + for i in range(n_element): + betx = betx_v[i]; + alfx = alfx_v[i]; + dx = dx_v[i]; + dpx = dpx_v[i]; + alfy = alfy_v[i]; + d_tld = alfx*dx+betx*dpx; + d_tld = alfx*dx+betx*dpx; + sigma_x = sqrt(sigma_xbet[i]*sigma_xbet[i]+dx*dx*self.dp_p*self.dp_p); + sigma_tmp = self.dp_p*sigma_xbet[i]/(gamma*sigma_x); + q = 2*beta*gamma*sqrt(sigma_y[i]/self.r_i); + a = sigma_tmp*sqrt(1+alfx*alfx)/sigma_xp[i]; + b2 = sigma_tmp*sqrt(1+alfy*alfy)/sigma_yp[i]; + b2 = b2**2; + c2 = q*sigma_tmp; + c2 = c2**2; + d2 = self.dp_p*dx/sigma_x; + d2 = d2**2; + dtld = self.dp_p*d_tld/sigma_x; + k1 = 1.0 / c2; + k2 = a * a * k1; + k3 = b2 * k1; + os=OpticalStorage(a,b2,c2,d2,dtld,k1,k2,k3); + os_all[i]=copy.deepcopy(os); + os=None + + #Coef f + + storage_u=[] + storage_v=[]; + dv = 2*pi/self.nv_y; + v = -0.5*dv; + + for iv in range(self.nv_y): + v += dv; + storage_v.append(TrigonometryStorageV(sin(v), cos(v))); + + du = pi/nu_x; + u = -0.5*du; + for iu in range(nu_x): + u += du; + sin_u = sin(u); + sin_u2 = sin_u **2; + cos_u2 = 1 - sin_u2; + g3 = 1 - 3 * cos_u2; + storage_uv=[] + for iv in range(self.nv_y): + sin_u2_cos_v2 = sin_u2 * storage_v[iv].cos_v * storage_v[iv].cos_v; + g1 = 1 - 3 * sin_u2_cos_v2; + g2_1 = 1 - 3 * sin_u2 * storage_v[iv].sin_v * storage_v[iv].sin_v; + g2_2 = 6 * sin_u * storage_v[iv].sin_v * storage_v[iv].cos_v; + tempUV = TrigonometryStorageUV(sin_u2_cos_v2, g1, g2_1, g2_2); + storage_uv.append(tempUV); + storage_u.append(TrigonometryStorageU(sin_u, sin_u2, cos_u2, g3, storage_uv)); + + + f1=[None]*n_element + f2=[None]*n_element + f3=[None]*n_element + + if (log_c > 0): + for ie in range(n_element): + duvTimes2Logc = 2*pi*pi/(self.nu_x*self.nv_y) * 2 * self.log_c; + tempf1=tempf2=tempf3=0; + for iu in range(nu_x): + tu = storage_u[iu]; + sum1=sum2=sum3=0; + for iv in range(self.nv_y): + tv = storage_v[iv]; + tuv = tu.storage_uv[iv]; + tmp = os_all[ie].a * tv.sin_v - os_all[ie].dtld * tv.cos_v; + tmp = tmp**2; + inv_int_z = (tuv.sin_u2_cos_v2 + tu.sin_u2 * tmp + os_all[ie].b2 * tu.cos_u2) * os_all[ie].k1; + sum1 += tuv.g1 / inv_int_z; + sum2 += (tuv.g2_1 + tuv.g2_2 * os_all[ie].dtld / os_all[ie].a) / inv_int_z; + sum3 += tu.g3 / inv_int_z; + + tempf1 += tu.sin_u * sum1; + tempf2 += tu.sin_u * sum2; + tempf3 += tu.sin_u * sum3; + + f1[ie]=tempf1 * os_all[ie].k1 * duvTimes2Logc; + f2[ie]=tempf2 * os_all[ie].k2 * duvTimes2Logc; + f3[ie]=tempf3 * os_all[ie].k3 * duvTimes2Logc; + + # if slicing by nz, longer but almost the same + + # f1=[None]*n_element + # f2=[None]*n_element + # f3=[None]*n_element + # for ie in range(10): + + # duv = 2*pi*pi/(nu_x*nv_y); + # tempf1 = tempf2 = tempf3 = 0; + # for iu in range(nu_x): + # tu = storage_u[iu]; + # sum1 = sum2 = sum3 = 0; + # for iv in range(nv_y): + # tv = storage_v[iv]; + # tuv = tu.storage_uv[iv]; + + # tmp = os_all[ie].a * tv.sin_v - os_all[ie].dtld * tv.cos_v; + # tmp = tmp**2; + # d_uv = (tuv.sin_u2_cos_v2 + tu.sin_u2 * tmp + os_all[ie].b2 * tu.cos_u2) * os_all[ie].k1; + # int_z = 0; + # dz = 20/(d_uv*nu_z); + # z = -0.5*dz; + # for iz in range(nu_z): + # z += dz; + # int_z += np.exp(-d_uv*z)*np.log(1+z*z)*dz; + + # sum1 += int_z * tuv.g1; + # sum2 += int_z * (tuv.g2_1 + tuv.g2_2 * os_all[ie].dtld / os_all[ie].a); + # sum3 += int_z * tu.g3; + + # tempf1 += tu.sin_u * sum1; + # tempf2 += tu.sin_u * sum2; + # tempf3 += tu.sin_u * sum3; + + # f1[ie] = tempf1 * os_all[ie].k1* duv; + # f2[ie] = tempf2 * os_all[ie].k2* duv; + # f3[ie] = tempf3 * os_all[ie].k3* duv; + # print(f1[1]) + # print(f2[1]) + # print(f3[1]) + + #koefficient a + liambda = self.N_ptcl/(2*sqrt(pi)*self.sigma_s0); + beta3 = beta**3; + gamma4 = gamma**4; + koef_a =c*self.r_i**2*liambda/(16*pi*sqrt(pi)*self.dp_p*beta3*gamma4)/(self.emit_x*self.emit_y); + + rsi=[None]*n_element + rxi=[None]*n_element + ryi=[None]*n_element + #Rates by element + a = koef_a + rx = 0; + ry = 0; + rs = 0; + n=1;#n=2 for coasting + inv_circ = 1/circ; + + for ie in range(n_element-1): + l_element = l_element_v[ie+1]; + rsi[ie] = n*a*(1-os_all[ie].d2)*f1[ie]*l_element*inv_circ; + rxi[ie] = a*(f2[ie]+f1[ie]*(os_all[ie].d2+os_all[ie].dtld*os_all[ie].dtld))*l_element*inv_circ; + ryi[ie] = a*f3[ie]*l_element*inv_circ; + rx += rxi[ie]; + ry += ryi[ie]; + rs += rsi[ie]; + + + rates=[rx,ry,rs] + rates_fix=[] + + for r in rates: + if np.abs([r])<1e-8: + rates_fix.append(0) + else: + rates_fix.append(r) + + self.rsx=rates_fix[0] + self.rsy=rates_fix[1] + self.rss=rates_fix[2] + +class IBS_BM_rates(object): + + def __init__(self, twiss,beam_dict,bunch,nu_x,nu_y,nu_z,log_c,fact): + + self.beam_dict=beam_dict + self.twiss=twiss + + self.nu_x=nu_x + self.nv_y=nu_y + self.nu_z=nu_z + + self.log_c=log_c + + + #Ion beam parameters + + self.Z = self.beam_dict['Z'] + self.mass =self.beam_dict['mass'] + self.KE=self.beam_dict['KE'] + + self.k_ke = 8.9875517873681764E9; #Coulomb constant + self.emit_x=self.beam_dict['emit_x'] + self.emit_y=self.beam_dict['emit_y'] + self.dp_p=self.beam_dict['dp_p'] + self.N_ptcl=self.beam_dict['N_ptcl'] + self.sigma_s0=self.beam_dict['sigma_s0'] + + gamma = 1+self.KE/self.mass; + beta = sqrt(gamma*gamma-1)/gamma; + dp_p2= self.dp_p**2 + inv_dp_p2=1/dp_p2 + gamma2=gamma**2 + + self.gamma=gamma + self.beta=beta + + self.r_i = self.k_ke*self.Z**2*e*1e-6/self.mass; #Classical ion radius + self.bunched = bunch; + + self.energy_spread = beta*beta*self.dp_p; + p0= gamma*self.mass*1e6*e*beta/c; + + ##Lattice functions + s_v=self.twiss['s'] + betx_v=self.twiss['betx'] + alfx_v=self.twiss['alfx'] + mux_v=self.twiss['mux'] + dx_v=self.twiss['dxst']*fact + dpx_v=self.twiss['dpxst']*fact + bety_v=self.twiss['bety'] + alfy_v=self.twiss['alfy'] + muy_v=self.twiss['muy'] + dy_v=self.twiss['dyst']*fact + dpy_v=self.twiss['dpyst']*fact + + l_element_v=[0] + + n_element=len(betx_v) + + for n in range(n_element-1): + l_element_v.append(s_v[n+1]-s_v[n]) + + circ=twiss.summary.length + + f0 = beta*c/circ; #revolution frequency. + w0= 2*pi*f0; #angular frequency. + + + + def rd( x, y, z): + lolim = 6.E-51 + uplim = 1.0E+48 + + if ( \ + x < 0.0 or \ + y < 0.0 or \ + x + y < lolim or \ + z < lolim or \ + uplim < x or \ + uplim < y or \ + uplim < z ): + print ( '' ) + print ( 'RD - Error!' ) + print ( ' Invalid input arguments.' ) + print ( ' X = %g' % ( x ) ) + print ( ' Y = %g' % ( y ) ) + print ( ' Z = %g' % ( z ) ) + print ( '' ) + ierr = 1 + value = 0.0 + return value + + ierr = 0 + xn = x + yn = y + zn = z + sigma = 0.0 + power4 = 1.0 + + while ( True ): + + mu = ( xn + yn + 3.0 * zn ) * 0.2 + xndev = ( mu - xn ) / mu + yndev = ( mu - yn ) / mu + zndev = ( mu - zn ) / mu + epslon = max ( abs ( xndev ), max ( abs ( yndev ), abs ( zndev ) ) ) + + if ( epslon < 0.01 ): + c1 = 3.0 / 14.0 + c2 = 1.0 / 6.0 + c3 = 9.0 / 22.0 + c4 = 3.0 / 26.0 + ea = xndev * yndev + eb = zndev * zndev + ec = ea - eb + ed = ea - 6.0 * eb + ef = ed + ec + ec + s1 = ed * ( - c1 + 0.25 * c3 * ed - 1.5 * c4 * zndev * ef ) + s2 = zndev * ( c2 * ef + zndev * ( - c3 * ec + zndev * c4 * ea ) ) + value = 3.0 * sigma + power4 * ( 1.0 + s1 + s2 ) / ( mu * np.sqrt ( mu ) ) + return value + + xnroot = np.sqrt ( xn ) + ynroot = np.sqrt ( yn ) + znroot = np.sqrt ( zn ) + lamda = xnroot * ( ynroot + znroot ) + ynroot * znroot + sigma = sigma + power4 / ( znroot * ( zn + lamda ) ) + power4 = power4 * 0.25 + xn = ( xn + lamda ) * 0.25 + yn = ( yn + lamda ) * 0.25 + zn = ( zn + lamda ) * 0.25 + + #BM optical storage + class Kernels : + psi=sx=sp=sxp=inv_sigma=0 + + def __init__(self,psi,sx,sp,sxp,inv_sigma): + self.psi=psi; + self.sx=sx; + self.sp=sp; + self.sxp=sxp; + self.inv_sigma=inv_sigma; + class OpticalStorageBM : #variables only depends on the TWISS parameters and the energy. + phi=dx2=dx_betax_phi_2=sqrt_betay=gamma_phi_2=0 + + def __init__(self,phi,dx2,dx_betax_phi_2,sqrt_betay,gamma_phi_2): + self.phi=phi; + self.dx2=dx2; #D_x * D_x + self.dx_betax_phi_2=dx_betax_phi_2; # D_x * D_x / (beta_x * beta_x) + phi * phi + self.sqrt_betay=sqrt_betay; # sqrt(beta_y) + self.gamma_phi_2=gamma_phi_2; # gamma * gamma * phi * phi + + ###Solver section### + + sigma_xbet=[None] * n_element + sigma_y=[None] * n_element + alfx2=[None] * n_element + sigma_xp=[None] * n_element + alfy2=[None] * n_element + sigma_yp=[None] * n_element + + #Bunch properties at each element + for i in range(n_element): + sigma_xbet[i]=sqrt(betx_v[i]*self.emit_x); + sigma_y[i]=sqrt(bety_v[i]*self.emit_y); + alfx2[i]=alfx_v[i]**2; + sigma_xp[i]=sqrt((1+alfx2[i])*self.emit_x/betx_v[i]); + alfy2[i]=alfy_v[i]**2; + sigma_yp[i]=sqrt((1+alfy2[i])*self.emit_y/bety_v[i]); + + optical_strage=[None] * n_element + for i in range(n_element): + + dx_betax = dx_v[i]/betx_v[i]; + dx2 = dx_v[i] * dx_v[i]; + phi = dpx_v[i]+alfx_v[i]*dx_betax; + dx_betax_phi_2 = dx_betax*dx_betax + phi*phi; + sqrt_betay = sqrt(bety_v[i]); + gamma_phi_2 = gamma*gamma*phi*phi; + optical_strage[i]= OpticalStorageBM(phi,dx2,dx_betax_phi_2,sqrt_betay,gamma_phi_2); + + kernels=[None] * n_element + + for i in range(n_element): + betx = betx_v[i]; + bety = bety_v[i]; + dx = dx_v[i]; + dpx = dpx_v[i]; + sigma_x = sqrt(sigma_xbet[i]*sigma_xbet[i]+dx*dx*self.dp_p*self.dp_p); + sigma_y = sqrt(bety_v[i]*self.emit_y); + inv_sigma = 1/(sigma_x*sigma_y); + + ax = betx/self.emit_x; + lamda_1 = bety/self.emit_y; #lamda_1 = ay. + asa = ax*optical_strage[i].dx_betax_phi_2 + inv_dp_p2; + a1 = gamma2*asa; + a2 = (ax-a1)/2; + a1 = (ax+a1)/2; + + lamda_2 = sqrt(a2*a2+ax*ax*optical_strage[i].gamma_phi_2); + lamda_3 = a1 - lamda_2; + tmp1 = 3/lamda_2; + lamda_2 = a1 + lamda_2; + + inv_lamda_1 = 1/lamda_1; + inv_lamda_2 = 1/lamda_2; + inv_lamda_3 = 1/lamda_3; + + r1 = rd(inv_lamda_2, inv_lamda_3, inv_lamda_1); + r2 = rd(inv_lamda_3, inv_lamda_1, inv_lamda_2); + r3 = 3*sqrt(lamda_1*lamda_2*lamda_3)-r1-r2; + + r1 *= inv_lamda_1*2; + r2 *= inv_lamda_2; + r3 *= inv_lamda_3; + + psi = -r1 + r2 + r3; + + sxp = tmp1*ax*optical_strage[i].gamma_phi_2*(r3-r2); + tmp1 = tmp1*a2; + tmp2 = 1 + tmp1; + tmp1 = 1 - tmp1; + sp = gamma2*(r1-r2*tmp1-r3*tmp2)/2; + sx = (r1-r2*tmp2-r3*tmp1)/2; + kernels[i] = Kernels(psi,sx,sp,sxp,inv_sigma); + + lambd = 1/(2*sqrt(pi)*self.sigma_s0); + + beta3 = beta*beta*beta; + gamma5 = gamma*gamma*gamma*gamma*gamma; + + lc=log_c + c_bm=lambd*self.N_ptcl*self.r_i**2*c/(circ*6*sqrt(pi)*beta3*gamma5); + c_bm=c_bm*lc; + + rsi=[None]*n_element + rxi=[None]*n_element + ryi=[None]*n_element + + rx = 0; + ry = 0; + rs = 0; + n=1; + + for ie in range(n_element-1): + l_element = l_element_v[ie+1]; + rxi[ie] = (betx_v[ie]*kernels[ie].inv_sigma*(kernels[ie].sx+optical_strage[ie].dx_betax_phi_2*kernels[ie].sp+kernels[ie].sxp)*l_element)*c_bm/self.emit_x; + ryi[ie] = (bety_v[ie]*kernels[ie].inv_sigma*kernels[ie].psi*l_element)*c_bm/self.emit_y; + rsi[ie] = (kernels[ie].inv_sigma*kernels[ie].sp*l_element)*n*c_bm/(dp_p2); + rx += rxi[ie]; + ry += ryi[ie]; + rs += rsi[ie]; + + rates=[rx,ry,rs] + rates_fix=[] + + for r in rates: + if np.abs([r])<1e-8: + rates_fix.append(0) + else: + rates_fix.append(r) + + self.rsx=rates_fix[0] + self.rsy=rates_fix[1] + self.rss=rates_fix[2] + +class IBS_kick(object): + + def __init__(self, growthrate_x, growthrate_y, growthrate_z, dt,Dx=0,Dy=0): + self.rxs=growthrate_x + self.rys=growthrate_y + self.rss=growthrate_z + self.dt=dt + self.Dx=Dx + self.Dy=Dy + + def track(self, beam): + + betagamma=beam.betagamma + + emit_x=beam.epsn_x()/betagamma + emit_y=beam.epsn_y()/betagamma + + num_sample=beam.macroparticlenumber + + betaX=1.971964#beam.beta_Twiss_x() + betaY=3.41751#2#beam.beta_Twiss_y() + SigDp=beam.sigma_dp() + beta_z=beam.sigma_z()/beam.sigma_dp() + + + if self.rxs>0: + beam.xp +=np.random.normal(0,np.sqrt(2*self.rxs*self.dt*emit_x/betaX),num_sample) + else: + beam.xp *=np.exp(self.rxs*self.dt) + + if self.rys>0: + beam.yp +=np.random.normal(0,np.sqrt(2*self.rys*self.dt*emit_y/betaY),num_sample) + else: + beam.yp *=np.exp(self.rys*self.dt) + + # subtract dispersion + beam.x -= self.Dx* beam.dp + beam.y -= self.Dy* beam.dp + #kick + if self.rss>0: + dp_delta=np.random.normal(0,np.sqrt(2*self.rss*self.dt*SigDp**2),num_sample) + beam.dp +=dp_delta + else: + beam.dp *=np.exp(self.rss*self.dt) + + + # re-add dispersion + beam.x += self.Dx* beam.dp + beam.y += self.Dy* beam.dp + \ No newline at end of file diff --git a/simulation_notebooks/IBS_mod/Quadrupole_settings b/simulation_notebooks/IBS_mod/Quadrupole_settings new file mode 100644 index 0000000..91832f7 --- /dev/null +++ b/simulation_notebooks/IBS_mod/Quadrupole_settings @@ -0,0 +1,373 @@ +!!! Three quadrupole families used in ELENA, each includes 4 equal members +!!! To complete optics definition, one more parameter has to be defined, it is edge angle of bending magnet + +!!! Quad setting for Qx=2.3,Qy=1.3, gap=90mm in BMs, E1=E2=Pi/10, FINT=0.5 +!KQ1:= 2.689573E+00; +!KQ2:=-1.658787E+00; +!KQ3:= 8.332620E-01; +!KSF= 22.19; +!KSD= -42.095; + +!!! Quad setting for Qx=2.3,Qy=1.3, gap=98mm in BMs, E1=E2=Pi/10, FINT=0.5 +!KQ1:= 2.70436E+00; +!KQ2:=-1.76217E+00; +!KQ3:= 8.43552E-01; + +!KSF= 20.9281; +!KSD=-39.9785*0; + +!!! Qx=2.3,Qy=1.3, gap=76mm, E1=E2=Pi*17/180, FINT=0.5 +!KQ1:= 2.35479E+00; +!KQ2:=-1.38968E+00; +!KQ3:= 7.14651E-01; + +!!! Qx=2.3,Qy=1.3, gap=76mm, E1=E2=Pi*17/180, FINT=0.424 +!KQ1:= 2.26868e+00; +!KQ2:=-1.19607e+00; +!KQ3:= 7.18355e-01; + +!!! For E1=E2=15 degrees +!KQ1:= 2.10E+00; +!KQ2:=-1.38E+00; +!KQ3:= 3.88E-01; + +!!! For E1=E2=19 degrees +!KQ1:= 2.60E+00; +!KQ2:=-1.44E+00; +!KQ3:= 8.37E-01; + +!!! For Qx=1.3,Qy=2.3, E1=E2=29 degrees +!KQ1:= 2.96953; +!KQ2:=-2.55010; +!KQ3:= 1.22255; +!KSF= 3.70028; +!KSD=-7.09711; + +!!! Qx=1.8,Qy=1.8, gap=76mm, E1=E2=17 degrees, FINT=0.5 +!KQ1:= 0.98 !2.23; +!KQ2:=-1.24; +!KQ3:= 0.35; + +/*********************************************************************************************** + ---------------------- Optics settings choice during design studies ------------------------- +***********************************************************************************************/ +! Qx=2.3,Qy=1.3, BM gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=Pi/3*0.927m ! Initial basic choice +!KQ1:= 2.27646e+00; +!KQ2:=-1.20793e+00; +!KQ3:= 7.19841e-01; +!KSF:=23.7315; +!KSD:=-45.4831; +!!! cooler: betx=2.32m, bety=2.39m, Dx=1.50m, betx(max)=11.8m, bety(max)=6.6m + +! Qx=2.309,Qy=1.309, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=Pi/3*0.927m +!KQ1:= 2.35527e+00; +!KQ2:=-1.30180e+00; +!KQ3:= 7.19724e-01; +!KSF:=23.7315; +!KSD:=-45.4831; +!!! cooler: betx=2.08m, bety=2.96m, Dx=1.43m, betx(max)=12.0m, bety(max)=6.4m + +! Qx=2.28,Qy=1.30, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! cooler: betx=2.47m, bety=2.48m, Dx=1.17m, betx(max)=14.1m, bety(max)=6.6m +!KQ1:=2.37320e+00; +!KQ2:=-1.18231e+00; +!KQ3:=6.48274e-01; + +! Qx=2.28,Qy=1.28, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! cooler: betx=2.36m, bety=2.21m, Dx=1.48m, betx(max)=12.5m, bety(max)=7.3m +!KQ1:=2.19284e+00; +!KQ2:=-1.02431e+00; +!KQ3:=6.87485e-01; +!KSF:=25.8331; +!KSD:=-49.6903; + +! Qx=2.27,Qy=1.27, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! cooler: betx=2.48m, bety=2.20m, Dx=1.49m, betx(max)=12.9m, bety(max)=7.6m +!KQ1:=2.13978e+00; +!KQ2:=-9.26563e-01; +!KQ3:=6.75039e-01; +!KSF:=26.5249; +!KSD:=-50.9945; + +! Qx=2.27,Qy=1.30, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! cooler: betx=2.67m, bety=2.64m, Dx=0.96m, betx(max)=15.5m, bety(max)=6.6m +!KQ1:=2.43083e+00; +!KQ2:=-1.18109e+00; +!KQ3:=6.16973e-01; +!!! Sol on, 100G, e-beam on: Qx=2.266, Qy=1.296 (dp=0), +!!! Qx=2.241, Qy=1.294 (dp=0.5%) +!!! Qx=2.292, Qy=1.299 (dp=-0.5%) +!!! Betx(max)=16.0m, Bety(max)=6.4m (dp=0) +!!! Betx,y(EC) = 2.60/2.55, Dx(EC) = 0.97 + +! Qx=2.26,Qy=1.29, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! cooler: betx=2.65m, bety=2.50m, Dx=1.05m, betx(max)=15.8m, bety(max)=6.6m +!KQ1:=2.36677e+00; +!KQ2:=-1.10585e+00; +!KQ3:=6.49360ee-01; + +! Qx=2.26,Qy=1.27, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! cooler: betx=2.65m, bety=2.29m, Dx=1.38m, betx(max)=13.8m, bety(max)=7.6m +!KQ1:=2.15794e+00; +!KQ2:=-8.99080e-01; +!KQ3:=6.49360e-01; +!KSF:=2.87864e+01; +!KSD:=-5.23605e+01; + +!!! Qx=2.46,Qy=1.46, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! cooler: betx=1.20 m, bety=2.16m +!KQ1:= 2.90619e+00; +!KQ2:= -2.43108e+00; +!KQ3:= 8.17551e-01; + +!!! Qx=2.46,Qy=1.46, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! cooler: betx=1.30 m, bety=2.39m, Dx=1.08m, Betmax=13.4/4.5, Dxmax=1.5m +!KQ1:= 3.06772e+00; +!KQ2:= -2.49675e+00; +!KQ3:= 7.39252e-01; + +!!! Qx=2.46,Qy=1.46, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! slighly bigger betx=1.4m, bety=2.4m needs bigger acceptances +!KQ1:= 3.27791e+00; +!KQ2:= -2.54332e+00; +!KQ3:= 7.67173e-01; + +!!! Qx=2.44,Qy=1.46, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! cooler: betx=1.45 m, bety=2.39 Dxm=1.38 m +!KQ1:= 3.19236e+00; +!KQ2:= -2.48155e+00; +!KQ3:= 7.65935e-01; + +!!! Qx=2.44,Qy=1.46, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! cooler: betx=1.45 m, bety=2.39 Dxm=1.38 m +!!!KQ1:= 3.844e+00; +!!!KQ2:= -2.355e+00; +!!!KQ3:= 1.275e+00; + +!!! Qx=2.44,Qy=1.46, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! cooler on: betx=1.36 m, bety=2.74 m, Dxm=1.03 m +!!! betx(max)/bety(max)/Dx(max) + +! MSOL and SCOMP on +!KQ1:= 2.98992e+00; +!KQ2:= -2.38468e+00; +!KQ3:= 7.10937e-01; + +/*********************************************************************************************** + ---------------------- Optics settings choice for commissioning ------------------------- +***********************************************************************************************/ +!!! Qx=2.44,Qy=1.46, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! cooler on: betx=1.38 m, bety=2.43 m, Dxm=1.05 m +!!! betx(max)/bety(max)/Dx(max) 10.1/4.51/1.50 + +! MSOL and SCOMP off +!!! cooler on: betx=1.38 m, bety=2.43 m, Dxm=1.05 m +!!! betx(max)/bety(max)/Dx(max) 10.1/4.51/1.50 +KQ1:= 3.040329655e+00; +KQ2:= -2.459580048e+00; +KQ3:= 0.7135949997e+00; +KSF:=0.; ! 23.9669; +KSD:=0.; !-37.7605; + +!!! for E1=E2=16.45 degrees, Qx=2.34, Qy=1.36 +!KQ1:= 2.601e+00; +!KQ2:= -1.725e+00; +!KQ3:= 6.252e-01; + +!!! Qx=2.43,Qy=1.41, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! cooler: betx=1.60m, bety=2.66m, Dx=0.79m, betx(max)=15.6m, bety(max)=4.6m +!KQ1:= 3.21723e+00; +!KQ2:= -2.30158e+00; +!KQ3:= 7.17592e-01; + +!!! Qx=2.41,Qy=1.43, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! cooler: betx=1.70m, bety=2.69m, Dx=0.71m, betx(max)=16.1m, bety(max)=4.6m +!KQ1:= 3.23872e+00; +!KQ2:= -2.35774e+00; +!KQ3:= 6.88295e-01; +!KSF= 31.6108; +!KSD=-41.640; +!!! Sol on, 100G, e-beam on: Qx=2.409, Qy=1.426 (dp=0), +!!! Qx=2.242, Qy=1.382 (dp=0.5%) +!!! Qx=2.185, Qy=1.185 (dp=-0.5%) +!!! Betx(max)=16.4m, Bety(max)=4.7m (dp=0) +!!! Betx,y(EC) = 1.58/2.43, Dx(EC) = 0.71 + +!!! Qx=2.40,Qy=1.42, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! cooler: betx=1.74m, bety=2.69m, Dx=0.75m, betx(max)=15.7m, bety(max)=4.7m +!KQ1:= 3.17017e+00; +!KQ2:= -2.27526e+00; +!KQ3:= 6.89292e-01; +!KSF= ; +!KSD=; +!!! Sol on, 100G, e-beam on: Qx=2.40, Qy=1 42(dp=0), +!!! Qx=2.242, Qy=1.382 (dp=0.5%) +!!! Qx=2.185, Qy=1.185 (dp=-0.5%) +!!! Betx(max)=m, Bety(max)=m (dp=0) +!!! Betx,y(EC) = /, Dx(EC) = + +!!! Qx=2.38,Qy=1.36, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! cooler: betx=1.70 m, bety=2.66 m, Dx=1.04 m, betx(max)=13.3 m, bety(max)=4.73 m +!KQ1:= 2.722872201e+00; +!KQ2:= -1.860329745e+00; +!KQ3:= 0.6634130943e+00; +!KSF=2.82361e+01 ; +!KSD=-4.30313e+01 ; + +!!! Qx=2.38,Qy=1.35, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! cooler: betx=1.81 m, bety=2.58 m, Dx=0.98 m, betx(max)=14.1 m, bety(max)=5.2 m +!KQ1:= 2.87132e+00; +!KQ2:= -1.84121e+00; +!KQ3:= 7.16335e-01; +!KSF=2.82361e+01 ; +!KSD=-4.30313e+01 ; + +!!! Qx=2.37,Qy=1.34, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! cooler: betx=1.84 m, bety=2.53 m, Dx=1.04 m, betx(max)=13.7 m, bety(max)=5.4 m +!KQ1:= 2.79186e+00; +!KQ2:= -1.748911e+00; +!KQ3:= 7.18483e-01; +!KSF=2.76013e+01 ; +!KSD=-4.32724e+01 ; + +!!! Qx=2.36,Qy=1.33, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! cooler: betx=1.84 m, bety=2.53 m, Dx=1.04 m, betx(max)=13.7 m, bety(max)=5.4 m +!KQ1:= 2.84381e+00; +!KQ2:= -1.70628e+00; +!KQ3:= 6.79473e-01; +!KSF=3.16846e+01 ; +!KSD=-4.42928e+01 ; + +!!! Qx=2.355,Qy=1.32, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! cooler: betx=2.02 m, bety=2.66 m, Dx=-.88 m, betx(max)=14.8 m, bety(max)=5.8 m +!KQ1:= 2.79400e+00; +!KQ2:=-1.63347e+00; +!KQ3:= 6.80319e-01; +!KSF= 3.14335e+01 ; +!KSD= -4.45486e+01 ; + +!!! Qx=2.345,Qy=1.31, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! cooler: betx=2.10 m, bety=2.66 m, Dx=0.87 m, betx(max)=15.0 m, bety(max)=6.0 m +!KQ1:= 2.74028e+00; +!KQ2:=-1.54680e+00; +!KQ3:= 6.73258e-01; +!KSF= 3.18740e+01 ; +!KSD= -4.50223e+01 ; + +!!! Qx=2.37,Qy=1.39, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! betx=1.80m, bety=2.55m, Dx=1.0m, betx(max)=13.8m, bety(max)=4.6m +!KQ1:= 2.89408e+00; +!KQ2:= -1.98930e+00; +!KQ3:= 7.11227e-01; + +!!! Qx=2.23,Qy=1.23, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! cooler: Betx(EC)=3.1/ 2.1, BETmax = 15.0 / 9.3 +!KQ1:= 1.94847e+00; +!KQ2:= -5.34602e-01; +!KQ3:= 6.22351e-01; + +!!! Qx=2.225,Qy=1.205, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! Betx(EC)=3.3/ 2.0, BETmax = 16.0 / 10.4, Dx (EC) = 1.4m +!KQ1:= 1.92190; +!KQ2:= -4.02044e-01; +!KQ3:= 6.09452e-01; +!! Solenoid off, e-beam off +!KSF= 32.3304; +!KSD=-56.241; + +!!! Qx=2.220,Qy=1.205, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! Betx(EC) = 3.4/2.0, BetMax = 16.1/10.5, Dx(EC) = 1.5, DQ1=-3.0,DQ2=-3.1 +!KQ1:= 1.89021; +!KQ2:= -3.70147e-01; +!KQ3:= 6.07140e-01; +!! Solenoid off, e-beam off +!KSF= 32.3304; +!KSD=-56.241; +!!! Sol on, 100G, : Qx=2.232, Qy=1.214 (dp=0), +!!! Sol on, 100G, e-beam on: Qx=2.215, Qy=1.203 (dp=0), +!!! Qx=2.241, Qy=1.209 (dp=-0.5%) +!!! Qx=2.186, Qy=1.195 (dp=0.5%) +!!! Betx(max)=16.7m, Bety(max)=10.2m (dp=0) +!!! Betx,y(EC) = 3.4/2.0, Dx(EC) = 1.5 + +!!! Qx=2.220,Qy=1.195 (bare), gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! Betx(EC) = 3.5/2.0, BetMax = 16.8/11.0, Dx(EC) = 1.3, DQ1=-3.1,DQ2=-3.1 +!KQ1:= 1.92299; +!KQ2:= -3.44745e-01; +!KQ3:= 5.97745e-01; +!! Solenoid off, e-beam off +!KSF= 32.3304; +!KSD=-56.241; +!!! Sol on, 100G, : Qx=2.232, Qy=1.204 (dp=0), +!!! Sol on, 100G, e-beam on: Qx=2.214, Qy=1.193 (dp=0), +!!! Qx=2.242, Qy=1.199 (dp=-0.5%) +!!! Qx=2.185, Qy=1.185 (dp=0.5%) +!!! Betx(max)=17.5m, Bety(max)=11.6m (dp=0) +!!! Betx,y(EC) = 3.5/2.0, Dx(EC) = 1.4 + +!!! Qx=2.205,Qy=1.225, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! Betx(EC)=3.5/ 2.0, BETmax = 15.8 / 9.7 , Dx (EC) = 1.7m +!KQ1:= 1.80789; +!KQ2:= -3.60641e-01; +!KQ3:= 6.04559e-01; + +!KSF= 23.1960; +!KSD=-42.7944; + +!!! Qx=2.21,Qy=1.23, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! Betx(EC)=3.5/ 2.2, BETmax = 16.6 / 9.3 , Dx (EC) = 1.4m +!KQ1:= 1.91979e+00; +!KQ2:= -4.52440e-01; +!KQ3:= 5.93859e-01; + +!KSF= 32.9738; +!KSD=--59.9510; + +!!! Qx=2.205,Qy=1.220, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! Betx(EC)=3.6 / 2.1, Dx (EC) = 1.6m, BETmax = 16.4 / 9.9 +!KQ1:= 1.83666; +!KQ2:= -3.53644e-01; +!KQ3:= 5.97472e-01; +!!! Sol on, 100G, : Qx=2.218, Qy=1.229 (dp=0), +!!! Sol on, 100G, e-beam on: Qx=2.199, Qy=1.218 (dp=0), +!!! Qx=2.226, Qy=1.224 (dp=-0.5%) +!!! Qx=2.170, Qy=1.209 (dp=0.5%) +!!! Betx(max)=17.1m, Bety(max)=9.7m (dp=0) +!!! Betx,y(EC) = 3.6/2.0, Dx(EC) = 1.7 + +!KSF= 31.4018; +!KSD=-60.4056; + +!!! Qx=2.20,Qy=1.22, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! Betx(EC)=3.7/ 2.1, BETmax = 16.9 / 9.9 , Dx (EC) = 1.57m +!KQ1:= 1.83215e+00; +!KQ2:= -3.34424e-01; +!KQ3:= 5.90297e-01; + +!KSF= 32.9738; +!KSD=--59.9510; + +!! Qx=2.20,Qy=1.20, gap=76mm, E1=E2=Pi*17/180, FINT=0.424, Lbm=0.927m +!!! Betx(EC)=3.7/ 2.1, BETmax = 16.9 / 9.9 , Dx (EC) = 1.57m +!KQ1:= 1.80039e+00; +!KQ2:= -2.38011e-01; +!KQ3:= 5.89587e-01; + +!KSF= 32.9738; +!KSD=--59.9510; + + +RETURN; + +KQ1:= 3.04032e+00; ! +KQ2:= -2.45958e+00; ! +KQ3:= 7.13597e-01; + +!KQ1:= 3.04032e+00; +!KQ2:= -KQ1*0.74; +!KQ3:= KQ1*0.25; + + + +RETURN; + diff --git a/simulation_notebooks/IBS_mod/elena_simplified_origin.madx b/simulation_notebooks/IBS_mod/elena_simplified_origin.madx new file mode 100644 index 0000000..f25788e --- /dev/null +++ b/simulation_notebooks/IBS_mod/elena_simplified_origin.madx @@ -0,0 +1,122 @@ +!/***************************************************** +! * +! * MADX file for ELENA low energy antiproton deaccelerator +! * simplified version than elena.madx +! * +! * July 2018 - davide.gamba@cern.ch +! *****************************************************/ + +set, format="-20s"; +set, format="8.5f"; +OPTION, ECHO, WARN=false, -INFO, VERBOSE=false; + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! /* Machine Elements definition */ +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +call, file = "ELENA_elements.def"; +call, file = "ELENA_ring_mod.seq"; + + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! /* Quadrupole strengths */ +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +call, file = "Quadrupole_settings"; + + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! BEAM, USE +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! do computation assuming 100 MeV/c beam momentum (injection) +!BEAM,PARTICLE=ANTIPROTON, PC=0.100; ! 100 MeV/c + +MREST0:=0.9383; +MREST2:=MREST0*MREST0; + +Ek:=0.0001; +E0:=MREST0+Ek; + +PC0:=0.0137; !13.7 MeV/c +gamma:=E0/MREST0; ! Lorentz factor +bc0:=PC0/(MREST0*gamma); ! relativistic beta + +beam_def: Beam, Particle = antiproton, PC=PC0, Npart = 4500000 +, BUNCHED,EX= 0.5e-6, EY= 0.1e-6, SIGE= 2.1312e-08, SIGT= 0.3284; + + +!USE, PERIOD=ELENA; + +seqedit, sequence=ELENA; +flatten; +endedit; + + +use, sequence=ELENA; + +!MATCH, SEQUENCE=ELENA; +!GLOBAL, Q1=2.454;//H-tune +!GLOBAL, Q2=1.416;//V-tune +!!GLOBAL, DQ1=1.416;//V-tune +!VARY, NAME= KQ1, STEP=0.000001; +!VARY, NAME = KQ2, STEP=0.000001; +!LMDIF, CALLS=100, TOLERANCE=1e-6;//method adopted +!ENDMATCH; + +!/******************************************************************************* +! ******************************************************************************* +! * ZERO CHROMATICITY. MATCHING SEXTUPOLES WITH PTC +! ******************************************************************************* +! *******************************************************************************/ +! option, -info; +! option, -echo; +! + + +! ******************************************************************************* +! * MAD-X macros +! ******************************************************************************* +maketwissptc : macro={ + + SELECT,flag=ptc_twiss,clear; + ! TIME=FALSE gives correct chromaticity calculation by default, i.e. no need to multiply dispersion by relativistic beta + PTC_CREATE_UNIVERSE; + PTC_CREATE_LAYOUT, MODEL=2, METHOD=6, NST=5, TIME=false, EXACT; ! TIME=FALSE gives correct chromaticity calculation + SELECT, flag=ptc_twiss,column=name,keyword,l,s,betx,alfx,mu1,disp1,disp2,bety,alfy,mu2,disp3,disp4; + PTC_TWISS, CLOSED_ORBIT, TABLE=PTC_TWISS, ICASE=5, NO=3, SUMMARY_TABLE; !, SUMMARY_FILE="summary", FILE="ptc_twiss.out" ; + PTC_END; +} +maketwiss : macro={ + + ! note that chromatic quantities (e.g. Dispersion) need to be multiplied by relativistic beta to obtain actual values. + dxst :=table(twiss,dx)*bc0; ! standard dispersion calculated from the MADX dispersion for low energies + dyst :=table(twiss,dy)*bc0; + dpyst :=table(twiss,dpy)*bc0; + dpxst :=table(twiss,dpx)*bc0; + sigma_x :=sqrt(table(twiss,betx)*1E-6)*1000; !mm + sigma_y :=sqrt(table(twiss,bety)*75E-8)*1000; !mm + select, flag=twiss, clear; + !select, flag=twiss,column=name,keyword,s,betx,alfx,mux,dxst,dpxst,bety,alfy,muy,dyst,dpyst,sigma_x,sigma_y,k1l,k3l; + select, flag=twiss,column=name,keyword,s,betx,alfx,mux,dxst,dpxst,bety,alfy,muy,dyst,dpyst; + twiss; +} + + +! ******************************************************************************* +! * do something +! ******************************************************************************* +exec, maketwissptc; +plot, noversion=true, table=ptc_twiss, ptc=true, haxis=s, vaxis1=betx,bety, vaxis2=disp1, interpolate, + hmin=0, hmax=30.4,title="PTC Twiss", colour=100, file=simpleTest; +write, table=ptc_twiss, file=NoEcoolTwissPTC.tfs; + + +exec, maketwiss; +write, table=twiss, file=NoEcoolTestTwiss.tfs; + +plot, noversion=true, table=twiss, haxis=s, vaxis1=betx,bety, vaxis2=dx, interpolate, + hmin=0, hmax=30.4,title="Normal Twiss", colour=100, file=simpleTest; +plot, noversion=true, table=twiss, haxis=s, vaxis1=sigma_x,sigma_y, vaxis2=dx, interpolate, + hmin=0, hmax=30.4,title="Normal Twiss", colour=100, file=simpleTest; + +IBS, FILE=ibs_elena_madx;