search

MethodicalAcceleratorDesign / MAD-X

MAD-X repository
Other
49 stars 39 forks source link

angle option in multipole not documented and new proposal #911

Open rdemaria opened 4 years ago

rdemaria commented 4 years ago

The documentation of the angle attribute is not implemented. The effect is a bit counter-intuitive

mm: multipole, angle=1e-20, tilt=0, knl:={1}, ksl:={0}; ! hkicker negative
mm: multipole, angle=1e-20, tilt=pi/2, knl:={1}, ksl:={0};  !  vkicker negative
mm: multipole, angle=1e-20, tilt=0, knl:={0}, ksl:={1}; ! vertical dipole
mm: multipole, angle=1e-20, tilt=0, knl:={1}, ksl:={1}; ! hkicker negative kick + vertical dipole

I wonder if one could change the implementation and simplify it in such a way:

mm: multipole, angle=1e-20, tilt=0, knl:={1}, ksl:={0}; ! hkicker negative
mm: multipole, angle=1e-20, tilt=pi/2, knl:={1}, ksl:={0};  ! vkicker negative
mm: multipole, angle=1e-20, tilt=0, knl:={0}, ksl:={1}; ! vkicker positive
mm: multipole, angle=1e-20, tilt=0, knl:={1}, ksl:={1}; ! hkicker negative + vkicker positive

In the second implementation angle and tilt fully specify the change of reference frame and knl,ksl only the multipole fields in the tilted frame.

I checked the implementation in twiss only, I have not checked whether track and survey are consistent with twiss.

rdemaria commented 4 years ago

Presently, the angle property of an sbend is not transferred by makethin. The proposal above also allows to pass the "angle" property into the slices in a consistent manner.

rdemaria commented 1 year ago

