DetunersTest.ipynb Open in SWAN Download

general imports

In [1]:
from __future__ import division, print_function

import numpy as np
np.random.seed(42)

from scipy.constants import m_p, c, e

import matplotlib.pyplot as plt
%matplotlib inline
In [2]:
# sets the PyHEADTAIL directory etc.
try:
    from settings import *
except:
    pass

PySussix imports

Sussix is a tune analysis tool: https://cds.cern.ch/record/702438/files/sl-note-98-017.pdf . If you do not have PySussix, its python interface, find it under https://github.com/PyCOMPLETE/PySussix .

In [3]:
try:
    from PySussix import Sussix
except ImportError:
    print ('ERROR: This interactive test needs the PySussix module. Trying PySUSSIX')
    from PySUSSIX import Sussix
    print('PySUSSIX found')

PyHEADTAIL imports

In [4]:
from PyHEADTAIL.trackers.transverse_tracking import TransverseMap
from PyHEADTAIL.trackers.detuners import Chromaticity, AmplitudeDetuning
import PyHEADTAIL.particles.generators as generators
PyHEADTAIL v1.10.5.269


Setting up the machine and functions

In [5]:
# Basic parameters.
n_turns = 1024
n_segments = 1
n_macroparticles = 500

Q_x = 64.28
Q_y = 59.31
Q_s = 0.0020443

C = 26658.883
R = C / (2.*np.pi)

alpha_x_inj = 0.
alpha_y_inj = 0.
beta_x_inj = 66.0064
beta_y_inj = 71.5376
alpha_0 = 0.0003225

# Parameters for transverse map.
s = np.arange(0, n_segments + 1) * C / n_segments

alpha_x = alpha_x_inj * np.ones(n_segments)
beta_x = beta_x_inj * np.ones(n_segments)
D_x = np.zeros(n_segments)

alpha_y = alpha_y_inj * np.ones(n_segments)
beta_y = beta_y_inj * np.ones(n_segments)
D_y = np.zeros(n_segments)
In [6]:
def calc_sussix_spec(x, xp, y, yp, p_idx, turn, window_width, q_x, q_y, n_lines=10):
    # Initialise Sussix object
    SX = Sussix()
    SX.sussix_inp(nt1=1, nt2=window_width, idam=2, ir=0, tunex=q_x, tuney=q_y)

    tunes_x = np.empty(n_lines)
    tunes_y = np.empty(n_lines)

    SX.sussix(x[p_idx,turn:turn+window_width], xp[p_idx,turn:turn+window_width],
              y[p_idx,turn:turn+window_width], yp[p_idx,turn:turn+window_width],
              x[p_idx,turn:turn+window_width], xp[p_idx,turn:turn+window_width]) # this line is not used by sussix!

    return SX.ox[:n_lines], SX.oy[:n_lines]

def track_n_save(bunch, map_):
    n_particles = bunch.macroparticlenumber
    x_i = np.empty((n_particles, n_turns))
    xp_i = np.empty((n_particles, n_turns))
    y_i = np.empty((n_particles, n_turns))
    yp_i = np.empty((n_particles, n_turns))

    for i in xrange(n_turns):
        x_i[:,i] = bunch.x[:]
        xp_i[:,i] = bunch.xp[:]
        y_i[:,i] = bunch.y[:]
        yp_i[:,i] = bunch.yp[:]

        for m_ in map_:
            m_.track(bunch)

    return x_i, xp_i, y_i, yp_i

def analyse_n_plot(bunch, x_i, xp_i, y_i, yp_i):
    ox = np.empty(bunch.macroparticlenumber)
    oy = np.empty(bunch.macroparticlenumber)
    print ('analysing particle spectra ... this may take some time.')
    for p_idx in range(bunch.macroparticlenumber):
        ox[p_idx], oy[p_idx] = calc_sussix_spec(x_i, xp_i, y_i, yp_i, p_idx,
                                   turn=0, window_width=512, q_x=Q_x%1, q_y=Q_y%1, n_lines=1)
        if (p_idx)%100 == 0:
            print ('particle', p_idx)

    fig = plt.figure(figsize=(20,20))
    ax1 = fig.add_subplot(311)
    ax2 = fig.add_subplot(312, sharex=ax1)
    ax3 = fig.add_subplot(313)

    ax1.scatter(ox, oy)
    ax1.set_ylabel(r'$Q_y$', fontsize=20)
    ax1.set_xlabel(r'$Q_x$', fontsize=20)
    ax2.hist(ox, bins=50, color='blue')
    ax2.set_xlabel(r'$Q_x$', fontsize=20)
    ax3.hist(oy, bins=50, color='red')
    ax3.set_xlabel(r'$Q_y$', fontsize=20)

    print ('std dev. Qx', np.std(ox))
    print ('std dev. Qy', np.std(oy))
    plt.show()

    os.remove('sussix.inp')

