Exact tracking of magnetic chicane line

A. Petrenko, N. Charitonidis (CERN, 2017)

This notebook can be executed at the MS Azure Cloud service (Exact_Chicate.ipynb file).

from math import *
import holoviews as hv
import numpy as np

hv.notebook_extension('bokeh')

def BLine(L=2, D=1, R=10):
    #    L    R
    return [
        [5*D, 0],
        [L  , R],
        [D  , 0],
        [L  ,-R],
        [2*D, 0],
        [L  ,-R],
        [D  , 0],
        [L  , R],
        [5*D, 0],
    ]

def track(R=3, yp=0.0):
    y  = 0.0;
    #yp = 0.05; # dy/dx
    x  = 0;

    x_values=[]
    y_values=[]

    split = 10 # split every element into that number of segments

    for itm in BLine(R=R):
        r=itm[1]; l=itm[0]

        for i in range(0,split): #splitting in segments with parameter "split"    

            x_values.append(x)
            y_values.append(y)

            if r==0:
                y=y+yp*(l/split) # exact solution for each segment
                x=x+(l/split)
            else:
                angle=atan(yp) #Angle of trajectory 

                nx=cos(angle-(pi/2)) # normal vector coordinates pointing to the center of bend. circle 
                ny=sin(angle-(pi/2))

                x0=x+nx*r
                y0=y+ny*r

                x=x+(l/split)

                ya=y0+sqrt(r*r-(x-x0)*(x-x0))  # positive sqrt solution
                yb=y0-sqrt(r*r-(x-x0)*(x-x0))  # negative sqrt solution

                if abs(y-ya) > abs(y-yb):
                    y=yb
                else: 
                    y=ya

                yp=(x-x0)/(y0-y)

    x_values.append(x)
    y_values.append(y)
    #print(x_values) ; print(y_values)
    
    return (x_values, y_values)

#hv.Curve( track(R=10, yp=0.02) ) # static plot

Magnet layout:

x = 0
positive_bends = []
negative_bends = []

def rectangle(x=0, y=0, width=1, height=12):
    return np.array([(x,y-height/2), (x+width, y-height/2), (x+width, y+height/2), (x, y+height/2)])

for itm in BLine():
    r=itm[1]; l=itm[0]

    if r>0:
        positive_bends.append(rectangle(x=x,width=l))
    if r<0:
        negative_bends.append(rectangle(x=x,width=l))
        
    x = x + l

p_bends = hv.Polygons(positive_bends)(style={'fill_color': '#0000ff', 'fill_alpha': 0.2, 'line_alpha': 0})
n_bends = hv.Polygons(negative_bends)(style={'fill_color': '#ff0000', 'fill_alpha': 0.2, 'line_alpha': 0})

#n_bends * p_bends

Interactive plot:

R_points  = np.linspace(2.3, 6.3, 9) # m
yp_points = np.linspace(0, 0.1, 6) # rad

keys = [(R, yp) for R in R_points for yp in yp_points]

items = [(k, hv.Curve(track(*k)) * n_bends * p_bends ) for k in keys]

%opts Curve [width=500 height=300] 
hv.HoloMap(items, kdims = ['Bending radius', "Incoming particle angle"])

Bending radius:

Incoming particle angle: