Twiss parameters output

A. Petrenko (Novosibirsk, 2019)

In [1]:
import numpy as np
import pandas as pd
import holoviews as hv
import os

hv.extension("bokeh")
In [2]:
%opts Curve Spread [width=700 height=300 show_grid=True \
            default_tools=['box_zoom','pan','wheel_zoom','reset']]

Elegant lattice file:

In [3]:
with open("machine.lte", "r") as f:
    print(f.read())

w0: watch,filename="results/w0.sdds",mode=coord
w1: watch,filename="results/w1.sdds",mode=coord
w2: watch,filename="results/w2.sdds",mode=coord
w3: watch,filename="results/w3.sdds",mode=coord

MA20:   MAXAMP, X_MAX=10.0e-3, Y_MAX=10.0e-3, ELLIPTICAL=1
MA16:   MAXAMP, X_MAX=8.00e-3, Y_MAX=8.00e-3, ELLIPTICAL=1

D20:   DRIFT, L=0.10
QF20:  QUAD,  L=0.10, K1 = 7.0
DQF20: DRIFT, L=0.10
RF10:  RFCA, L=3.0, VOLT=-60.0e6, PHASE=90, FREQ=2856e6, change_p0=1

D30:   DRIFT, L=0.10
QD30:  QUAD,  L=0.10, K1 =-3.5
DQD30: DRIFT, L=0.10

D40:   DRIFT, L=0.10
QF40:  QUAD,  L=0.10, K1 = 3.5
DQF40: DRIFT, L=0.10

D50:   DRIFT, L=0.10
QD50:  QUAD,  L=0.10, K1 =-3.5
DQD50: DRIFT, L=0.10


FODO1: LINE=(MA16,&
    D20,QF20,DQF20,MA20,RF10,MA16,&
    D30,QD30,DQD30,MA20,RF10)

FODO: LINE=(MA16,&
    D40,QF40,DQF40,MA20,RF10,MA16,&
    D50,QD50,DQD50,MA20,RF10)

BLINE: LINE=(w0,&
    FODO1,&
    20*FODO,&
w3)

machine: LINE=(w0,-BLINE,w3)
In [4]:
with open("twiss.ele", "r") as f:
    print(f.read())

&transmute_elements name=*, type=WATCH, new_type=DRIF &end
!&transmute_elements name=*, type=RFCA, new_type=DRIF &end

&divide_elements
    name = *
    maximum_length = 0.05
&end

&run_setup
    lattice = machine.lte
    magnets = results/beamline.mag
    centroid = results/beam.cen
    parameters = results/parameters.sdds
    p_central_mev = 2528.0
    use_beamline = machine
    default_order=2
&end

&twiss_output
    filename = results/twiss.twi
    matched = 0
    beta_x = 10.5
    alpha_x = +1.8
    beta_y = 4.0
    alpha_y = -0.75
    eta_x = 0
    eta_y = 0
    etap_x = 0
    etap_y = 0
&end

&run_control &end

&bunched_beam &end
&track &end

Run Elegant

In [5]:
!elegant twiss.ele

Running elegant at Thu Apr 25 15:47:29 2019

This is elegant 35.1.0, Dec 17 2018, by M. Borland, J. Calvey, M. Carla', N. Carmignani, M. Ehrlichman, L. Emery, W. Guo, R. Lindberg, V. Sajaev, R. Soliday, Y.-P. Sun, C.-X. Wang, Y. Wang, Y. Wu, and A. Xiao.
=====================================================================================
Thanks for using elegant.  Please cite the following reference in your publications:
  M. Borland, "elegant: A Flexible SDDS-Compliant Code for Accelerator Simulation,"
  Advanced Photon Source LS-287, September 2000.
If you use a modified version, please indicate this in all publications.
=====================================================================================
Link date: Dec 17 2018 16:14:43, SVN revision: 25872M
statistics:    ET:     00:00:00 CP:    0.00 BIO:0 DIO:0 PF:0 MEM:0
&transmute_elements
    name = *,
    type = WATCH,
    exclude = {NULL},
    new_type = DRIF,
    disable = 0,
    clear_all = 0,
&end
statistics:    ET:     00:00:00 CP:    0.00 BIO:0 DIO:0 PF:0 MEM:0
&divide_elements
    name = *,
    type = {NULL},
    exclude = {NULL},
    divisions = 0,
    maximum_length = 5.000000000000000e-02,
    clear = 0,