def generate_bunch(n_macroparticles, alpha_x, alpha_y, beta_x, beta_y, alpha_0, Q_s, R):
    intensity = 1.05e11
    sigma_z = 0.059958
    gamma = 3730.26
    eta = alpha_0 - 1. / gamma**2
    gamma_t = 1. / np.sqrt(alpha_0)
    p0 = np.sqrt(gamma**2 - 1) * m_p * c

    beta_z = eta * R / Q_s

    epsn_x = 3.75e-6 # [m rad]
    epsn_y = 3.75e-6 # [m rad]
    epsn_z = 4 * np.pi * sigma_z**2 * p0 / (beta_z * e)

    bunch = generators.generate_Gaussian6DTwiss(
        macroparticlenumber=n_macroparticles, intensity=intensity, charge=e,
        gamma=gamma, mass=m_p, circumference=C,
        alpha_x=alpha_x, beta_x=beta_x, epsn_x=epsn_x,
        alpha_y=alpha_y, beta_y=beta_y, epsn_y=epsn_y,
        beta_z=beta_z, epsn_z=epsn_z)

    return bunch

Let's go

Amplitude detuning

Amplitude detuning (such as from octupole magnets) can be modelled with the AmplitudeDetuning detuner class. Check out this LHC example with currents of $I=400~A$ where we expect $\approx 5\times 10^{-4}$ detuning in both transverse planes:

In [7]:
bunch = generate_bunch(
    n_macroparticles, alpha_x_inj, alpha_y_inj, beta_x_inj, beta_y_inj,
    alpha_0, Q_s, R)

ampl_det = AmplitudeDetuning.from_octupole_currents_LHC(i_focusing=400, i_defocusing=-400)
trans_map = TransverseMap(
    s, alpha_x, beta_x, D_x, alpha_y, beta_y, D_y, Q_x, Q_y, [ampl_det])

map_ = [m for m in trans_map]

x_i, xp_i, y_i, yp_i = track_n_save(bunch, map_)
analyse_n_plot(bunch, x_i, xp_i, y_i, yp_i)
analysing particle spectra ... this may take some time.
particle 0
particle 100
particle 200
particle 300
particle 400
std dev. Qx 0.000493943454511
std dev. Qy 0.000504673628496
No description has been provided for this image

Chromatic detuning

Chromaticity can be modelled with the Chromaticity detuner class. First-order chromaticity for values of $Q'_x=6$ and $Q'_y=3$ with LHC tunes yields:

In [8]:
bunch = generate_bunch(
    n_macroparticles, alpha_x_inj, alpha_y_inj, beta_x_inj, beta_y_inj,
    alpha_0, Q_s, R)

chroma = Chromaticity(Qp_x=[6], Qp_y=[3]) # note the lists!
trans_map = TransverseMap(
    s, alpha_x, beta_x, D_x, alpha_y, beta_y, D_y, Q_x, Q_y, [chroma])

map_ = [m for m in trans_map]

x_i, xp_i, y_i, yp_i = track_n_save(bunch, map_)
analyse_n_plot(bunch, x_i, xp_i, y_i, yp_i)
analysing particle spectra ... this may take some time.
particle 0
particle 100
particle 200
particle 300
particle 400
std dev. Qx 0.000561082708596
std dev. Qy 0.00028045654874
No description has been provided for this image

Higher-order Chromatic Detuning

The Chromaticity detuner class supports higher orders in chromatic detuning. Check out this LHC example with horizontal chromaticities of $Q'_x=6, Q''_x = 4\times10^4$ and vertical chromaticities of $Q'_y=3, Q''_y=0, Q'''_y=2\times10^8$:

In [9]:
bunch = generate_bunch(
    n_macroparticles, alpha_x_inj, alpha_y_inj, beta_x_inj, beta_y_inj,
    alpha_0, Q_s, R)

chroma = Chromaticity(Qp_x=[6., 4e4], Qp_y=[3., 0., 2e8]) 
# here the lists are extended beyond first order! ^^^
trans_map = TransverseMap(
    s, alpha_x, beta_x, D_x, alpha_y, beta_y, D_y, Q_x, Q_y, [chroma])

map_ = [m for m in trans_map]

x_i, xp_i, y_i, yp_i = track_n_save(bunch, map_)
analyse_n_plot(bunch, x_i, xp_i, y_i, yp_i)
analysing particle spectra ... this may take some time.
particle 0
particle 100
particle 200
particle 300
particle 400
std dev. Qx 0.000568932172386
std dev. Qy 0.000322448107092
No description has been provided for this image

The corresponding detuning curves versus $\delta = \frac{p-p_0}{p_0}$ show a parabola in the horizontal and a cubic polynomial in the vertical:

In [10]:
fig = plt.figure(figsize=(12,8))

chroma = Chromaticity(Qp_x=[6., 4e4], Qp_y=[3., 0., 2e8])
chroma.generate_segment_detuner(dmu_x=1, dmu_y=1)

dp = np.linspace(-5e-4, 5e-4, 500)

dQx = chroma.segment_detuners[0].calc_detuning_x(dp)
dQy = chroma.segment_detuners[0].calc_detuning_y(dp)

plt.plot(dp, dQx, label='Detuning x')
plt.plot(dp, dQy, label='Detuning y')
plt.grid('on')
plt.xlabel('dp')
plt.ylabel('Detuning')
plt.legend(loc='upper left')
plt.show()
No description has been provided for this image