BBDiagnostics: Computer Class MatlabAT

1
BB-Diagnostics Computer Class (Matlab/AT)
Using Accelerator Tracking Code AT, Middle Layer
Software, and Response Matrix Fitting

• Greg Portmann, Christoph Steier, and James
  Safranek

• Monday Intro to the Middle Layer and AT
• Orbit errors
• Tuesday Gradient errors, beta function, phase
  advance
• Wednesday Gradient errors, LOCO
• Thursday Detuning with amplitude, decoherence,
  resonances
• (Christoph Steier)
• Friday no computer class

2
Starting Matlab

• Start Matlab 6.5 (AT and LOCO require at least
  Matlab version 6.0)
• Matlab 7.1 is the latest version.
• Set the matlab path for AT, MiddleLayer, and
  LOCO
• cd c\bbdiag\acceleratorcontrol
• setpathals
• or
• setpathspear3
• setpathxray
• setpathvuv
• Change to the class examples directory
• cd docs
• Note When in doubt, setpathals will also
  restore the default lattice.

3
Monday AT Lattice Model Orbits

• Plot lattice (AT)
• intlat
• Edit the AT lattice file
• edit alslat
• Plot the orbit with plotfamily (Middle Layer)
  
• plotfamily
• Menu Graph 1 -gt Add to plot -gt Draw Lattice
• Switch synchrotron radiation and RF cavity on
  off (Middle Layer)
• setcavity on
• setradiation on
• Update plotfamily
• setcavity off
• setradiation off
• Update plotfamily
• Change the first corrector magnet in sector 6
  by 1 amp (Middle Layer)

4
Monday Closed Orbits

• Get the orbit, beta, and phase at all elements
  in the AT model (Middle Layer)
• x, y, BPMs modeltwiss('ClosedOrbit')
• BetaX, BetaY, BPMs modeltwiss('Beta')
• MuX, MuY modeltwiss('Phase')
• Plot the orbit vs. s-position
• figure
• subplot(2,1,1)
• plot(BPMs, x)
• subplot(2,1,2)
• plot(BPMs, y)
• Plot the orbit vs. phase
• figure
• subplot(2,1,1)
• plot(MuX, x)
• subplot(2,1,2)
• plot(MuY, y)

5
Monday Closed Orbits

• Try these commands
• modelbeta
• modelbeta('BPMx')
• modelbeta('QF')
• modeldisp
• modeldisp('BEND')
• FractionalTune, IntegerTune modeltune
• Get closed orbit directly from AT
• help findorbit4
• global THERING
• x findorbit4(THERING, 0.0, 1length(THERING))
• s findspos(THERING, 1length(THERING))
  
• Plot the orbit vs. s-position
• subplot(2,1,1)

6
Tuesday Gradient Errors
Add the perturbed beta to figure(1) figure(1)
subplot(2,2,1) hold on plot(MuX0, BetaX1,
'r') subplot(2,2,3) hold on plot(MuY0,
BetaY1, 'r') subplot(2,2,2) hold on
plot(MuX0, MuX1/(2pi) , 'r') subplot(2,2,4)
hold on plot(MuY0, MuY1/(2pi) , 'r')
Plot beta beat figure(2) subplot(2,1,1)
plot(s, BetaX1./BetaX0) title('Beta Beat from
the Nominal Model') subplot(2,1,2) plot(s,
BetaY1./BetaY0) Restore the
lattice setsp('QF', qf, 7 1)

• Plot the beta function and phase advance for
  the
• nominal model and a model with a gradient
  error
• Get beta phase at all elements in the AT
  model (Middle Layer)
• Tune0 gettune
• BetaX0, BetaY0, s modeltwiss('Beta')
• MuX0, MuY0 modeltwiss('Phase')
• Plot beta vs. position
• figure(1) clf
• subplot(2,2,1) plot(MuX0, BetaX0, 'b')
  ylabel('Beta X')
• subplot(2,2,3) plot(MuY0, BetaY0, 'b')
  ylabel('Beta Y')