&end
statistics:    ET:     00:00:00 CP:    0.00 BIO:0 DIO:0 PF:0 MEM:0
&run_setup
    lattice = machine.lte,
    use_beamline = machine,
    rootname = {NULL},
    output = {NULL},
    centroid = results/beam.cen,
    sigma = {NULL},
    final = {NULL},
    acceptance = {NULL},
    losses = {NULL},
    magnets = results/beamline.mag,
    semaphore_file = {NULL},
    parameters = results/parameters.sdds,
    combine_bunch_statistics = 0,
    wrap_around = 1,
    final_pass = 0,
    default_order = 2,
    concat_order = 0,
    print_statistics = 0,
    show_element_timing = 0,
    monitor_memory_usage = 0,
    random_number_seed = 987654321,
    correction_iterations = 1,
    echo_lattice = 0,
    p_central = 0.000000000000000e+00,
    p_central_mev = 2.528000000000000e+03,
    always_change_p0 = 0,
    load_balancing_on = 0,
    random_sequence_No = 1,
    expand_for = {NULL},
    tracking_updates = 1,
    search_path = {NULL},
    element_divisions = 0,
&end
Seeding random number generators
length of beamline MACHINE per pass: 1.385999999999973e+02 m
statistics:    ET:     00:00:00 CP:    0.03 BIO:0 DIO:0 PF:0 MEM:0
&twiss_output
    filename = results/twiss.twi,
    matched = 0,
    output_at_each_step = 0,
    output_before_tune_correction = 0,
    final_values_only = 0,
    statistics = 0,
    radiation_integrals = 0,
    beta_x = 1.050000000000000e+01,
    alpha_x = 1.800000000000000e+00,
    eta_x = 0.000000000000000e+00,
    etap_x = 0.000000000000000e+00,
    beta_y = 4.000000000000000e+00,
    alpha_y = -7.500000000000000e-01,
    eta_y = 0.000000000000000e+00,
    etap_y = 0.000000000000000e+00,
    reference_file = {NULL},
    reference_element = {NULL},
    reference_element_occurrence = 0,
    reflect_reference_values = 0,
    concat_order = 3,
    higher_order_chromaticity = 0,
    higher_order_chromaticity_points = 5,
    higher_order_chromaticity_range = 4.000000000000000e-04,
    quick_higher_order_chromaticity = 0,
    chromatic_tune_spread_half_range = 0.000000000000000e+00,
    cavities_are_drifts_if_matched = 1,
    compute_driving_terms = 0,
    leading_order_driving_terms_only = 0,
    s_dependent_driving_terms_file = {NULL},
    local_dispersion = 1,
    n_periods = 1,
&end
final Twiss parameters (chromaticity valid for fully second-order calculation only!):
         beta          alpha           nu           eta          eta'       dnu/d(dp/p)   dbeta/(dp/p)     accept.
          m                          1/2pi           m                         1/2pi            m          mm-mrad
--------------------------------------------------------------------------------------------------------------------
  x:  7.557469e+00 -1.087911e+00  4.120081e+00  0.000000e+00  0.000000e+00 -3.799099e+02 -3.459526e+04  5.920224e+00
  y:  1.516758e+00 -1.412802e+00  4.149525e+00  0.000000e+00  0.000000e+00 -8.186052e+02 -2.595400e+03  4.792310e+00
x acceptance limited to 5.920224e-06 by D40 ending at 7.260000e+01 m
y acceptance limited to 4.792310e-06 by D50 ending at 1.023000e+02 m
statistics:    ET:     00:00:01 CP:    0.86 BIO:0 DIO:0 PF:0 MEM:0
&run_control
    n_steps = 1,
    bunch_frequency = 0.000000000000000e+00,
    n_indices = 0,
    n_passes = 1,
    n_passes_fiducial = 0,
    reset_rf_for_each_step = 1,
    first_is_fiducial = 0,
    restrict_fiducialization = 0,
    reset_scattering_seed = 0,
