Helmut Weick, 12. Dec 2014

GICOSY was not orginally made for storage rings but if you know how you can calculate and design rings. All parameters can be derived from the system transfer matrix and the system variables.

BRHO = 13.0 ; [Tm]

PC = BRHO * CHARGE *2.99792458E8 *1E-6 ; [MeV]

M0 = MASS * 931.494061 ; [MeV]

ENERGY = SQRT(PC*PC+M0*M0) - M0 ; [MeV]

GAMMA = 1 + ENERGY/M0 ;

ENE = ENERGY/MASS ; [MeV/u]

BETA = SQRT(1-(1/GAMMA^2)) ;

;

R P ENERGY MASS CHARGE ;

The usual definition of the beam is with the Twiss parameters.

T X "α_{x}" "β_{x}" "emittance_{x}" ;

T Y "α_{y}" "β_{y}" "emittance_{y}" ;

In GICOSY the system is described by the transfer matrix of the full system starting with an upright ellipse.
The case with a non-upright ellipse (Twiss parameter α <> 0) is introduced by an object lens at the beginning.
In this case the system matrix after the start of the system will not be unity but have in addition a coefficient
[A,X] = -α_{x} / β_{x} or [B,Y]=-α_{y} / β_{y}.

The standad calculation is always for a beam in the middle of the dipoles, but one can define a size of
the momentum deviation, that is used for plots later.

D P 0 "dp/p" ;

Sometimes you also want an initial dispersion. You can define it one position by simply fixing the system
matrix manually

[X,D] = 0.55 ; [m]

[A,D] = 0.0 ; (slope of dispersion)

There are a number of system variables which are updated along the path. They all end with $.
The beam parameters are one of them

TXA$ = Twiss parameter α in x plane

TYB$ = Twiss parameter β in Y plane

...

The full list is available from the
reference manual.

Plots along the ring are created in a separate loop at the end after running the rest of the program.
The system variables are also updated in the plot loops, whereas the other variables stay constant
as assigned before.

Example: plot of Y-envelope along the ring. The emittance EX and the momentum, deviation D need to be defined before.

D E 0.1 0.070 100 (SQRT(EX*TXB$)+([X,D]*D)) ; (plot of beam envelope including dispersion)

D E 0.20 50.0 20 ([X,D]) ;
(plot of dispersion function)

For text output to file (P F) or scree (P S) the matrix coeffients and system variables have to be transferred
first to a normal variable.

DISPERSION = [X,D] ;

P S '(X,D)=',DISPERSION ;

B B "name" ; (block begin) and

B E ; (block end)

and then combine many blocks with

I B "name" ; (include block)

also in reversed order

I B "name" R ;

A while loop can be created with a counter variable

"counter" = 10 ;
W B "counter";

"counter" = "counter" -1 ;
(...)

W E ;

A first condition could be to demand an upright ellipse that can be put into fit variable

MIN = TXA$ ;

F V MIN ;

In reality the fit function will contain the quadratic sum of many requirements.

A typical task is to find the beam for which the beam is matched. In case of initial α=0 it simply becomes:

BETAX = SQRT(ABS([X,A]^2 / (1-[X,X]^2))) ;

In case you want to do a fit to find the matching beam be sure to start your calculation with α=0,
otherwise you include the object lens in your system matrix. Then you can do the multiplication of the matrix
above in a separate calculation in GICOSY and use a fit to match initial and final α and β .

TXB = TXBI$ ;

TXA = TXAI$ ;

F B TXB TXA ;

TXG = (1+TXA^2)/TXB ;

HELP = [X,A]*[A,X] + [X,X]*[A,A] ;

TXB2 = [X,X]^2 *TXB -2*[X,X]*[X,A] *TXA +[X,A]^2 *TXG ;

TXA2 =-[X,X]*[A,X] *TXB +HELP *TXA -[X,A]*[A,A] *TXG ;

MIN = (TXB-TXB2)^2*999 + (TXA-TXA2)^2*999 ;