• subplot(2,2,2) plot(MuX0, MuX0/(2pi) , 'b')
  ylabel('Phase X')
• subplot(2,2,4) plot(MuY0, MuY0/(2pi) , 'b')
  ylabel('Phase Y')
• Perturb the lattice at 1 quadrupole, QF(7,1)
• qf getsp('QF', 7 1)
• setsp('QF', 1.05qf, 7 1)

7
Tuesday Quadrupole Centers

• Beam based measurement of quadrupole centers
• quadcenter
• ? QF
• ? QF(7,1)
• ? Vertical only
• Note your data got put in c\bbdiag\machine\als
  data\HighTune\QMS

8
Wednesday FFT Analyze

• Simulate a phase advance measurement
• Get the orbit, beta, and phase at all elements
  in the AT model (Middle Layer)
• Tune0 gettune
• BetaX0, BetaY0, s modeltwiss('Beta',
  'BPMx')
• MuX0, MuY0 modeltwiss('Phase', 'BPMx')
• Starting condition or tracking (0.1mm)
• X0 0.0001 0 0.0001 0 0 0
• Track for 1024 turns
• global THERING
• X1 ringpass(THERING, X0, 1024)
• size(X1)
• Track coordinates for every turn along the
  ring (to all BPMs)
• BPMindex findcells(THERING, 'FamName', 'BPM')
• BPM findorbit4(THERING, 0.0, BPMindex)
• X2 linepass(THERING, X1, BPMindex)
• size(X2)
• Recover matrix structure (turns x BPM)

9
Wednesday FFT Analyze

• Calculate fractional tunes (interpolating FFT,
  sine window)
• nux, nuy, ax, ay findfreq(BPMx, BPMy)
• Calculate phase at every BPM
• (integral convolution with sine and cosine
  trajectories)
• MuX, MuY calcphase(nux, nuy, BPMx, BPMy)
• Calcphase asks for a frequency, typically just
  accepting
• the precalculated result is fine.
• Compare the 'measured' phase advance with the
  computed nominal one.
• DeltaMuX MuX()-MuX0/(2pi)
• DeltaMuY MuY()-MuY0/(2pi)
• figure(2)
• subplot(2,1,1)
• plot(MuX0, DeltaMuX-DeltaMuX(1), '.-b')
• subplot(2,1,2)
• plot(MuY0, DeltaMuY-DeltaMuY(1), '.-b')

10
Wednesday FFT Analyze

• Add noise to the BPM data and recalculation the
  phase
• BPMxNoise BPMx 5e-6randn(size(BPMx))
• BPMyNoise BPMy 5e-6randn(size(BPMy))
• Calculate fractional tunes (interpolating FFT,
  sine window)
• nux, nuy, ax, ay findfreq(BPMxNoise ,
  BPMyNoise)
• Calculate phase at every BPM
• (integral convolution with sine and cosine
  trajectories)
• MuXnoise, MuYnoise calcphase(nux, nuy,
  BPMxNoise, BPMyNoise)
• Calcphase asks for a frequency, typically just
  accepting
• the precalculated result is fine.
• Compare the 'measured' phase advance with the
  computed nominal one.
• DeltaMuX MuXnoise()-MuX0/(2pi)
• DeltaMuY MuYnoise()-MuY0/(2pi)
• figure(2)
• subplot(2,1,1) hold on

11
Wednesday LOCO Optics Analysis

• Calibrate/control optics using orbit response
  matrix
• Determine quadrupole gradients
• Correcte coupling
• Calibrate BPM gains, steering magnets
• Measure local chromaticity and transverse
  impedance

• New MATLAB version of code
• rewritten from FORTRAN
• linked to control system
• linked to AT simulator

12
Wednesday LOCO - locogui

• Start the LOCO graphical user interface (GUI)
• Look at an ALS LOCO run
• locogui
• ? open
• ? c\bbdiag\machine\alsdata\HighTune\LOCO\20
  05-10-16

13
Wednesday LOCO measure the input data