&end
statistics:    ET:     00:00:01 CP:    0.86 BIO:0 DIO:0 PF:0 MEM:0
&bunched_beam
    bunch = {NULL},
    n_particles_per_bunch = 1,
    time_start = 0.000000000000000e+00,
    matched_to_cell = {NULL},
    emit_x = 0.000000000000000e+00,
    emit_nx = 0.000000000000000e+00,
    beta_x = 1.000000000000000e+00,
    alpha_x = 0.000000000000000e+00,
    eta_x = 0.000000000000000e+00,
    etap_x = 0.000000000000000e+00,
    emit_y = 0.000000000000000e+00,
    emit_ny = 0.000000000000000e+00,
    beta_y = 1.000000000000000e+00,
    alpha_y = 0.000000000000000e+00,
    eta_y = 0.000000000000000e+00,
    etap_y = 0.000000000000000e+00,
    use_twiss_command_values = 0,
    use_moments_output_values = 0,
    Po = 4.947171527086567e+03,
    sigma_dp = 0.000000000000000e+00,
    sigma_s = 0.000000000000000e+00,
    dp_s_coupling = 0.000000000000000e+00,
    emit_z = 0.000000000000000e+00,
    beta_z = 0.000000000000000e+00,
    alpha_z = 0.000000000000000e+00,
    momentum_chirp = 0.000000000000000e+00,
    one_random_bunch = 1,
    save_initial_coordinates = 1,
    limit_invariants = 0,
    symmetrize = 0,
    halton_sequence[0] = 0, 0, 0,
    halton_radix[0] = 0, 0, 0, 0, 0, 0,
    optimized_halton = 0,
    randomize_order[0] = 0, 0, 0,
    limit_in_4d = 0,
    enforce_rms_values[0] = 0, 0, 0,
    distribution_cutoff[0] = 2.000000000000000e+00, 2.000000000000000e+00, 2.000000000000000e+00,
    distribution_type[0] = gaussian, gaussian, gaussian,
    centroid[0] = 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00, 
        0.000000000000000e+00, 0.000000000000000e+00,
    first_is_fiducial = 0,
&end
statistics:    ET:     00:00:01 CP:    0.86 BIO:0 DIO:0 PF:0 MEM:0
&track
    center_on_orbit = 0,
    center_momentum_also = 1,
    offset_by_orbit = 0,
    offset_momentum_also = 1,
    soft_failure = 1,
    use_linear_chromatic_matrix = 0,
    longitudinal_ring_only = 0,
    ibs_only = 0,
    stop_tracking_particle_limit = -1,
    check_beam_structure = 0,
    interrupt_file = %s.interrupt,
&end
tracking step 1
generating bunch 1
tracking 1 particles
25 Apr 19 15:47:30: This step establishes energy profile vs s (fiducial beam).
25 Apr 19 15:47:30: Rf phases/references reset.
1 particles transmitted, total effort of 1 particle-turns
0 multipole kicks done

Dumping centroid data...done.
Post-tracking output completed.
Tracking step completed   ET:     00:00:01 CP:    0.88 BIO:0 DIO:0 PF:0 MEM:0


Saving lattice parameters to results/parameters.sdds...done.
Finished tracking.
End of input data encountered.
statistics:    ET:     00:00:01 CP:    0.88 BIO:0 DIO:0 PF:0 MEM:0
=====================================================================================
Thanks for using elegant.  Please cite the following reference in your publications:
  M. Borland, "elegant: A Flexible SDDS-Compliant Code for Accelerator Simulation,"
  Advanced Photon Source LS-287, September 2000.
If you use a modified version, please indicate this in all publications.
=====================================================================================
In [6]:
os.listdir("results")
Out[6]:
['beam.cen', 'beamline.mag', 'parameters.sdds', 'twiss.twi']

Reading beamline.mag

In [7]:
out = !sdds2stream results/beamline.mag -col=ElementName,s,Profile -pipe=out
In [8]:
out[:3]
Out[8]:
['_BEGIN_ 0.000000000000000e+00 0.000000000000000e+00',
 'W0 0.000000000000000e+00 0.000000000000000e+00',
 'W0 0.000000000000000e+00 0.000000000000000e+00']
In [9]:
print("\n".join(out[:3]))

_BEGIN_ 0.000000000000000e+00 0.000000000000000e+00
W0 0.000000000000000e+00 0.000000000000000e+00
W0 0.000000000000000e+00 0.000000000000000e+00
In [10]:
from io import StringIO
In [11]:
DATA = StringIO("\n".join(out))
In [12]:
df = pd.read_csv(DATA, names=['ElementName', 's', 'Profile'], delim_whitespace=True)
In [13]:
df.head(3)
Out[13]:
ElementName s Profile
0 _BEGIN_ 0.0 0.0
1 W0 0.0 0.0
2 W0 0.0 0.0
In [14]:
df.tail(3)
Out[14]:
ElementName s Profile
26130 W0 138.6 0.0
26131 W3 138.6 0.0
26132 W3 138.6 0.0
In [15]:
dim_s = hv.Dimension('s', unit='m', label="s") # range=(-80,0))
In [16]:
#hv.help(hv.Curve)
In [17]:
from bokeh.models import HoverTool
#hover = HoverTool(tooltips=[("Name", "@ElementName")])
hover = HoverTool(tooltips="@ElementName")