Related to this, in case an sbend has k0=!0 and angle=0, makethin should add angle=1e-30 (whatever small number testing !=0 to disable the assumption angle==k0.

rdemaria commented 1 year ago

Reviewing the code in track I noticed an inconsistency when angle!=0.

The equations for the thin dipole are:

\begin{aligned}
H & = p_\tau - (1+h_x x -h_y y) (1+\delta) + k_0 \left(x + \frac{h_x x^2}{2}\right) + \hat k_0 \left(y + \frac{h_y y^2}{2}\right)\\
\Delta p_x &= - l ( (k_0 - h_x) - h_x\delta + k_0 h_x x)  \\
\Delta p_y &= +l ( (\hat k_0 - h_y) -h_y\delta + \hat k_0 h_y y)  \\
\Delta \tau &= l \frac{h_x x - h_y y}{\beta}
 \end{aligned}

In track

if (an .ne. 0) f_errors(0) = f_errors(0) + normal(0) - an
dbr = bvk * f_errors(0) !field(1,0)
dbi = bvk * f_errors(1) !field(2,0)
dipr = bvk * normal(0) !vals(1,0)
dipi = bvk * skew(0)   !vals(2,0)

In the weak focusing is calculated using 

DXT(:ktrack) = DXT(:ktrack) + (dipr-dbr)*dipr*TRACK(1,:ktrack)/elrad
DYT(:ktrack) = DYT(:ktrack) + (dipi-dbi)*dipi*TRACK(3,:ktrack)/elrad

but (dipr-dbr) != an when there are errors. In addition:

     track(2,jtrk) = track(2,jtrk) - (dbr + dxt(jtrk) - dipr * (ttt - one))
     track(4,jtrk) = track(4,jtrk) + (dbi + dyt(jtrk) - dipi * (ttt - one))
     track(5,jtrk) = track(5,jtrk) - &
          (dipr*track(1,jtrk) - dipi*track(3,jtrk)) *   &
          ((one + betas*track(6,jtrk))/ttt) * beti

where (ttt-one)=$\delta$ which is not correct when angle!=0.

rdemaria commented 1 year ago

Also on twiss, the weak focusing does not scale with h*k0, but with h**2

https://github.com/MethodicalAcceleratorDesign/MAD-X/blob/master/src/twiss.f90#L4155

xksq = h**2 + sk1
rdemaria commented 1 year ago

Adding script from @giadarol reporting convergence issue with sbend and multipoles

import numpy as np
import matplotlib.pyplot as plt

config = {
    'on_translation': 1,
    'on_rotation': 1,
    'on_k1': 1,
    'on_shifts': 1,
    'kill_orbit': False,
    'n_slices_bend': 100
}

sequence_src = '''
none = 0;
start_check: marker,l= 0,kmax= 0,kmin= 0,calib= 0,polarity= 0,type,apertype="circle",aperture={ 0},aper_offset={ 0},aper_tol={ 0},aper_vx={ -1},aper_vy={ -1},slot_id= 0,assembly_id= 0,mech_sep= 0,v_pos= 0,magnet= 0,model= -1,method= -1
,exact= -1,nst= -1,fringe= 0,bend_fringe=false,kill_ent_fringe=false,kill_exi_fringe=false,dx= 0,dy= 0,ds= 0,dtheta= 0,dphi= 0,dpsi= 0,aper_tilt= 0,comments;
ptb_b0pf: translation,dx= -0.03022642196;
prb_b0pf: yrotation,angle= 0.02670439316;
b0pf: sbend,l= 1.2,k0= -0.008660083518, k1:= -0.05940545468*on_k1, angle=1e-30;
pte_b0pf: translation,dx= -0.0007490490899;
pre_b0pf: yrotation,angle= -0.025;
ptb_sol_f: translation,dx= -0.04999479183;
prb_sol_f: yrotation,angle= 0.025;
sol_half: solenoid,l= 2;
star_detect_f: sol_half;
pre_sol_f: yrotation,angle= -0.025;
ip6w: marker;
hsr41: sequence, l = 8;
start_check, at = 0;
ptb_b0pf, at = 0.8851317104, dx = -0.03022642196 ;
prb_b0pf, at = 0.8851317104;
b0pf, at = 1.48513171;
pte_b0pf, at = 2.08513171, dx = -0.0007490490899 ;
pre_b0pf, at = 2.08513171;
ptb_sol_f, at = 5.88756965, dx = -0.04999479183 ;
prb_sol_f, at = 5.88756965;
star_detect_f, at = 6.88756965;
pre_sol_f, at = 7.88756965;
ip6w, at = 7.88756965;
endsequence;
'''

new_lines = []
for ll in sequence_src.split('\n'):
    if 'translation' in ll:
        ll = ll.replace(';', "*on_translation;")
        ll = ll.replace('dx=', "dx:=")
    elif 'yrotation' in ll:
        ll = ll.replace(';', "*on_rotation;")
        ll = ll.replace('angle=', "angle:=")
    elif 'dx = ' in ll:
        ll = ll.replace(';', "*on_shifts;")
        ll = ll.replace('dx = ', "dx:=")
    new_lines.append(ll)
new_lines.append(f"on_translation={config['on_translation']};")
new_lines.append(f"on_rotation={config['on_rotation']};")
new_lines.append(f"on_k1={config['on_k1']};")
new_lines.append(f"on_shifts={config['on_shifts']};")
sequence_src = '\n'.join(new_lines)

tw_init = {'betx': 84.02557752667167,
            'alfx': 14.23087869857309,
            'bety': 737.8650722816444,
            'alfy': 56.03568671202535,
            'x': 0.014288215914637994,
            'px': -0.009838644402766588,
            'y': 0.0,
            'py': 0.0}

if config['kill_orbit']:
    tw_init['x'] = 0
    tw_init['px'] = 0
    tw_init['y'] = 0
    tw_init['py'] = 0

mad_thick = Madx()
mad_thick.input(sequence_src)
mad_thick.beam(energy=41.0,particle='antiproton')
mad_thick.use(sequence='hsr41')
tw_thick = mad_thick.twiss(**tw_init)

mad_thin = Madx()
mad_thin.input(sequence_src)
mad_thin.beam(energy=41.0,particle='antiproton')
mad_thin.use(sequence='hsr41')

mad_thin.input(f'''
select, flag=makethin, clear;
select, flag=makethin, class=quadrupole, slice = 20, thick=false;
select, flag=makethin, class=sextupole, slice = 1, thick=true;
select, flag=makethin, class=sbend, slice=0;

select, flag=makethin, pattern=b0pf, slice={config['n_slices_bend']}, thick=false;

makethin, sequence=hsr41, style=teapot, makedipedge=false;
''')
mad_thin.use('hsr41')

tw_thin = mad_thin.twiss(**tw_init)

bety_thin_on_thick_check = np.interp(tw_thick.s, tw_thin['s'], np.sqrt(tw_thin['bety']))**2
alfy_thin_on_thick_check = np.interp(tw_thick.s, tw_thin['s'], tw_thin['alfy'])

plt.close('all')
plt.figure()
ax1 = plt.subplot(311)
plt.plot(tw_thick['s'], tw_thick['bety']/bety_thin_on_thick_check -1, '.-')
plt.ylabel(r'$\Delta \beta_y / \beta_y$')

plt.subplot(312, sharex=ax1)
plt.plot(tw_thick['s'], tw_thick['alfy']/alfy_thin_on_thick_check -1, '.-')
plt.ylabel(r'$\Delta \alpha_y / \alpha_y$')

ax2 = plt.subplot(313, sharex=ax1)
plt.plot(tw_thick['s'], tw_thick.x)
plt.plot(tw_thin['s'], tw_thin.x)
plt.ylabel(r'$x$')

s_b0pf = tw_thick.dframe().loc[tw_thick.name=='b0pf:1', 's'].values[0]

plt.axvline(s_b0pf, color='r')
plt.subplots_adjust(left=.15)
plt.show()
HelmutCERN commented 1 year ago

Dear Riccardo

from a first look : There is some misuse of makethin here, ( not respecting the makethin documentation ) : Thick slicing sextuples does not exist. If turning off sextuple slicing is wanted do with "slice = 0”. sbend,slice=0. turns off sbend slicing. Why makedipedge is turned off ? Rather special short sequence : The start_check marker has all kind of attributes that have no effect. The sequence has a solenoid, that will be sliced here with a single thin slice, is that intended ?

Cheers, Helmut

On 1 Mar, 2023, at 10:27, Riccardo De Maria @.***> wrote:

Adding script from @giadarol https://github.com/giadarol reporting convergence issue with sbend and multipoles

import numpy as np import matplotlib.pyplot as plt

config = { 'on_translation': 1, 'on_rotation': 1, 'on_k1': 1, 'on_shifts': 1, 'kill_orbit': False, 'n_slices_bend': 100 }

sequence_src = ''' none = 0; start_check: marker,l= 0,kmax= 0,kmin= 0,calib= 0,polarity= 0,type,apertype="circle",aperture={ 0},aper_offset={ 0},aper_tol={ 0},aper_vx={ -1},aper_vy={ -1},slot_id= 0,assembly_id= 0,mech_sep= 0,v_pos= 0,magnet= 0,model= -1,method= -1 ,exact= -1,nst= -1,fringe= 0,bend_fringe=false,kill_ent_fringe=false,kill_exi_fringe=false,dx= 0,dy= 0,ds= 0,dtheta= 0,dphi= 0,dpsi= 0,aper_tilt= 0,comments; ptb_b0pf: translation,dx= -0.03022642196; prb_b0pf: yrotation,angle= 0.02670439316; b0pf: sbend,l= 1.2,k0= -0.008660083518, k1:= -0.05940545468*on_k1, angle=1e-30; pte_b0pf: translation,dx= -0.0007490490899; pre_b0pf: yrotation,angle= -0.025; ptb_sol_f: translation,dx= -0.04999479183; prb_sol_f: yrotation,angle= 0.025; sol_half: solenoid,l= 2; star_detect_f: sol_half; pre_sol_f: yrotation,angle= -0.025; ip6w: marker; hsr41: sequence, l = 8; start_check, at = 0; ptb_b0pf, at = 0.8851317104, dx = -0.03022642196 ; prb_b0pf, at = 0.8851317104; b0pf, at = 1.48513171; pte_b0pf, at = 2.08513171, dx = -0.0007490490899 ; pre_b0pf, at = 2.08513171; ptb_sol_f, at = 5.88756965, dx = -0.04999479183 ; prb_sol_f, at = 5.88756965; star_detect_f, at = 6.88756965; pre_sol_f, at = 7.88756965; ip6w, at = 7.88756965; endsequence; '''

new_lines = [] for ll in sequence_src.split('\n'): if 'translation' in ll: ll = ll.replace(';', "on_translation;") ll = ll.replace('dx=', "dx:=") elif 'yrotation' in ll: ll = ll.replace(';', "on_rotation;") ll = ll.replace('angle=', "angle:=") elif 'dx = ' in ll: ll = ll.replace(';', "*on_shifts;") ll = ll.replace('dx = ', "dx:=") new_lines.append(ll) new_lines.append(f"on_translation={config['on_translation']};") new_lines.append(f"on_rotation={config['on_rotation']};") new_lines.append(f"on_k1={config['on_k1']};") new_lines.append(f"on_shifts={config['on_shifts']};") sequence_src = '\n'.join(new_lines)

tw_init = {'betx': 84.02557752667167, 'alfx': 14.23087869857309, 'bety': 737.8650722816444, 'alfy': 56.03568671202535, 'x': 0.014288215914637994, 'px': -0.009838644402766588, 'y': 0.0, 'py': 0.0}

if config['kill_orbit']: tw_init['x'] = 0 tw_init['px'] = 0 tw_init['y'] = 0 tw_init['py'] = 0

mad_thick = Madx() mad_thick.input(sequence_src) mad_thick.beam(energy=41.0,particle='antiproton') mad_thick.use(sequence='hsr41') tw_thick = mad_thick.twiss(**tw_init)

mad_thin = Madx() mad_thin.input(sequence_src) mad_thin.beam(energy=41.0,particle='antiproton') mad_thin.use(sequence='hsr41')

mad_thin.input(f''' select, flag=makethin, clear; select, flag=makethin, class=quadrupole, slice = 20, thick=false; select, flag=makethin, class=sextupole, slice = 1, thick=true; select, flag=makethin, class=sbend, slice=0;

select, flag=makethin, pattern=b0pf, slice={config['n_slices_bend']}, thick=false;

makethin, sequence=hsr41, style=teapot, makedipedge=false; ''') mad_thin.use('hsr41')

tw_thin = mad_thin.twiss(**tw_init)

bety_thin_on_thick_check = np.interp(tw_thick.s, tw_thin['s'], np.sqrt(tw_thin['bety']))**2 alfy_thin_on_thick_check = np.interp(tw_thick.s, tw_thin['s'], tw_thin['alfy'])

plt.close('all') plt.figure() ax1 = plt.subplot(311) plt.plot(tw_thick['s'], tw_thick['bety']/bety_thin_on_thick_check -1, '.-') plt.ylabel(r'$\Delta \beta_y / \beta_y$')

plt.subplot(312, sharex=ax1) plt.plot(tw_thick['s'], tw_thick['alfy']/alfy_thin_on_thick_check -1, '.-') plt.ylabel(r'$\Delta \alpha_y / \alpha_y$')

ax2 = plt.subplot(313, sharex=ax1) plt.plot(tw_thick['s'], tw_thick.x) plt.plot(tw_thin['s'], tw_thin.x) plt.ylabel(r'$x$')

s_b0pf = tw_thick.dframe().loc[tw_thick.name=='b0pf:1', 's'].values[0]

plt.axvline(s_b0pf, color='r') plt.subplots_adjust(left=.15) plt.show() — Reply to this email directly, view it on GitHub https://github.com/MethodicalAcceleratorDesign/MAD-X/issues/911#issuecomment-1449679439, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEUCRLRA5WECHMAOY4MJNGTWZ4I67ANCNFSM4ODJ54RQ. You are receiving this because you are subscribed to this thread.

giadarol commented 1 year ago

Hi Helmut, thanks for your comments. I cleaned up the script (which indeed was a bit messy). In particular:

Here the cleaned up version:


from cpymad.madx import Madx
import numpy as np
import matplotlib.pyplot as plt

config = {
    'on_translation': 1,
    'on_rotation': 1,
    'on_k1': 1,
    'on_shifts': 1,
    'kill_orbit': False,
    'n_slices_bend': 1000
}

sequence_src = '''

option, thin_cf=false;

none = 0;
start_check: marker,l= 0,kmax= 0,kmin= 0,calib= 0,polarity= 0,type,apertype="circle",aperture={ 0},aper_offset={ 0},aper_tol={ 0},aper_vx={ -1},aper_vy={ -1},slot_id= 0,assembly_id= 0,mech_sep= 0,v_pos= 0,magnet= 0,model= -1,method= -1
,exact= -1,nst= -1,fringe= 0,bend_fringe=false,kill_ent_fringe=false,kill_exi_fringe=false,dx= 0,dy= 0,ds= 0,dtheta= 0,dphi= 0,dpsi= 0,aper_tilt= 0,comments;
ptb_b0pf: translation,dx= -0.03022642196;
prb_b0pf: yrotation,angle= 0.02670439316;
b0pf: sbend,l= 1.2,k0= -0.008660083518, k1:= -0.05940545468*on_k1 + corr_k1, angle=1e-30;
pte_b0pf: translation,dx= -0.0007490490899;
pre_b0pf: yrotation,angle= -0.025;
ptb_sol_f: translation,dx= -0.04999479183;
prb_sol_f: yrotation,angle= 0.025;
sol_half: drift,l= 2;
star_detect_f: sol_half;
pre_sol_f: yrotation,angle= -0.025;
ip6w: marker;
hsr41: sequence, l = 8;
start_check, at = 0;
ptb_b0pf, at = 0.8851317104, dx = -0.03022642196 ;
prb_b0pf, at = 0.8851317104;
b0pf, at = 1.48513171;
pte_b0pf, at = 2.08513171, dx = -0.0007490490899 ;
pre_b0pf, at = 2.08513171;
ptb_sol_f, at = 5.88756965, dx = -0.04999479183 ;
prb_sol_f, at = 5.88756965;
star_detect_f, at = 6.88756965;
pre_sol_f, at = 7.88756965;
ip6w, at = 7.88756965;
endsequence;
'''

new_lines = []
for ll in sequence_src.split('\n'):
    if 'translation' in ll:
        ll = ll.replace(';', "*on_translation;")
        ll = ll.replace('dx=', "dx:=")
    elif 'yrotation' in ll:
        ll = ll.replace(';', "*on_rotation;")
        ll = ll.replace('angle=', "angle:=")
    elif 'dx = ' in ll:
        ll = ll.replace(';', "*on_shifts;")
        ll = ll.replace('dx = ', "dx:=")
    new_lines.append(ll)
new_lines.append(f"on_translation={config['on_translation']};")
new_lines.append(f"on_rotation={config['on_rotation']};")
new_lines.append(f"on_k1={config['on_k1']};")
new_lines.append(f"on_shifts={config['on_shifts']};")
sequence_src = '\n'.join(new_lines)

tw_init = {'betx': 84.02557752667167,
            'alfx': 14.23087869857309,
            'bety': 737.8650722816444,
            'alfy': 56.03568671202535,
            'x': 0.014288215914637994,
            'px': -0.009838644402766588,
            'y': 0.0,
            'py': 0.0}

if config['kill_orbit']:
    tw_init['x'] = 0
    tw_init['px'] = 0
    tw_init['y'] = 0
    tw_init['py'] = 0

mad_thick = Madx()
mad_thick.input(sequence_src)
mad_thick.beam(energy=41.0,particle='antiproton')
mad_thick.use(sequence='hsr41')
tw_thick = mad_thick.twiss(**tw_init)

mad_thin = Madx()
mad_thin.input(sequence_src)
mad_thin.beam(energy=41.0,particle='antiproton')
mad_thin.use(sequence='hsr41')

mad_thin.input(f'''
select, flag=makethin, clear;
select, flag=makethin, class=quadrupole, slice = 20, thick=false;
select, flag=makethin, pattern=b0pf, slice={config['n_slices_bend']}, thick=false;

makethin, sequence=hsr41, style=teapot, makedipedge=true;
''')
mad_thin.use('hsr41')

tw_thin = mad_thin.twiss(**tw_init)

bety_thin_on_thick_check = np.interp(tw_thick.s, tw_thin['s'], np.sqrt(tw_thin['bety']))**2
alfy_thin_on_thick_check = np.interp(tw_thick.s, tw_thin['s'], tw_thin['alfy'])

plt.close('all')
plt.figure()
ax1 = plt.subplot(211)
plt.plot(tw_thick['s'][:-2], (tw_thick['bety']/bety_thin_on_thick_check -1)[:-2], '.-')
plt.ylabel(r'$\Delta \beta_y / \beta_y$')

ax2 = plt.subplot(212, sharex=ax1)
plt.plot(tw_thick['s'], tw_thick.x)
plt.plot(tw_thin['s'], tw_thin.x)
plt.ylabel(r'$x$')

s_b0pf = tw_thick.dframe().loc[tw_thick.name=='b0pf:1', 's'].values[0]

plt.axvline(s_b0pf, color='r')
plt.subplots_adjust(left=.15)
plt.show()
giadarol commented 1 year ago

I forgot to mention that these changes have no effect on the observed problem