F V MIN ;

F E 1E-10 1000 1 ;

P S 'fit: TXA=',TXA,' TXB=',TXB,' MIN=',MIN ;

A good verification is the plot of the beta function around the ring, beginning and end must match.

D E 0.1 80. 100 (TYB$) ;

PHIX = ASIN(([X,A]/SQRT(TXBI$*TXB$))) ;

The number still has ambiguity regarding multiples of π because of the asin(). It helps to look at the plot of the phase advance along the ring and to count how often (n) the curve touches π/2 or 0.

The tune (Q) then is Q = n + Δ Φ / 2π . Once you know the number it usually does not change so much in integers.

The phase advance is also available in the output of the P E command.

For a periodic symmetric lattice we can also write:

PHIX = ACOS(([X,X]+[A,A])/2) ;

Chromaticity can be calculated by the same expression but adding chromatic terms.
Just replace the coefficient by one for shifted momentum.

PHIX2 = ACOS(0.5*ABS([X,X]+[A,A]+[X,XD]*D+[A,AD]*D)) ;

The difference of phases divided by the tune gives the usual first order chromaticity.

C M 3 ;

;F T H ;

When selecting C M 3 the fring-field integral table (F T) must not be used otherwise you will still get C M 1 ; The definition of path length is a bit different in both modes. The variable L in calculation mode (CM) 2+3 means path length difference, but in mode 1 it is the time difference divided by the reference velocity.

Only for output of single element matrices (P A) the variable L really is the path length difference also in C M 1.

So for the usual definition use C M 2 or 3 even if it is slower. Time deviations are the same in both modes.

Very important are the coordinate types for output, with this the meaning of some variables will change. For example D then means relative momentum or energy deviation.

O C N E T S ; energy and time scaled (dT/T_0) best for isoc. mode

O C N M L N ; momentum and path length mode needed for matrix output to MOCADI

You can also choose between time (T) or path length deviation (L), both scaled (S) to relative deviation or absolute (N) in seconds or meters.

In the isochronous CR for example the output of [L,D] after one turn then has different meaning and values: The letter will also be "L" but sometimes it is actually time.

O C N M L N (momentum, path length dL) [L,D] = 79.77 m

O C N M L S (momentum, path length dL/L) [L,D] = 0.360 (as the circumference is 221.45 m)

O C N M T N (momentum, time dT) [L,D] = 1.39E-11 s

O C N M T S (momentum, time dT/T) [L,D] = 1.51E-5 (one turn takes 0.924E-6 s)

So in mode (O C N M L S) the coefficient [L,D] actually is alpha_p (momentum compaction factor),
whereas in mode (O C N M T S) it is [L,D] = eta (slip factor).

When choosing energy as dispersive coordinate the numbers are a bit different, here at gamma=1.67

O C N E T S (energy, path length dL/L) [L,D] = 0.225 (as delta_E > delta_p)

A simple fit of isochronicity also in higher order simply means setting [L,D^n]=0 in time mode.

For faster calculation you can still use C M 1 with GIOS style fringe field integrals,
just the path length deviation is different. As the slip factor is given correctly you can still easily
compute alpha_p.

O C N M T S ; (momentum, time dT/T)

(...)

ALPHAP = [L,D] + 1/GAMMA^2 ;

When you compare calculations in CM1 or CM3 they will differ slightly. This is mainly due to the different description
of the fringe fields. For the GICOSY standard fringe fields the Integrals in CM1 were derived for the same Enge functions
like in CM2 or CM3. But for new fringe fields one must check carefully that one really calculates the same fields.

Watch out, as it is old style FORTRAN and the line for an arithmetic expression must not be longer than 72 characters.

One full example for GICOSY for the CR @ FAIR can be found on the MOCADI web page.

(isochronous CR).

*
This page was last updated by Dr. Helmut Weick, 25 ^{th} May 2018,
contact h.weick(at)gsi.de,
Imprint (Impressum),
Privacy Policy (Datenschutzerklärung)
*