%opts Curve.mag [height=70 show_frame=False show_title=False \
                 xaxis=None yaxis=None tools=['xbox_zoom, xpan', hover]] (color='black', alpha=0.3)

mag = hv.Curve(df, kdims=dim_s, vdims=['Profile', 'ElementName'], group='mag')
mag
Out[17]:

Reading Twiss parameters file

In [18]:
#!sddsquery results/twiss.twi
In [19]:
out = !sdds2stream results/twiss.twi \
    -col=ElementName,s,betax,betay,alphax,alphay,etax,etay,pCentral0,xAperture,yAperture -pipe=out
DATA = StringIO("\n".join(out))
df = pd.read_csv(DATA, names=['ElementName','s','betax','betay','alphax','alphay',
                              'etax','etay','p','xAperture','yAperture'],
                   delim_whitespace=True)
In [20]:
df.tail(3)
Out[20]:
ElementName s betax betay alphax alphay etax etay p xAperture yAperture
2858 MA16 138.6 7.557469 1.516758 -1.087911 -1.412802 0.0 0.0 15.623737 0.008 0.008
2859 W0 138.6 7.557469 1.516758 -1.087911 -1.412802 0.0 0.0 15.623737 0.008 0.008
2860 W3 138.6 7.557469 1.516758 -1.087911 -1.412802 0.0 0.0 15.623737 0.008 0.008
In [21]:
mc = 0.511 # MeV/c
In [22]:
df['p'] = mc*df['p'] # MeV/c
df['xAperture'] = 1e3*df['xAperture'] # mm
df['yAperture'] = 1e3*df['yAperture'] # mm
In [23]:
%opts Curve.fx (color='red', alpha=0.7, line_width=3)
%opts Curve.fy (color='blue', alpha=0.7, line_width=3)
dim_p    = hv.Dimension('p', unit='MeV/c', range=(0,None))
dim_beta = hv.Dimension('beta', unit='m', label="β", range=(0,None))

s_p    = hv.Curve((df.s, df.p), kdims=dim_s, vdims=dim_p)
beta_x = hv.Curve((df.s, df.betax), label='βx', kdims=dim_s, vdims=dim_beta, group='fx')
beta_y = hv.Curve((df.s, df.betay), label='βy', kdims=dim_s, vdims=dim_beta, group='fy')
In [24]:
(s_p + beta_x*beta_y + mag).cols(1)
Out[24]:
In [25]:
#Dx = hv.Curve((df.s, df.etax), label='Dx', kdims='s (m)', vdims='D (m)', group='fx')
#Dy = hv.Curve((df.s, df.etay), label='Dy', group='fy')

#(Dx*Dy + mag).cols(1)
In [26]:
e_nx = 9.2 # normalized emittance [mm*mrad]
e_ny = e_nx

#e_x = e_nx * 1e-6 / gamma # geometric emittance [m]
#e_y = e_ny * 1e-6 / gamma

df['e_x'] = e_nx * 1e-6 / (df['p']/mc) # geometric emittance [m]
df['dp/p'] = 0
In [27]:
df['sigma_x'] = 1e3*np.sqrt(df['betax']*df['e_x'] + df['etax']*df['etax']*df['dp/p']*df['dp/p']) # mm
df['sigma_y'] = 1e3*np.sqrt(df['betay']*df['e_x'] + df['etay']*df['etay']*df['dp/p']*df['dp/p']) # mm
df.tail(3)
Out[27]:
ElementName s betax betay alphax alphay etax etay p xAperture yAperture e_x dp/p sigma_x sigma_y
2858 MA16 138.6 7.557469 1.516758 -1.087911 -1.412802 0.0 0.0 7.98373 8.0 8.0 5.888476e-07 0 2.109549 0.94506
2859 W0 138.6 7.557469 1.516758 -1.087911 -1.412802 0.0 0.0 7.98373 8.0 8.0 5.888476e-07 0 2.109549 0.94506
2860 W3 138.6 7.557469 1.516758 -1.087911 -1.412802 0.0 0.0 7.98373 8.0 8.0 5.888476e-07 0 2.109549 0.94506
In [28]:
%opts Curve.Aper (color='gray', alpha=0.5, line_width=3)

