Some General Comments:
Synchrotron Tune vs. Amplitude
rm(list=ls())
E0 = 938 # MeV
eV = 0.3 # MeV
h = 12
Es = E0+200 # MeV
gamt = 8.5
Nturns = 2000
The Alternating Gradient Synchrotron (AGS) at Brookhaven National Laboratory operates at a top energy of 24 GeV (Kinetic Energy) and is used as the injector to the Relativistic Heavy Ion Collider (RHIC). Its injection energy was designed to be 200 MeV, below its transition energy. We can use the following parameters for the AGS at its design injection energy:
\(V\) = 0.3 MV; \(h\) = 12; \(E_s\) = 1.138 GeV; \(\gamma_t\) = 8.5
- For beam contained within a stationary bucket at this energy, what synchronous phase must be chosen?
# below transition: cos(phi_s) > 0
phis = 0
- Create a function which accepts a vector of \((\phi,\Delta E/E_s)_i\) and returns the vector \((\phi,\Delta E/E_s)_f\) corresponding to one turn about the synchrotron, say at the exit of the RF acceleration system, given the above parameters.
gam = Es/E0
vonc = sqrt(1-1/gam^2)
eta = 1/gamt^2 - 1/gam^2
G = 2*pi*h*eta/vonc^2
K = eV/Es
SynchTrk = function(U){
U[1] = U[1] + G*U[2]
U[2] = U[2] + K*(sin(U[1]) - sin(phis) )
return(U)
}
- Take a proton that starts with initial \(\Delta E/E_s\) = 0 and an initial phase \(\phi = \pi/20\). Make a plot in \(\phi-\Delta E/E_s\) phase space of the particle motion over 1000 turns.
phi0 = pi/20
dEonEmax = 0.001
plot(phi0, 0, type="n", xlab="phase", ylab="dE/E",
xlim=c(-1,1)*pi, ylim=c(-1,1)*dEonEmax )
# track the particle...
U = c(phi0,0)
i = 0
while(i<Nturns*8){
i=i+1
U = SynchTrk(U)
points(U[1], U[2], type="p",pch=".")
}
- R has a built-in Fast Fourier Transform (FFT) function that allows one to analyze data for frequency content. The frequency (or, in this case, the “synchrotron tune”) of the particle’s motion can be estimated by looking at what tune value corresponds to the peak of the FFT coefficients. The following is a code chunck that creates a function to perform an FFT on a set of data (
x
), and finds the value of the peak coefficient:
# Fast Fourier analyze the data to extract a "tune"
fftpeak = function(x){
zz_x = abs(fft(x))
f_x = c(1:length(x))/length(x)
f_xfft = f_x[ max.col( t( cbind( f_x, zz_x) ) )[2] ]
nu_x = min(f_xfft, (1-f_xfft))
nu_x
}
- As a trial, let’s take an oscillating set of data and find its frequency using the FFT function:
ntr = c(1:500)
xtrial = cos(2*pi*0.3812*ntr) # ie., cos(2pi * tune * n)
fftpeak(xtrial)
[1] 0.384
(Try varying the number of points used to see how close the FFT frequency comes to the original value.)
Again, track a particle for 1000 turns (Note: FFT calculations often perform best when using \(2^N\) data points) which has initial conditions \(\phi_0 = \pi/20\) and \(\Delta E_0/E_s\) = 0. This time, keep track of \(\phi_n\) each turn.
What synchrotron tune does an FFT analysis of the phase data (\(\phi_n\)) show?
phi0 = pi/20
rm(U)
U = c(phi0,0)
phi=phi0
i = 0
while(i<Nturns){
i=i+1
U = SynchTrk(U)
phi[i] = U[1]
}
# Fast Fourier analyze the data to extract a "tune"
nu_u = fftpeak(phi)
nu_u
[1] 0.032
- How does this compare with the analytically predicted value?
nus_theory = sqrt(-h*eta*eV*cos(phis)/(2*pi*vonc^2*Es))
nus_theory
[1] 0.03232919
Na = 250*2^c(0:5)
- Try varying the number of turns tracked, such as 250, 500, 1000, 2000, 4000, 8000, etc. For what number of turns is the FFT result within 1% of the predicted result?
# track the particle Na turns...
nu_s = 0
i = 0
while(i<length(Na)){
i = i+1
phi0 = pi/20
U = c(phi0,0)
phi=phi0
j=0
while(j<Na[i]){
j=j+1
U = SynchTrk(U)
phi[j] = U[1]
}
nu_s[i] = fftpeak(phi)
}
plot(Na,nu_s,ylim=c(0.025,.035))
abline(h=nus_theory*c(0.99,1,1.01),lty=c(2,1,2))
- Next, to get an accurate result, track a particle through \(1.6384\times 10^{4}\) turns. Track for various initial \(\phi_0\)’s from \(5^\circ\) toward \(180^\circ\) and keep a tally of the synchrotron tune.
- How does the tune behave? Make a plot of \(\nu_s\) vs. \(\phi_0\) for the range within the stable region (separatrix).
Nturns = 16000
phi_init = c(1:36)*5*pi/180
nu_s = 0
i = 0
while(i<length(phi_init)){
i = i+1
phi0 = phi_init[i]
U = c(phi0,0)
phi=phi0
j=0
while(j<Nturns){
j=j+1
U = SynchTrk(U)
phi[j] = U[1]
}
nu_s[i] = fftpeak(phi)
}
plot(phi_init/pi*180,nu_s,ylim=c(0,nus_theory),
xlab="Phase Amplitude (deg.)",ylab="Synchrotron Tune")
abline(h=c(0,nus_theory),lty=c(1,2),col=c("black","red"))
text(150,0.030,"Small Ampl. Limit",col="red")
- What can you say about the particle motion on or very near the separatrix?
# The motion slows down; oscillations would cease *on* the separatrix.
Island Hopping
nu = 0.4208090190
Nturns = 1000
Let’s model the trajectory of a particle in an otherwise ideal storage ring but with a single sextupole present. The ideal betatron tune is set to \(\nu\) = 0.420809. As shown in class, the mapping of phase space variables for a single turn around the ring can be expressed with the following function:
cn = cos(2*pi*nu)
sn = sin(2*pi*nu)
SextTrack = function(U){
ac = U[1]
as = U[2]
U[1] = ac*cn + (as-ac^2)*sn
U[2] = -ac*sn + (as-ac^2)*cn
return(U)
}
- Augment the above code so as to track a single particle for 1000 turns and plot the trajectory in \(x,p\) phase space (where \(p\equiv \beta x+\alpha x'\)). Try the following initial coordinates:
x0 = 0
p0 = c(0.05,0.1,0.2,0.5,1.0,1.4)
Hint: For the sake of comparison, make a series of plots, with equal axes (use asp=1
) and arrange in a grid using the par(mfrow=c(3,2))
function in R prior to plotting, where c(3,2)
indicates there are three rows of two plots each. Naturally, you may wish to select different numbers for these. (If you use ggplot2, etc., there are other ways of setting up a grid; but par
works in base R.)
p0 = c(0.05,0.1,0.2,0.5,1.0,1.4)
par(mfrow=c(1,2))
i=0
while(i<length(p0)){
i=i+1
U = c(0,p0[i])
x=0
p=0
j=0
while(j<Nturns){
j=j+1
U = SextTrack(U)
x[j] = U[1]
p[j] = U[2]
}
plot(x,p,pch=".",ylim=c(-2,2),asp=1)
}
- We now wish to simulate an experiment, where we take a beam of particles and “kick” the beam with a kicker magnet to explore the phase space:
Create a distribution of 2000 particles each with coordinates \(x\) and \(p\), with mean(x) = mean(p) = 0
, and sd(x) = sd(p) = 0.05
.
Give all the particles of the entire distribution a “kick” in the \(p\) coordinate by same amount, say, kck = 0.05
.
Now track the entire distribution for 1000 turns, keeping tally of xave[n] = mean(x)
for each of the n
turns. We can imagine that the data xave
are results of average beam position measurements at a monitor near the sextupole, taken once per revolution.
Make plots of xave
vs. n
for the following conditions:
1 |
0.05 |
0.05 |
2 |
0.05 |
0.10 |
3 |
0.05 |
0.20 |
4 |
0.10 |
0.20 |
5 |
0.10 |
0.50 |
6 |
0.10 |
1.00 |
7 |
0.20 |
0.50 |
8 |
0.20 |
1.00 |
9 |
0.20 |
1.40 |
sigx = c(0.05, 0.05, 0.05, 0.10, 0.10, 0.10, 0.20, 0.20, 0.20)
kck = c(0.05, 0.10, 0.20, 0.20, 0.50, 1.00, 0.50, 1.00, 1.40)
Npar = 2000
par(mfrow=c(2,2))
m = 0
while(m<length(sigx)){
m=m+1
x0 = rnorm(Npar,0,sigx[m])
p0 = rnorm(Npar,0,sigx[m])
xf=x0
pf=p0 + kck[m]
xave=0
i=0
while(i<Nturns){
i=i+1
j=0
while(j<Npar){
j=j+1
U = c(xf[j],pf[j])
U = SextTrack(U)
xf[j] = ifelse(U[1]<10,U[1],NA)
pf[j] = ifelse(U[2]<10,U[2],NA)
}
xave[i] = mean(xf,na.rm=TRUE)
}
plot(xave,typ="l")
}
- For condition (5) above, why does the signal die away? For condition (9), why does the signal persist?
- For condition (9), take the vector of data containing the mean as a function of turn number (
xave
above) and select only the last several hundred turns. (Note: This can be done in R by the following command: xend = xave[200:1000]
, for example.) Perform an FFT on this subset of data and comment on the answer. How does the frequency measured coincide with the phase space for this set of conditions?
xend = xave[200:1000]
fftpeak(xend)
[1] 0.3982522
We see here the “island tune”; particles in phase space near the center of an island, will “hop” from one island to another, skipping every other island, until it returns to the original island after 5 turns – the tune is thus 2/5 = 0.4.
As a check, look at the following:
x= 0.35
p=-1.5
j=0
U=c(x,p)
while(j<50){
j=j+1
U = SextTrack(U)
x[j] = U[1]
p[j] = U[2]
}
plot(x,p,typ="l",xlim=c(-2,2),ylim=c(-2,2),asp=1)
LS0tCnRpdGxlOiAiSG9tZXdvcmsgNyBEaXNjdXNzaW9uIgphdXRob3I6IE1pa2UgU3lwaGVyc15bTm9ydGhlcm4gSWxsaW5vaXMgVW5pdmVyc2l0eSBhbmQgRmVybWkgTmF0aW9uYWwgQWNjZWxlcmF0b3IKICBMYWJvcmF0b3J5XQpkYXRlOiAnbGFzdCB1cGRhdGU6IGByIGZvcm1hdCgoU3lzLkRhdGUoKSksICIlZCAlYiAlWSIpYCcKb3V0cHV0OiAgaHRtbF9ub3RlYm9vawotLS0KCioqU29tZSBHZW5lcmFsIENvbW1lbnRzOioqCgotIGdyYWRlOiAgMjAgcG9zc2libGUgIAogICAgICArIDEwIHB0cyBlYWNoIHByb2JsZW0gIAogICAgICArIGF2ZSBzY29yZTogIDE3LzIwICAKCi0ga2ljayB3YXMgIm1lYW50IiB0byBiZSBvbmx5IG9uY2UgLS0gdG8gc2V0IHRoZSBpbml0aWFsIGNvbmRpdGlvbnMgIAotIHVzZSAiYXNwPTEiIHdoZW4gd2FudCB0aGUgKnNjYWxlKiBvbiB0aGUgYXhlcyB0byBoYXZlIGVxdWFsIGxlbmd0aHMgIAoKCgojIyAqKlN5bmNocm90cm9uIFR1bmUgdnMuIEFtcGxpdHVkZSoqCgpgYGB7ciwgZWNobz1UUlVFfQpybShsaXN0PWxzKCkpCgpFMCA9ICA5MzggICMgTWVWCmVWID0gIDAuMyAgIyBNZVYKaCAgPSAgMTIKRXMgPSBFMCsyMDAgIyBNZVYKZ2FtdCA9IDguNQpOdHVybnMgPSAyMDAwCmBgYAoKVGhlIEFsdGVybmF0aW5nIEdyYWRpZW50IFN5bmNocm90cm9uIChBR1MpIGF0IEJyb29raGF2ZW4gTmF0aW9uYWwgTGFib3JhdG9yeSBvcGVyYXRlcyBhdCBhIHRvcCBlbmVyZ3kgb2YgMjQgR2VWIChLaW5ldGljIEVuZXJneSkgYW5kIGlzIHVzZWQgYXMgdGhlIGluamVjdG9yIHRvIHRoZSBSZWxhdGl2aXN0aWMgSGVhdnkgSW9uIENvbGxpZGVyIChSSElDKS4gIEl0cyBpbmplY3Rpb24gZW5lcmd5IHdhcyBkZXNpZ25lZCB0byBiZSAyMDAgTWVWLCBiZWxvdyBpdHMgdHJhbnNpdGlvbiBlbmVyZ3kuIFdlIGNhbiB1c2UgdGhlIGZvbGxvd2luZyBwYXJhbWV0ZXJzIGZvciB0aGUgQUdTIGF0IGl0cyBkZXNpZ24gaW5qZWN0aW9uIGVuZXJneToKCj4gJFYkID0gYHIgZVZgIE1WOyAkaCQgPSBgciBoYDsgJEVfcyQgPSBgciBFcy8xMDAwYCBHZVY7ICRcZ2FtbWFfdCQgPSBgciBnYW10YAoKCmEuIEZvciBiZWFtIGNvbnRhaW5lZCB3aXRoaW4gYSAqc3RhdGlvbmFyeSBidWNrZXQqIGF0IHRoaXMgZW5lcmd5LCB3aGF0IHN5bmNocm9ub3VzIHBoYXNlIG11c3QgYmUgY2hvc2VuPwoKYGBge3IsIGVjaG89VFJVRX0KIyBiZWxvdyB0cmFuc2l0aW9uOiAgY29zKHBoaV9zKSA+IDAKcGhpcyA9IDAKYGBgCgpiLiBDcmVhdGUgYSBmdW5jdGlvbiB3aGljaCBhY2NlcHRzIGEgdmVjdG9yIG9mICQoXHBoaSxcRGVsdGEgRS9FX3MpX2kkIGFuZCByZXR1cm5zIHRoZSB2ZWN0b3IgJChccGhpLFxEZWx0YSBFL0VfcylfZiQgY29ycmVzcG9uZGluZyB0byBvbmUgdHVybiBhYm91dCB0aGUgc3luY2hyb3Ryb24sIHNheSBhdCB0aGUgZXhpdCBvZiB0aGUgUkYgYWNjZWxlcmF0aW9uIHN5c3RlbSwgZ2l2ZW4gdGhlIGFib3ZlIHBhcmFtZXRlcnMuCgpgYGB7ciwgZWNobz1UUlVFfQpnYW0gID0gRXMvRTAKdm9uYyA9IHNxcnQoMS0xL2dhbV4yKQpldGEgID0gMS9nYW10XjIgLSAxL2dhbV4yCkcgPSAyKnBpKmgqZXRhL3ZvbmNeMgpLID0gZVYvRXMKU3luY2hUcmsgPSBmdW5jdGlvbihVKXsKICAgVVsxXSA9IFVbMV0gKyBHKlVbMl0KICAgVVsyXSA9IFVbMl0gKyBLKihzaW4oVVsxXSkgLSBzaW4ocGhpcykgKQogICByZXR1cm4oVSkKfQpgYGAKCmMuIFRha2UgYSBwcm90b24gdGhhdCBzdGFydHMgd2l0aCBpbml0aWFsICRcRGVsdGEgRS9FX3MkID0gMCBhbmQgYW4gaW5pdGlhbCBwaGFzZSAkXHBoaSA9IFxwaS8yMCQuICBNYWtlIGEgcGxvdCBpbiAkXHBoaS1cRGVsdGEgRS9FX3MkIHBoYXNlIHNwYWNlIG9mIHRoZSBwYXJ0aWNsZSBtb3Rpb24gb3ZlciBgciBOdHVybnNgIHR1cm5zLgoKYGBge3IgZWNobz1UUlVFLGV2YWw9VFJVRX0KcGhpMCA9IHBpLzIwCgpkRW9uRW1heCA9IDAuMDAxCnBsb3QocGhpMCwgMCwgdHlwZT0ibiIsIHhsYWI9InBoYXNlIiwgeWxhYj0iZEUvRSIsCiAgIHhsaW09YygtMSwxKSpwaSwgeWxpbT1jKC0xLDEpKmRFb25FbWF4ICkKCiMgdHJhY2sgdGhlIHBhcnRpY2xlLi4uClUgPSBjKHBoaTAsMCkKaSA9IDAKd2hpbGUoaTxOdHVybnMqOCl7CiAgIGk9aSsxCiAgIFUgPSBTeW5jaFRyayhVKQogIAlwb2ludHMoVVsxXSwgVVsyXSwgdHlwZT0icCIscGNoPSIuIikJCn0KYGBgCgpkLiAqKlIqKiBoYXMgYSBidWlsdC1pbiBGYXN0IEZvdXJpZXIgVHJhbnNmb3JtIChGRlQpIGZ1bmN0aW9uIHRoYXQgYWxsb3dzIG9uZSB0byBhbmFseXplIGRhdGEgZm9yIGZyZXF1ZW5jeSBjb250ZW50LiBUaGUgZnJlcXVlbmN5IChvciwgaW4gdGhpcyBjYXNlLCB0aGUgInN5bmNocm90cm9uIHR1bmUiKSBvZiB0aGUgcGFydGljbGUncyBtb3Rpb24gY2FuIGJlIGVzdGltYXRlZCBieSBsb29raW5nIGF0IHdoYXQgdHVuZSB2YWx1ZSBjb3JyZXNwb25kcyB0byB0aGUgcGVhayBvZiB0aGUgRkZUIGNvZWZmaWNpZW50cy4gIFRoZSBmb2xsb3dpbmcgaXMgYSBjb2RlIGNodW5jayB0aGF0IGNyZWF0ZXMgYSBmdW5jdGlvbiB0byBwZXJmb3JtIGFuIEZGVCBvbiBhIHNldCBvZiBkYXRhIChgeGApLCBhbmQgZmluZHMgdGhlIHZhbHVlIG9mIHRoZSBwZWFrIGNvZWZmaWNpZW50OgoKYGBge3IgRkZUZmNuLCBlY2hvPVRSVUV9CiMgRmFzdCBGb3VyaWVyIGFuYWx5emUgdGhlIGRhdGEgdG8gZXh0cmFjdCBhICJ0dW5lIgpmZnRwZWFrID0gZnVuY3Rpb24oeCl7CiAgIHp6X3ggPSAgYWJzKGZmdCh4KSkKICAgZl94ICA9IGMoMTpsZW5ndGgoeCkpL2xlbmd0aCh4KQogICBmX3hmZnQgPSBmX3hbIG1heC5jb2woIHQoIGNiaW5kKCBmX3gsIHp6X3gpICkgKVsyXSBdCiAgIG51X3ggPSBtaW4oZl94ZmZ0LCAoMS1mX3hmZnQpKQogICBudV94Cn0KYGBgCgogICAtIEFzIGEgdHJpYWwsIGxldCdzIHRha2UgYW4gb3NjaWxsYXRpbmcgc2V0IG9mIGRhdGEgYW5kIGZpbmQgaXRzIGZyZXF1ZW5jeSB1c2luZyB0aGUgRkZUIGZ1bmN0aW9uOiAgCgpgYGB7cn0KbnRyID0gYygxOjUwMCkKeHRyaWFsID0gY29zKDIqcGkqMC4zODEyKm50cikgICMgIGllLiwgY29zKDJwaSAqIHR1bmUgKiBuKQpmZnRwZWFrKHh0cmlhbCkKYGBgCihUcnkgdmFyeWluZyB0aGUgbnVtYmVyIG9mIHBvaW50cyB1c2VkIHRvIHNlZSBob3cgY2xvc2UgdGhlIEZGVCBmcmVxdWVuY3kgY29tZXMgdG8gdGhlIG9yaWdpbmFsIHZhbHVlLikgIAoKCiAgIC0gQWdhaW4sIHRyYWNrIGEgcGFydGljbGUgZm9yIGByIE50dXJuc2AgdHVybnMgKE5vdGU6IEZGVCBjYWxjdWxhdGlvbnMgb2Z0ZW4gcGVyZm9ybSBiZXN0IHdoZW4gdXNpbmcgJDJeTiQgZGF0YSBwb2ludHMpIHdoaWNoIGhhcyBpbml0aWFsIGNvbmRpdGlvbnMgJFxwaGlfMCA9IFxwaS8yMCQgYW5kICRcRGVsdGEgRV8wL0VfcyQgPSAwLiAgVGhpcyB0aW1lLCBrZWVwIHRyYWNrIG9mICRccGhpX24kIGVhY2ggdHVybi4gICAKCiAgIC0gV2hhdCBzeW5jaHJvdHJvbiB0dW5lIGRvZXMgYW4gRkZUIGFuYWx5c2lzIG9mIHRoZSBwaGFzZSBkYXRhICgkXHBoaV9uJCkgc2hvdz8gIAoKYGBge3IsIGVjaG89VFJVRSwgZXZhbD1UUlVFfQpwaGkwID0gcGkvMjAKcm0oVSkKVSA9IGMocGhpMCwwKQpwaGk9cGhpMAppID0gMAp3aGlsZShpPE50dXJucyl7CiAgIGk9aSsxCiAgIFUgPSBTeW5jaFRyayhVKQogIAlwaGlbaV0gPSBVWzFdCQp9CgojIEZhc3QgRm91cmllciBhbmFseXplIHRoZSBkYXRhIHRvIGV4dHJhY3QgYSAidHVuZSIKbnVfdSA9IGZmdHBlYWsocGhpKQpudV91CmBgYAoKICAgLSBIb3cgZG9lcyB0aGlzIGNvbXBhcmUgd2l0aCB0aGUgYW5hbHl0aWNhbGx5IHByZWRpY3RlZCB2YWx1ZT8gIAoKYGBge3IsIGVjaG89VFJVRSwgZXZhbD1UUlVFfQpudXNfdGhlb3J5ID0gc3FydCgtaCpldGEqZVYqY29zKHBoaXMpLygyKnBpKnZvbmNeMipFcykpCm51c190aGVvcnkKYGBgCgpgYGB7ciwgZWNobz1UUlVFfQpOYSA9IDI1MCoyXmMoMDo1KQpgYGAKCiAgIC0gVHJ5IHZhcnlpbmcgdGhlIG51bWJlciBvZiB0dXJucyB0cmFja2VkLCBzdWNoIGFzIGByIE5hYCwgKmV0Yy4qICBGb3Igd2hhdCBudW1iZXIgb2YgdHVybnMgaXMgdGhlIEZGVCByZXN1bHQgd2l0aGluIDElIG9mIHRoZSBwcmVkaWN0ZWQgcmVzdWx0PyAgCgpgYGB7ciwgZWNobz1UUlVFLCBldmFsPVRSVUUsIHJlc3VsdHM9J21hcmt1cCd9CgojIHRyYWNrIHRoZSBwYXJ0aWNsZSBOYSB0dXJucy4uLgpudV9zID0gMAppID0gMAp3aGlsZShpPGxlbmd0aChOYSkpewogICBpID0gaSsxCiAgIHBoaTAgPSBwaS8yMAogICBVID0gYyhwaGkwLDApCiAgIHBoaT1waGkwCiAgIGo9MAogICB3aGlsZShqPE5hW2ldKXsKICAgICAgaj1qKzEKICAgICAgVSA9IFN5bmNoVHJrKFUpCiAgICAgCXBoaVtqXSA9IFVbMV0KICAgfQogIG51X3NbaV0gPSBmZnRwZWFrKHBoaSkKfQoKcGxvdChOYSxudV9zLHlsaW09YygwLjAyNSwuMDM1KSkKYWJsaW5lKGg9bnVzX3RoZW9yeSpjKDAuOTksMSwxLjAxKSxsdHk9YygyLDEsMikpCmBgYAoKICAgZS4gTmV4dCwgdG8gZ2V0IGFuIGFjY3VyYXRlIHJlc3VsdCwgdHJhY2sgYSBwYXJ0aWNsZSB0aHJvdWdoICRgciAyXjE0YCQgdHVybnMuIFRyYWNrIGZvciB2YXJpb3VzIGluaXRpYWwgJFxwaGlfMCQncyBmcm9tICQ1XlxjaXJjJCB0b3dhcmQgJDE4MF5cY2lyYyQgYW5kIGtlZXAgYSB0YWxseSBvZiB0aGUgc3luY2hyb3Ryb24gdHVuZS4gIAogICAKICAgLSBIb3cgZG9lcyB0aGUgdHVuZSBiZWhhdmU/ICBNYWtlIGEgcGxvdCBvZiAkXG51X3MkICBfdnMuXyAkXHBoaV8wJCAgZm9yIHRoZSByYW5nZSB3aXRoaW4gdGhlIHN0YWJsZSByZWdpb24gKCpzZXBhcmF0cml4KikuICAKCmBgYHtyLCBlY2hvPVRSVUUsIGV2YWw9VFJVRSwgcmVzdWx0cz0nbWFya3VwJ30KTnR1cm5zID0gMTYwMDAKcGhpX2luaXQgPSBjKDE6MzYpKjUqcGkvMTgwCgpudV9zID0gMAppID0gMAp3aGlsZShpPGxlbmd0aChwaGlfaW5pdCkpewogICBpID0gaSsxCiAgIHBoaTAgPSBwaGlfaW5pdFtpXQogICBVID0gYyhwaGkwLDApCiAgIHBoaT1waGkwCiAgIGo9MAogICB3aGlsZShqPE50dXJucyl7CiAgICAgIGo9aisxCiAgICAgIFUgPSBTeW5jaFRyayhVKQogICAgIAlwaGlbal0gPSBVWzFdCiAgIH0KICBudV9zW2ldID0gZmZ0cGVhayhwaGkpCn0KCnBsb3QocGhpX2luaXQvcGkqMTgwLG51X3MseWxpbT1jKDAsbnVzX3RoZW9yeSksCiAgIHhsYWI9IlBoYXNlIEFtcGxpdHVkZSAoZGVnLikiLHlsYWI9IlN5bmNocm90cm9uIFR1bmUiKQphYmxpbmUoaD1jKDAsbnVzX3RoZW9yeSksbHR5PWMoMSwyKSxjb2w9YygiYmxhY2siLCJyZWQiKSkKdGV4dCgxNTAsMC4wMzAsIlNtYWxsIEFtcGwuIExpbWl0Iixjb2w9InJlZCIpCmBgYAoKICAgLSBXaGF0IGNhbiB5b3Ugc2F5IGFib3V0IHRoZSBwYXJ0aWNsZSBtb3Rpb24gb24gb3IgdmVyeSBuZWFyIHRoZSBzZXBhcmF0cml4PwoKYGBge3IsIGVjaG89VFJVRSwgZXZhbD1UUlVFfQojICBUaGUgbW90aW9uIHNsb3dzIGRvd247IG9zY2lsbGF0aW9ucyB3b3VsZCBjZWFzZSAqb24qIHRoZSBzZXBhcmF0cml4LgpgYGAKCgojIyAqKklzbGFuZCBIb3BwaW5nKioKCmBgYHtyLCBlY2hvPVRSVUUsIGV2YWw9VFJVRX0KbnUgPSAwLjQyMDgwOTAxOTAKTnR1cm5zID0gMTAwMApgYGAKCkxldCdzIG1vZGVsIHRoZSB0cmFqZWN0b3J5IG9mIGEgcGFydGljbGUgaW4gYW4gb3RoZXJ3aXNlIGlkZWFsIHN0b3JhZ2UgcmluZyBidXQgd2l0aCBhIHNpbmdsZSBzZXh0dXBvbGUgcHJlc2VudC4gIFRoZSBpZGVhbCBiZXRhdHJvbiB0dW5lIGlzIHNldCB0byAkXG51JCA9IGByIG51YC4gIEFzIHNob3duIGluIGNsYXNzLCB0aGUgbWFwcGluZyBvZiBwaGFzZSBzcGFjZSB2YXJpYWJsZXMgZm9yIGEgc2luZ2xlIHR1cm4gYXJvdW5kIHRoZSByaW5nIGNhbiBiZSBleHByZXNzZWQgd2l0aCB0aGUgZm9sbG93aW5nIGZ1bmN0aW9uOgoKYGBge3J9CmNuID0gY29zKDIqcGkqbnUpCnNuID0gc2luKDIqcGkqbnUpClNleHRUcmFjayA9IGZ1bmN0aW9uKFUpewogICBhYyA9IFVbMV0KICAgYXMgPSBVWzJdCiAgIFVbMV0gPSAgYWMqY24gKyAoYXMtYWNeMikqc24KICAgVVsyXSA9IC1hYypzbiArIChhcy1hY14yKSpjbgogICByZXR1cm4oVSkKfQpgYGAKCmEuIEF1Z21lbnQgdGhlIGFib3ZlIGNvZGUgc28gYXMgdG8gdHJhY2sgYSBzaW5nbGUgcGFydGljbGUgZm9yIGByIE50dXJuc2AgdHVybnMgYW5kIHBsb3QgdGhlIHRyYWplY3RvcnkgaW4gJHgscCQgcGhhc2Ugc3BhY2UgKHdoZXJlICRwXGVxdWl2IFxiZXRhIHgrXGFscGhhIHgnJCkuICBUcnkgdGhlIGZvbGxvd2luZyBpbml0aWFsIGNvb3JkaW5hdGVzOgpgYGB7cn0KeDAgPSAwCnAwID0gYygwLjA1LDAuMSwwLjIsMC41LDEuMCwxLjQpCmBgYAoKSGludDogIEZvciB0aGUgc2FrZSBvZiBjb21wYXJpc29uLCBtYWtlIGEgc2VyaWVzIG9mIHBsb3RzLCB3aXRoIGVxdWFsIGF4ZXMgKHVzZSBgYXNwPTFgKSBhbmQgYXJyYW5nZSBpbiBhIGdyaWQgdXNpbmcgdGhlIGBwYXIobWZyb3c9YygzLDIpKWAgZnVuY3Rpb24gaW4gKipSKiogcHJpb3IgdG8gcGxvdHRpbmcsIHdoZXJlIGBjKDMsMilgIGluZGljYXRlcyB0aGVyZSBhcmUgdGhyZWUgcm93cyBvZiB0d28gcGxvdHMgZWFjaC4gIE5hdHVyYWxseSwgeW91IG1heSB3aXNoIHRvIHNlbGVjdCBkaWZmZXJlbnQgbnVtYmVycyBmb3IgdGhlc2UuICAoSWYgeW91IHVzZSBnZ3Bsb3QyLCBldGMuLCB0aGVyZSBhcmUgb3RoZXIgd2F5cyBvZiBzZXR0aW5nIHVwIGEgZ3JpZDsgYnV0IGBwYXJgIHdvcmtzIGluIGJhc2UgKipSKiouKQoKYGBge3IsIGVjaG89VFJVRSwgZXZhbD1UUlVFfQpwMCA9IGMoMC4wNSwwLjEsMC4yLDAuNSwxLjAsMS40KQpwYXIobWZyb3c9YygxLDIpKQppPTAKd2hpbGUoaTxsZW5ndGgocDApKXsKICAgaT1pKzEKICAgVSA9IGMoMCxwMFtpXSkKICAgeD0wCiAgIHA9MAogICBqPTAKICAgd2hpbGUoajxOdHVybnMpewogICAgICBqPWorMQogICAgICBVID0gU2V4dFRyYWNrKFUpCiAgICAgCXhbal0gPSBVWzFdCiAgICAgCXBbal0gPSBVWzJdCiAgIH0KcGxvdCh4LHAscGNoPSIuIix5bGltPWMoLTIsMiksYXNwPTEpCn0KYGBgCgpiLiAgV2Ugbm93IHdpc2ggdG8gc2ltdWxhdGUgYW4gZXhwZXJpbWVudCwgd2hlcmUgd2UgdGFrZSBhIGJlYW0gb2YgcGFydGljbGVzIGFuZCAia2ljayIgdGhlIGJlYW0gd2l0aCBhIGtpY2tlciBtYWduZXQgdG8gZXhwbG9yZSB0aGUgcGhhc2Ugc3BhY2U6ICAKCi0gQ3JlYXRlIGEgZGlzdHJpYnV0aW9uIG9mIDIwMDAgcGFydGljbGVzIGVhY2ggd2l0aCBjb29yZGluYXRlcyAkeCQgYW5kICRwJCwgd2l0aCBgbWVhbih4KSA9IG1lYW4ocCkgPSAwYCwgYW5kIGBzZCh4KSA9IHNkKHApID0gMC4wNWAuCgotIEdpdmUgYWxsIHRoZSBwYXJ0aWNsZXMgb2YgdGhlIGVudGlyZSBkaXN0cmlidXRpb24gYSAia2ljayIgaW4gdGhlICRwJCBjb29yZGluYXRlIGJ5IHNhbWUgYW1vdW50LCBzYXksIGBrY2sgPSAwLjA1YC4KCi0gTm93IHRyYWNrIHRoZSBlbnRpcmUgZGlzdHJpYnV0aW9uIGZvciBgciBOdHVybnNgIHR1cm5zLCBrZWVwaW5nIHRhbGx5IG9mIGB4YXZlW25dID0gbWVhbih4KWAgZm9yIGVhY2ggb2YgdGhlIGBuYCB0dXJucy4gV2UgY2FuIGltYWdpbmUgdGhhdCB0aGUgZGF0YSBgeGF2ZWAgYXJlIHJlc3VsdHMgb2YgYXZlcmFnZSBiZWFtIHBvc2l0aW9uIG1lYXN1cmVtZW50cyBhdCBhIG1vbml0b3IgbmVhciB0aGUgc2V4dHVwb2xlLCB0YWtlbiBvbmNlIHBlciByZXZvbHV0aW9uLgoKTWFrZSBwbG90cyBvZiBgeGF2ZWAgKnZzLiogYG5gIGZvciB0aGUgZm9sbG93aW5nIGNvbmRpdGlvbnM6Cgp8IHBsb3QgfCBgc2QoeCk9c2QocClgfCBga2NrYCB8IAp8LS0tLS0tfC0tLS0tLS18LS0tLS0tLS0tLS0tLS18CnwgICAgMSB8IDAuMDUgfCAwLjA1IHwgCnwgICAgMiB8IDAuMDUgfCAwLjEwIHwgCnwgICAgMyB8IDAuMDUgfCAwLjIwIHwgCnwgICAgNCB8IDAuMTAgfCAwLjIwIHwgCnwgICAgNSB8IDAuMTAgfCAwLjUwIHwgCnwgICAgNiB8IDAuMTAgfCAxLjAwIHwgCnwgICAgNyB8IDAuMjAgfCAwLjUwIHwgCnwgICAgOCB8IDAuMjAgfCAxLjAwIHwgCnwgICAgOSB8IDAuMjAgfCAxLjQwIHwgCgoKYGBge3IsIGVjaG89VFJVRSwgZXZhbD1UUlVFfQpzaWd4ID0gYygwLjA1LCAwLjA1LCAwLjA1LCAwLjEwLCAwLjEwLCAwLjEwLCAwLjIwLCAwLjIwLCAwLjIwKQprY2sgID0gYygwLjA1LCAwLjEwLCAwLjIwLCAwLjIwLCAwLjUwLCAxLjAwLCAwLjUwLCAxLjAwLCAxLjQwKQpOcGFyID0gMjAwMApgYGAKCmBgYHtyLCBlY2hvPVRSVUUsIGV2YWw9VFJVRX0KcGFyKG1mcm93PWMoMiwyKSkKbSA9IDAKd2hpbGUobTxsZW5ndGgoc2lneCkpewogICBtPW0rMQogICB4MCA9IHJub3JtKE5wYXIsMCxzaWd4W21dKQogICBwMCA9IHJub3JtKE5wYXIsMCxzaWd4W21dKQogICB4Zj14MAogICBwZj1wMCArIGtja1ttXQogICB4YXZlPTAKICAgaT0wCiAgIHdoaWxlKGk8TnR1cm5zKXsKICAgICAgaT1pKzEKICAgICAgaj0wCiAgICAgIHdoaWxlKGo8TnBhcil7CiAgICAgICAgIGo9aisxCiAgICAgICAgIFUgPSBjKHhmW2pdLHBmW2pdKQogICAgICAgICBVID0gU2V4dFRyYWNrKFUpCiAgICAgICAgIHhmW2pdID0gaWZlbHNlKFVbMV08MTAsVVsxXSxOQSkKICAgICAgICAgcGZbal0gPSBpZmVsc2UoVVsyXTwxMCxVWzJdLE5BKQogICAgICB9CiAgICAgIHhhdmVbaV0gPSBtZWFuKHhmLG5hLnJtPVRSVUUpCiAgIH0KICAgcGxvdCh4YXZlLHR5cD0ibCIpCn0KYGBgCgotIEZvciBjb25kaXRpb24gKDUpIGFib3ZlLCB3aHkgZG9lcyB0aGUgc2lnbmFsIGRpZSBhd2F5PyAgRm9yIGNvbmRpdGlvbiAoOSksIHdoeSBkb2VzIHRoZSBzaWduYWwgcGVyc2lzdD8KCmMuICBGb3IgY29uZGl0aW9uICg5KSwgdGFrZSB0aGUgdmVjdG9yIG9mIGRhdGEgY29udGFpbmluZyB0aGUgbWVhbiBhcyBhIGZ1bmN0aW9uIG9mIHR1cm4gbnVtYmVyIChgeGF2ZWAgYWJvdmUpIGFuZCBzZWxlY3Qgb25seSB0aGUgbGFzdCBzZXZlcmFsIGh1bmRyZWQgdHVybnMuICAoTm90ZTogIFRoaXMgY2FuIGJlIGRvbmUgaW4gKipSKiogYnkgdGhlIGZvbGxvd2luZyBjb21tYW5kOiAgYHhlbmQgPSB4YXZlWzIwMDoxMDAwXWAsIGZvciBleGFtcGxlLikgIFBlcmZvcm0gYW4gRkZUIG9uIHRoaXMgc3Vic2V0IG9mIGRhdGEgYW5kIGNvbW1lbnQgb24gdGhlIGFuc3dlci4gIEhvdyBkb2VzIHRoZSBmcmVxdWVuY3kgbWVhc3VyZWQgY29pbmNpZGUgd2l0aCB0aGUgcGhhc2Ugc3BhY2UgZm9yIHRoaXMgc2V0IG9mIGNvbmRpdGlvbnM/CgpgYGB7ciwgZWNobz1UUlVFLCBldmFsPVRSVUV9CnhlbmQgPSB4YXZlWzIwMDoxMDAwXQpmZnRwZWFrKHhlbmQpCmBgYAoKV2Ugc2VlIGhlcmUgdGhlICJpc2xhbmQgdHVuZSI7IHBhcnRpY2xlcyBpbiBwaGFzZSBzcGFjZSBuZWFyIHRoZSBjZW50ZXIgb2YgYW4gaXNsYW5kLCB3aWxsICJob3AiIGZyb20gb25lIGlzbGFuZCB0byBhbm90aGVyLCBza2lwcGluZyBldmVyeSBvdGhlciBpc2xhbmQsIHVudGlsIGl0IHJldHVybnMgdG8gdGhlIG9yaWdpbmFsIGlzbGFuZCBhZnRlciA1IHR1cm5zIC0tIHRoZSB0dW5lIGlzIHRodXMgMi81ID0gMC40LgoKQXMgYSBjaGVjaywgbG9vayBhdCB0aGUgZm9sbG93aW5nOgpgYGB7cn0KeD0gMC4zNQpwPS0xLjUKaj0wClU9Yyh4LHApCndoaWxlKGo8NTApewogICBqPWorMQogICBVID0gU2V4dFRyYWNrKFUpCiAgCXhbal0gPSBVWzFdCiAgCXBbal0gPSBVWzJdCn0KcGxvdCh4LHAsdHlwPSJsIix4bGltPWMoLTIsMikseWxpbT1jKC0yLDIpLGFzcD0xKQpgYGAKCg==