• Use the ALS ring
• setpathals
• Perturb the lattice at 1 quadrupole, QF(1,1)
• qf getsp('QF', 1 1)
• qf_k0 getsp('QF', 1 1, 'Physics')
  2.2371
• setsp('QF', 1.05qf, 1 1)
• qf_k1 getsp('QF', 1 1, 'Physics')
  2.3432
• DeltaK qf_k1 - qf_k0
  0.1060
• Note 1 micron of noise to the model BPMs
• Get LOCO data
• measlocodata
• ? Use the default model
• ? Use the default directory
• ? Dispersion
• ? BPM response matrix
• ? BPM sigma
• ? Dont correct the orbit

14
Wednesday LOCO Find a quadrupole error

• Open the GUI run
• locogui
• ? Inputs
• ? Turn off BPM coupling
• ? Turn off Correcor Magnet coupling
• ? Turn off Include off diagonal
  response matrix terms
• Start
• Result After 3 iterations
• FitValues 2.3433
• 2.2370
• 2.2372
• 2.2370

15
Thursday Detuning with amplitude

• Start Matlab 6.5 (newer versions 7.0 and 7.1
  will work as well)
• cd c\bbdiag\acceleratorcontrol
• setpathals
• global THERING
• Starting condition for tracking (100 micron x,
  100 micron y)
• X00.000100.0001000
• Track for 1024 turns
• X1ringpass(THERING,X0,1024)
• Plot tracking results
• figuresubplot(2,1,1)plot(X1(1,))subplot(2,1,2)
  plot(X1(3,))
• Calculate FFT
• fftxabs(fft(X1(1,))) fftyabs(fft(X1(3,)))
• Frequency/Tune vector
• nu01023/1023
• Plot FFT (x/y) - tunes
• fffiguresubplot(2,1,1)plot(nu,fftx/max(fftx))h
  old on
• subplot(2,1,2)plot(nu,ffty/max(ffty))hold on
• Find tunes and print result to screen
• maxx,kxmax(fftx(2512))

16
Thursday Detuning with amplitude

• Plot in tunespace
• ff2 figureplot(nux,nuy,'o')axis(0 0.26 0
  0.21)hold on
• Repeat with larger amplitude (1 mm x, 100
  micron y)
• X00.00100.0001000X1ringpass(THERING,X0,1
  024)
• fftxabs(fft(X1(1,)))fftyabs(fft(X1(3,)))
• figure(ff)subplot(2,1,1)plot(nu,fftx/max(fftx),'
  r')
• subplot(2,1,2)plot(nu,ffty/max(ffty),'r')
• maxx,kxmax(fftx(2512))maxy,kymax(ffty(251
  2))
• nux(2)nu(kx1)nuy(2)nu(ky1)
• fprintf('nu_x g, nu_y g \n',nu(kx1),nu(ky1
  ))
• figure(ff2)plot(nux,nuy,'o')axis(0 0.26 0
  0.21)
• Repeat again with larger amplitude (5 mm x, 100
  micron y) tuneshift with
• amplitude in the simplest case scales like
  initial position2
• X00.00500.0001000
• X1ringpass(THERING,X0,1024)
• fftxabs(fft(X1(1,)))fftyabs(fft(X1(3,)))
• figure(ff)subplot(2,1,1)plot(nu,fftx/max(fftx),'
  g')
• subplot(2,1,2)plot(nu,ffty/max(ffty),'g')
• maxx,kxmax(fftx(2512))maxy,kymax(ffty(251
  2))

17
Thursday Detuning with amplitude

• fprintf('nu_x g, nu_y g \n',nu(kx1),nu(ky1
  ))
• figure(ff2)plot(nux,nuy,'o')axis(0 0.26 0
  0.21)
• Now with even larger horizontal amplitude (10
  mm x, 100 mciron y)
• X00.0100.0001000
• X1ringpass(THERING,X0,1024)
• fftxabs(fft(X1(1,)))fftyabs(fft(X1(3,)))
• figure(ff)subplot(2,1,1)plot(nu,fftx/max(fftx),'
  m')
• subplot(2,1,2)plot(nu,ffty/max(ffty),'m')
• maxx,kxmax(fftx(2512))maxy,kymax(ffty(251
  2))
• nux(4)nu(kx1)nuy(4)nu(ky1)
• fprintf('nu_x g, nu_y g \n',nu(kx1),nu(ky1
  ))
• figure(ff2)plot(nux,nuy,'o')axis(0 0.26 0
  0.21)
• Note that the fundamental tune lines shift in
  frequency and change their
• shape (the plot normalizes their amplitude to
  one, since otherwise you
• would not see the small initial tunelines on
  the final plot). Also note that
• with increasing amplitude more and more small
  resonance lines besides
• the main tune lines show up (compare lecture on
  resonance driving
• terms).

18
Thursday Decoherence

• To show how the tuneshift with amplitude
  together with the finite beamsize
• will create decoherence of the coherent bunch
  oscillations, we track an
• ensemble of particles (100 particles for 100
  turns)
• Initial conditions (100 random coordinates with
  300 micron horizontal
• beamsize, 20 micron vertical beamsize and
  corresponding divergences
• ALS, 1.9 GeV, insertion device straights, plus
  1 mm x-offset)
• X0 300e-6randn(1,100)1.0e-3300e-6/13.5randn
  (1,100)
• 20e-6randn(1,100)1.0e-320e-6/3.65randn(1,100)
  
• zeros(1,100)zeros(1,100)
• Track for 100 turns, use 'reuse' flag to speed
  up tracking (can be used if
• lattice does not change, after every change of
  a lattice parameters, one
• has to track at least once without the 'reuse'
  flag)
• X1ringpass(THERING,X0,100,'reuse')
• Bring matrices into the right shape (we are
  only interested in x-x
• phasespace)
• xreshape(X1(1,),100,100)xpreshape(X1(2,),100,
  100)
• Plot results (upper plot shows x-x' phasespace,
  lower plot beam center)
• f1figurefigure(f1)subplot(2,1,1)plot(X0(1,),X
  0(2,),'.')
• axis(min(min(x)) max(max(x)) min(min(xp))
  max(max(xp)))

19
Thursday Decoherence

• subplot(2,1,2)plot(0,mean(X0(1,)),'.')
• axis(0 length(x) min(min(x)) max(max(x)))hold
  on
• Now plot information turn-by-turn (movie)
• for loop1100
• figure(f1)subplot(2,1,1)plot(x(,loop),xp(,loop
  ),'.')
• axis(min(min(x)) max(max(x)) min(min(xp))
  max(max(xp)))
• subplot(2,1,2)plot(loop,mean(x(,loop)),'.')
• pause(1)
• end
• At this amplitude (1 mm) one hardly sees any
  decoherence
• Now we repeat this with another 100 random
  particles with the same
• beamsize, but a 5 mm x-offset
• X0 300e-6randn(1,100)5.0e-3300e-6/13.5randn
  (1,100)
• 20e-6randn(1,100)1.0e-320e-6/3.65randn(1,100)
  
• zeros(1,100)zeros(1,100)
• X1ringpass(THERING,X0,100,'reuse')
• xreshape(X1(1,),100,100)xpreshape(X1(2,),100,
  100)
• f1figurefigure(f1)subplot(2,1,1)plot(X0(1,),X
  0(2,),'.')

20
Thursday Decoherence

• subplot(2,1,2)plot(0,mean(X0(1,)),'.')
• axis(0 length(x) min(min(x)) max(max(x)))hold
  on
• for loop1100
• figure(f1)subplot(2,1,1)plot(x(,loop),xp(,loop
  ),'.')
• axis(min(min(x)) max(max(x)) min(min(xp))
  max(max(xp)))
• subplot(2,1,2)plot(loop,mean(x(,loop)),'.')
• pause(1)
• end
• The decoherence is visible, the particles with
  small amplitudes complete
• little over ½ additional rotation in
  phasespace, compared to the ones with
• large amplitude.
• Now we repeat this with another 100 random
  particles with the same
• beamsize, but an 11 mm x-offset
• X0 300e-6randn(1,100)11.0e-3300e-6/13.5rand
  n(1,100)
• 20e-6randn(1,100)1.0e-320e-6/3.65randn(1,100)
  
• zeros(1,100)zeros(1,100)
• X1ringpass(THERING,X0,100,'reuse')
• xreshape(X1(1,),100,100)xpreshape(X1(2,),100,
  100)

21
Thursday Decoherence

• f1figurefigure(f1)subplot(2,1,1)plot(X0(1,),X
  0(2,),'.')
• axis(min(min(x)) max(max(x)) min(min(xp))
  max(max(xp)))
• subplot(2,1,2)plot(0,mean(X0(1,)),'.')
• axis(0 length(x) min(min(x)) max(max(x)))hold
  on
• for loop1100
• figure(f1)subplot(2,1,1)plot(x(,loop),xp(,loop
  ),'.')
• axis(min(min(x)) max(max(x)) min(min(xp))
  max(max(xp)))
• subplot(2,1,2)plot(loop,mean(x(,loop)),'.')
• pause(1)
• end
• The decoherence is fairly strong, the particles
  with small amplitudes complete
• more than one additional rotation in
  phasespace, compared to the ones with
• large amplitude, the oscillation amplitude of
  the beam centroid got much
• smaller. Even though there is no damping in the
  simulation, the particle seem
• to follow a spiral, which is created by the
  fact that particles with smaller
• amplitudes (in the ALS) have a higher tune.

22
Thursday Decoherence/Resonances

• Now to illustrate, how a single machine error
  can significantly change the
• nonlinear dynamicd, we put one gradient and one
  skew gradient error into
• our model machine
• QFI findcells(THERING,'FamName','QF')
• QF getcellstruct(THERING,'PolynomB',QFI,2)
• THERING setcellstruct(THERING,'PolynomB',QFI(9),
  QF(9)1.01,2)
• THERING setcellstruct(THERING,'PolynomA',QFI(9),
  QF(9)0.002,2)
• And we track again with the same particle
  distribution and 11 mm x-offset
• X0 300e-6randn(1,100)11.0e-3300e-6/13.5rand
  n(1,100)
• 20e-6randn(1,100)1.0e-320e-6/3.65randn(1,100)
  
• zeros(1,100)zeros(1,100)
• Now one first has to track once (for one turn)
  without the 'reuse' flag (we
• changed some magnet strengths
• X1ringpass(THERING,X0,1)
• X1ringpass(THERING,X0,100,'reuse')
• xreshape(X1(1,),100,100)xpreshape(X1(2,),100,
  100)

23
Thursday Decoherence/Resonances

• f1figurefigure(f1)subplot(2,1,1)plot(X0(1,),X
  0(2,),'.')
• axis(min(min(x)) max(max(x)) min(min(xp))
  max(max(xp)))
• subplot(2,1,2)plot(0,mean(X0(1,)),'.')
• axis(0 length(x) min(min(x)) max(max(x)))hold
  on
• for loop1100
• figure(f1)subplot(2,1,1)plot(x(,loop),xp(,loop
  ),'.')
• axis(min(min(x)) max(max(x)) min(min(xp))
  max(max(xp)))
• subplot(2,1,2)plot(loop,mean(x(,loop)),'.')
• pause(1)
• end
• The machine errors created a resonance which
  actually captures
• part of the particle distribution in a
  resonance island and therefore
• slows down the decoherence compared to the case
  without errors.
• Just in case you want to continue doing
  something else, remove the
• error from the lattice
• THERING setcellstruct(THERING,'PolynomB',QFI(9)
  ,QF(9),2)
• THERING setcellstruct(THERING,'PolynomA',QFI(9)
  ,0,2)

24
Computer Homework
25
Motivation for Homework Problems
Storage Ring Setup

• Dispersion ?x ?x / (?p/p)
• Chromaticity ? ? ? / (?p/p)
• Momentum compaction factor
• a (?L/L) / (?p/p)

26
Computer Homework 1

• Computer Homework 1 Dispersion Chromaticity
• Given the following middle layer commands
• x getam('BPMx') Horizontal orbit at
  BPMs(122x1)
• RF getrf Get the RF
  frequency (MHz)
• setrf(RF) Set the RF
  frequency (MHz)
• Tune gettune Tune NuX NuY
• Alpha getmcf Momentum
  compaction factor
• Write a Matlab script to
• 1. Plot the horizontal dispersion function ?x
  ?x / (?p/p)
• Note (?p/p) - (1/a)(?f/f)
• Use .1 energy change
• 2. Plot the horizontal tune vs energy change
• (compute at -.4 -.2 0 .2 .4
  energy changes)
• 3. Compute the horizontal chromaticity by
  using a numerical derivative.
• (I.e., the slope of horizontal tune vs
  energy change at ?p/p0)

27
Computer Homework 2 Tune Correction

• Given the following middle layer commands
• Tune gettune Tune NuX NuY
• qf getsp(QF') Get all the QF setpoints
• setsp(QF', qf.5) Change QF by .5 amp
• Write a Matlab script to
• 1. Record the starting tunes
• 2. Measure the 2x2 tune response matrix.
  I.e., the response
• from QF QD to horizontal and vertical
  tune. There are
• 24 magnets in each of the families,
  however, treat the entire
• quadrupole family as a single tune
  corrector. Use a .5 amp
• change in the quadrupole family to get
  the response matrix.
• 3. Change the tunes by -.1 .05 using the
  QF QD families.
• How much change in the QF QD families
  did it take?

28
Computer Homework 3

• Computer Homework 3 Orbit Correction
• Given the following middle layer commands
• x getam('BPMx') Horizontal orbit at
  BPMs(122x1)
• R getbpmresp Get the coupled BPM
  response matrix
• Rx R(1122,197) Just the horizontal
  part of the matrix
• hcm getsp('HCM') Get the horizontal
  corrector setpoints (97x1)
• setsp('HCM', hcm) Set the horizontal
  correctors
• 1. Create a random orbit error with the following
  commands
• hcm .05 randn(97,1) 97 horizontal
  correctors at the ALS
• setsp('HCM', hcm)
• 2. Write a Matlab script to correct the orbit
• a. Plot the singular values of the response
  matrix (semilogy)
• b. Plot the remaining RMS orbit error vs.
  of the singular values.
• Ie., orbit error when 1 sv is used, 12
  sv, 13 sv, 14 sv, etc.
• c. Plot the total corrector strength change
  (sum(abs(DeltaHCM))) vs.
• of the singular values used.

29
Computer Homework 4 Local Bumps

• Use the equations from Problem Set 2 and
  the optics data below to compute the kicks for a
  3-magnet bump.

?2
?1
Electron Beam
?3
s2
s1
s3
Corrector Magnets
30
Computer Homework 4 Local Bumps (cont.)
B. Put the 3 correction magnet changes in the
model and plot the horizontal orbit. The bump
should look basically closed at this point. C.
Since the dispersion function is not zero at the
corrector magnets, the correctors changed the
energy of the beam. Can you change the RF
frequency by the classic knob turning method
(ie, trial-and-error) to make the bump more
local? How much RF change did it take? D. In
matlab Energy, HCMEnergy getenergy will
return the energy change due to the HCMs.
Compute the required RF change to counter the
energy change due to the corrector magnet change.

• Matlab Help
• Get the horizontal orbit
• x getam('BPMx')
• Change HCM6 1 by .1 mrad (Physcis units
  are radians for a corrector magnet)
• hcm getsp('HCM', 6 1, 'Physics')
• setsp('HCM', hcm.0001, 6 1, 'Physics')
• Get the beta phase radian advance of the
  model
• beta1 modeltwiss('Beta', 'HCM', 6 1)
• phase1 modeltwiss('Phase', 'HCM', 6 1)
• Get Set the RF frequency
• RF getrf
• setrf(RF)