dim_sigma = hv.Dimension('sigma', unit='mm', label="2σ", range=(0,None))
dim_aper  = hv.Dimension('aper', unit='mm', label="Aperture")#, range=(0,None))

Ax =     hv.Curve((df.s,  df.xAperture), label='Aperture', kdims=dim_s, vdims=dim_sigma, group='Aper')
Ax_pos = Ax
Ax_neg = hv.Curve((df.s, -df.xAperture), label='Aperture', kdims=dim_s, vdims=dim_sigma, group='Aper')

Sx = hv.Curve((df.s, 2*df.sigma_x), label='2σx', kdims=dim_s, vdims=dim_sigma, group='fx')
Sy = hv.Curve((df.s, 2*df.sigma_y), label='2σy', kdims=dim_s, vdims=dim_sigma, group='fy')
In [29]:
(Sx*Sy*Ax + mag).cols(1)
Out[29]:

Beam centroid:

In [30]:
!sddsquery results/beam.cen

file results/beam.cen is in SDDS protocol version 1
description: centroid output--input: twiss.ele  lattice: machine.lte
contents: centroid output
data is little-endian binary

13 columns of data:
NAME            UNITS           SYMBOL          FORMAT          TYPE    FIELD  DESCRIPTION
                                                                        LENGTH
s               m               NULL            NULL            double  0       Distance
ElementName     NULL            NULL            %10s            string  0       Element name
ElementOccurence NULL            NULL            %6ld            long    0       Occurence of element
ElementType     NULL            NULL            %10s            string  0       Element-type name
Cx              m               <x>             NULL            double  0       x centroid
Cxp             NULL            <x'>            NULL            double  0       x' centroid
Cy              m               <y>             NULL            double  0       y centroid
Cyp             NULL            <y'>            NULL            double  0       y' centroid
Cs              m               <s>             NULL            double  0       mean distance traveled
Cdelta          NULL            <$gd$r>         NULL            double  0       delta centroid
Particles       NULL            NULL            NULL            long    0       Number of particles
pCentral        m$be$nc         p$bcen$n        NULL            double  0       Reference beta*gamma
Charge          C               NULL            NULL            double  0       Charge in the beam

2 parameters:
NAME                UNITS               SYMBOL              TYPE                DESCRIPTION
Step                NULL                NULL                long                Simulation step
SVNVersion          NULL                NULL                string              SVN version number
In [31]:
out = !sdds2stream results/beam.cen -col=ElementName,s,Cx,Cy,pCentral -pipe=out
DATA = StringIO("\n".join(out))
df_cen = pd.read_csv(DATA, names=['ElementName','s','Cx','Cy','p'],
                   delim_whitespace=True)
In [32]:
df_cen.tail(3)
Out[32]:
ElementName s Cx Cy p
2858 MA16 138.6 0.0 0.0 15.623737
2859 W0 138.6 0.0 0.0 15.623737
2860 W3 138.6 0.0 0.0 15.623737
In [33]:
df_cen['Cx'] = df_cen['Cx']*1e3 # m -> mm
df_cen['Cy'] = df_cen['Cy']*1e3 # m -> mm
df_cen['p']  = df_cen['p']*0.511 # gamma*beta -> MeV/c
In [34]:
dim_x = hv.Dimension('x', unit='mm', label="x", range=(-16,+16))
dim_y = hv.Dimension('y', unit='mm', label="y", range=(-16,+16))
dim_p = hv.Dimension('p', unit='MeV/c', label="p", range=(0,None))
In [35]:
s_p = hv.Curve((df_cen.s, df_cen.p), kdims=dim_s, vdims=dim_p)

Cx  = hv.Spread((df_cen.s, df_cen.Cx, 2*df.sigma_x), label='x', \
               kdims=dim_s, vdims=[dim_x,dim_aper])
Cy  = hv.Spread((df_cen.s, df_cen.Cy, 2*df.sigma_y), label='y', \
               kdims=dim_s, vdims=[dim_y,dim_aper])
In [36]:
(s_p + Ax_pos*Ax_neg*Cx*Cy + mag).cols(1)
Out[36]:

Misc tests

In [ ]:
In [ ]: