diff --git a/tracy/tracy/src/t2elem.cc b/tracy/tracy/src/t2elem.cc
index 6aadc6705f852092ef57d950fa2a9396bbd90565..c38f5ca20b0e831d29470b65670dba51bedbee9d 100644
--- a/tracy/tracy/src/t2elem.cc
+++ b/tracy/tracy/src/t2elem.cc
@@ -11,7 +11,7 @@
double c_1, d_1, c_2, d_2, cl_rad, q_fluct;
double I2, I4, I5, dcurly_H, dI4;
ElemFamType ElemFam[Elem_nFamMax];
-CellType Cell[Cell_nLocMax + 1];
+CellType Cell[Cell_nLocMax+1];
// for IBS
int i_, j_;
@@ -26,24 +26,23 @@ double **C_;
Purpose:
Global to local coordinates
- ****************************************************************************/
+****************************************************************************/
template<typename T>
-void GtoL(ss_vect<T> &X, Vector2 &S, Vector2 &R, const double c0,
- const double c1, const double s1) {
+void GtoL(ss_vect<T> &X, Vector2 &S, Vector2 &R,
+ const double c0, const double c1, const double s1)
+{
ss_vect<T> x1;
/* Simplified rotated p_rot */
- X[px_] += c1;
- X[3] += s1;
+ X[px_] += c1; X[3] += s1;
/* Translate */
- X[x_] -= S[X_];
- X[y_] -= S[Y_];
+ X[x_] -= S[X_]; X[y_] -= S[Y_];
/* Rotate */
x1 = X;
- X[x_] = R[X_] * x1[x_] + R[Y_] * x1[y_];
- X[px_] = R[X_] * x1[px_] + R[Y_] * x1[py_];
- X[y_] = -R[Y_] * x1[x_] + R[X_] * x1[y_];
- X[py_] = -R[Y_] * x1[px_] + R[X_] * x1[py_];
+ X[x_] = R[X_]*x1[x_] + R[Y_]*x1[y_];
+ X[px_] = R[X_]*x1[px_] + R[Y_]*x1[py_];
+ X[y_] = -R[Y_]*x1[x_] + R[X_]*x1[y_];
+ X[py_] = -R[Y_]*x1[px_] + R[X_]*x1[py_] ;
/* Simplified p_rot */
X[px_] -= c0;
}
@@ -55,20 +54,21 @@ void GtoL(ss_vect<T> &X, Vector2 &S, Vector2 &R, const double c0,
Purpose:
Local to global coordinates
- ****************************************************************************/
+****************************************************************************/
template<typename T>
-void LtoG(ss_vect<T> &X, Vector2 &S, Vector2 &R, double c0, double c1,
- double s1) {
+void LtoG(ss_vect<T> &X, Vector2 &S, Vector2 &R,
+ double c0, double c1, double s1)
+{
ss_vect<T> x1;
/* Simplified p_rot */
X[px_] -= c0;
/* Rotate */
x1 = X;
- X[x_] = R[X_] * x1[x_] - R[Y_] * x1[y_];
- X[px_] = R[X_] * x1[px_] - R[Y_] * x1[py_];
- X[y_] = R[Y_] * x1[x_] + R[X_] * x1[y_];
- X[py_] = R[Y_] * x1[px_] + R[X_] * x1[py_];
+ X[x_] = R[X_]*x1[x_] - R[Y_]*x1[y_];
+ X[px_] = R[X_]*x1[px_] - R[Y_]*x1[py_];
+ X[y_] = R[Y_]*x1[x_] + R[X_]*x1[y_];
+ X[py_] = R[Y_]*x1[px_] + R[X_]*x1[py_];
/* Translate */
X[x_] += S[X_];
X[y_] += S[Y_];
@@ -81,16 +81,17 @@ void LtoG(ss_vect<T> &X, Vector2 &S, Vector2 &R, double c0, double c1,
Purpose:
Get the longitudinal momentum ps
- **********************************************************/
+**********************************************************/
template<typename T>
-inline T get_p_s(ss_vect<T> &x) {
+inline T get_p_s(ss_vect<T> &x)
+{
T p_s, p_s2;
if (!globval.H_exact)
- p_s = 1.0 + x[delta_];
+ p_s = 1.0+x[delta_];
else {
- p_s2 = sqr(1.0 + x[delta_]) - sqr(x[px_]) - sqr(x[py_]);
+ p_s2 = sqr(1.0+x[delta_]) - sqr(x[px_]) - sqr(x[py_]);
if (p_s2 >= 0.0)
p_s = sqrt(p_s2);
else {
@@ -100,7 +101,7 @@ inline T get_p_s(ss_vect<T> &x) {
exit_(1);
}
}
- return (p_s);
+ return(p_s);
}
/****************************************************************************/
@@ -111,9 +112,10 @@ inline T get_p_s(ss_vect<T> &x) {
If H_exact = false, using approximation Hamiltonian:
- px^2+py^2
- H(x,y,l,px,py,delta) = -------------
- 2(1+delta)
+ px^2+py^2
+ H(x,y,l,px,py,delta) = -------------
+ 2(1+delta)
+
Otherwise, using exact Hamiltonian
@@ -135,16 +137,17 @@ inline T get_p_s(ss_vect<T> &x) {
****************************************************************************/
+
template<typename T>
void Drift(double L, ss_vect<T> &x) {
T u;
if (!globval.H_exact) {
- u = L / (1.0 + x[delta_]);
- x[ct_] += u * (sqr(x[px_]) + sqr(x[py_])) / (2.0 * (1.0 + x[delta_]));
+ u = L/(1.0+x[delta_]);
+ x[ct_] += u*(sqr(x[px_])+sqr(x[py_]))/(2.0*(1.0+x[delta_]));
} else {
- u = L / get_p_s(x);
- x[ct_] += u * (1.0 + x[delta_]) - L;
+ u = L/get_p_s(x);
+ x[ct_] += u*(1.0+x[delta_]) - L;
}
x[x_] += x[px_] * u;
x[y_] += x[py_] * u;
@@ -617,36 +620,36 @@ static double get_psi(double irho, double phi, double gap) {
/* template<typename T>
void thin_kick(int Order, double MB[], double L, double h_bend, double h_ref,ss_vect<T> &x)
- Purpose:
- Calculate multipole kick. The Hamiltonian is
+ Purpose:
+ Calculate multipole kick. The Hamiltonian is
- H = A + B where A and B are the kick part defined by
- 2 2
- px + py
- A(x,y,-l,px,py,dP) = ---------
- 2(1+dP)
- 2 2
- B(x,y,-l,px,py,dP) = -h*x*dP + 0.5 h x + int(Re(By+iBx)/Brho)
+ H = A + B where A and B are the kick part defined by
+ 2 2
+ px + py
+ A(x,y,-l,px,py,dP) = ---------
+ 2(1+dP)
+ 2 2
+ B(x,y,-l,px,py,dP) = -h*x*dP + 0.5 h x + int(Re(By+iBx)/Brho)
- so this is the appproximation for large ring
- the chromatic term is missing hx*A
+ so this is the appproximation for large ring
+ the chromatic term is missing hx*A
- The kick is given by
+ The kick is given by
- e L L delta L x e L
- Dp_x = - --- B_y + ------- - ----- , Dp_y = --- B_x
- p_0 rho rho^2 p_0
+ e L L delta L x e L
+ Dp_x = - --- B_y + ------- - ----- , Dp_y = --- B_x
+ p_0 rho rho^2 p_0
- where
- e 1
- --- = -----
- p_0 B rho
- ====
- \
- (B_y + iB_x) = B rho > (ia_n + b_n ) (x + iy)^n-1
- /
- ====
+ where
+ e 1
+ --- = -----
+ p_0 B rho
+ ====
+ \
+ (B_y + iB_x) = B rho > (ia_n + b_n ) (x + iy)^n-1
+ /
+ ====
Input:
Order maximum non zero multipole component
@@ -718,13 +721,13 @@ void thin_kick(int Order, double MB[], double L, double h_bend, double h_ref,
Compute edge focusing for a dipole
There is no radiation coming from the edge
The standard formula used is :
- irho
- px = px0 + ------ tan(phi) *x0
- 1 + dP
+ irho
+ px = px0 + ------ tan(phi) *x0
+ 1 + dP
- irho
- pz = pz0 - ------ tan(phi - psi) *z0
- 1 + dP
+ irho
+ pz = pz0 - ------ tan(phi - psi) *z0
+ 1 + dP
for psi definition see its function
@@ -819,28 +822,29 @@ void bend_fringe(double hb, ss_vect<T> &x) {
form. If keep up to the 4-th order nonlinear terms, the formula with goes to the
following:
(see E. Forest's book (Beam Dynamics: A New Attitude and Framework), p390):
- b2
- x = x0 +/- ---------- (x0^3 + 3*z0^2*x0)
- 12(1 + dP)
+ b2
+ x = x0 +/- ---------- (x0^3 + 3*z0^2*x0)
+ 12(1 + dP)
- b2
- px = px0 +/- ---------- (2*x0*y0*py0 - x0^2*px0 - y0^2*py0)
- 4(1 + dP)
+ b2
+ px = px0 +/- ---------- (2*x0*y0*py0 - x0^2*px0 - y0^2*py0)
+ 4(1 + dP)
- b2
- y = y0 -/+ ---------- (y0^3 + 3*x0^2*y0)
- 12(1 + dP)
+ b2
+ y = y0 -/+ ---------- (y0^3 + 3*x0^2*y0)
+ 12(1 + dP)
- b2
- py = py0 -/+ ---------- (2*x0*y0*px0 - y0^2*py0 - x0^2*py0)
- 4(1 + dP)
+ b2
+ py = py0 -/+ ---------- (2*x0*y0*px0 - y0^2*py0 - x0^2*py0)
+ 4(1 + dP)
- dP = dP0;
+ dP = dP0;
+
+ b2
+ tau = tau0 -/+ ----------- (y0^3*px - x0^3*py + 3*x0^2*y*py - 3*y0^2*x0*px)
+ 12(1 + dP)^2
- b2
- tau = tau0 -/+ ----------- (y0^3*px - x0^3*py + 3*x0^2*y*py - 3*y0^2*x0*px)
- 12(1 + dP)^2
Input:
b2 = quadrupole strength
@@ -1862,22 +1866,10 @@ void GtoL_M(Matrix &X, Vector2 &T) {
Matrix R;
/* Rotate */
- R[0][0] = T[0];
- R[0][1] = 0.0;
- R[0][2] = T[1];
- R[0][3] = 0.0;
- R[1][0] = 0.0;
- R[1][1] = T[0];
- R[1][2] = 0.0;
- R[1][3] = T[1];
- R[2][0] = -T[1];
- R[2][1] = 0.0;
- R[2][2] = T[0];
- R[2][3] = 0.0;
- R[3][0] = 0.0;
- R[3][1] = -T[1];
- R[3][2] = 0.0;
- R[3][3] = T[0];
+ R[0][0] = T[0]; R[0][1] = 0.0; R[0][2] = T[1]; R[0][3] = 0.0;
+ R[1][0] = 0.0; R[1][1] = T[0]; R[1][2] = 0.0; R[1][3] = T[1];
+ R[2][0] = -T[1]; R[2][1] = 0.0; R[2][2] = T[0]; R[2][3] = 0.0;
+ R[3][0] = 0.0; R[3][1] = -T[1]; R[3][2] = 0.0; R[3][3] = T[0];
MulLMat(4L, R, X);
}
@@ -1885,22 +1877,10 @@ void LtoG_M(Matrix &X, Vector2 &T) {
Matrix R;
/* Rotate */
- R[0][0] = T[0];
- R[0][1] = 0.0;
- R[0][2] = -T[1];
- R[0][3] = 0.0;
- R[1][0] = 0.0;
- R[1][1] = T[0];
- R[1][2] = 0.0;
- R[1][3] = -T[1];
- R[2][0] = T[1];
- R[2][1] = 0.0;
- R[2][2] = T[0];
- R[2][3] = 0.0;
- R[3][0] = 0.0;
- R[3][1] = T[1];
- R[3][2] = 0.0;
- R[3][3] = T[0];
+ R[0][0] = T[0]; R[0][1] = 0.0; R[0][2] = -T[1]; R[0][3] = 0.0;
+ R[1][0] = 0.0; R[1][1] = T[0]; R[1][2] = 0.0; R[1][3] = -T[1];
+ R[2][0] = T[1]; R[2][1] = 0.0; R[2][2] = T[0]; R[2][3] = 0.0;
+ R[3][0] = 0.0; R[3][1] = T[1]; R[3][2] = 0.0; R[3][3] = T[0];
MulLMat(4, R, X);
}
@@ -1979,15 +1959,9 @@ static void make3by3(Matrix &A, double a11, double a12, double a13, double a21,
*/
UnitMat(ss_dim, A); /* set matrix to unit 3x3 matrix */
- A[0][0] = a11;
- A[0][1] = a12;
- A[0][2] = a13;
- A[1][0] = a21;
- A[1][1] = a22;
- A[1][2] = a23;
- A[2][0] = a31;
- A[2][1] = a32;
- A[2][2] = a33;
+ A[0][0] = a11; A[0][1] = a12; A[0][2] = a13;
+ A[1][0] = a21; A[1][1] = a22; A[1][2] = a23;
+ A[2][0] = a31; A[2][1] = a32; A[2][2] = a33;
}
/****************************************************************************
@@ -2024,20 +1998,13 @@ static void make4by5(Matrix &A, double a11, double a12, double a15, double a21,
double a22, double a25, double a33, double a34, double a35, double a43,
double a44, double a45) {
UnitMat(ss_dim, A); /* Initializes A to identity matrix */
- A[0][0] = a11;
- A[0][1] = a12;
- A[0][4] = a15;
- A[1][0] = a21;
- A[1][1] = a22;
- A[1][4] = a25;
- A[2][2] = a33;
- A[2][3] = a34;
- A[2][4] = a35;
- A[3][2] = a43;
- A[3][3] = a44;
- A[3][4] = a45;
+ A[0][0] = a11; A[0][1] = a12; A[0][4] = a15;
+ A[1][0] = a21; A[1][1] = a22; A[1][4] = a25;
+ A[2][2] = a33; A[2][3] = a34; A[2][4] = a35;
+ A[3][2] = a43; A[3][3] = a44; A[3][4] = a45;
}
+
static void mergeto4by5(Matrix &A, Matrix &AH, Matrix &AV) {
/*
merges two 3x3 matrices AH (H-plane) and AV (V-plane) into one
@@ -2084,13 +2051,17 @@ void Drift_SetMatrix(int Fnum1, int Knum1) {
cellp = &Cell[ElemFam[Fnum1 - 1].KidList[Knum1 - 1]];
elemp = &cellp->Elem;
D = elemp->D;
- L = elemp->PL / (1.0 + globval.dPparticle); /* L = L / (1 + dP) */
- make4by5(D->D55, 1.0, L, 0.0, 0.0, 1.0, 0.0, 1.0, L, 0.0, 0.0, 1.0, 0.0);
+ L = elemp->PL/(1.0+globval.dPparticle); /* L = L / (1 + dP) */
+ make4by5(D->D55,
+ 1.0, L, 0.0, 0.0, 1.0, 0.0,
+ 1.0, L, 0.0, 0.0, 1.0, 0.0);
}
static void driftmat(Matrix &ah, double L) {
L /= 1 + globval.dPparticle;
- make4by5(ah, 1.0, L, 0.0, 0.0, 1.0, 0.0, 1.0, L, 0.0, 0.0, 1.0, 0.0);
+ make4by5(ah,
+ 1.0, L, 0.0, 0.0, 1.0, 0.0,
+ 1.0, L, 0.0, 0.0, 1.0, 0.0);
}
static void quadmat(Matrix &ahv, double L, double k) {
@@ -2153,49 +2124,49 @@ static void quadmat(Matrix &ahv, double L, double k) {
/****************************************************************************/
/* static void bendmat(vector *M, double L, double irho, double phi1,
- double phi2, double gap, double k)
+ double phi2, double gap, double k)
- Purpose: called by Mpole_Setmatrix
+ Purpose: called by Mpole_Setmatrix
- For a quadrupole see quadmat routine for explanation
+ For a quadrupole see quadmat routine for explanation
- For a dipole
+ For a dipole
- (1 0 0)
- Edge(theta) = (h*tan(theta) 1 0)
- (0 0 1)
+ (1 0 0)
+ Edge(theta) = (h*tan(theta) 1 0)
+ (0 0 1)
- (1 0 0)
- Edge(theta) = (-h*tan(theta-psi) 1 0)
- (0 0 1)
+ (1 0 0)
+ Edge(theta) = (-h*tan(theta-psi) 1 0)
+ (0 0 1)
- 2
- K1*gap*h*(1 + sin phi)
- psi = -----------------------, K1 = 1/2
- cos phi
+ 2
+ K1*gap*h*(1 + sin phi)
+ psi = -----------------------, K1 = 1/2
+ cos phi
- Input:
- L : length [m]
- irho: 1/rho [1/m]
- phi1: entrance edge angle [degres]
- phi2: exit edge angle [degres]
- K : gradient = n/Rho
+ Input:
+ L : length [m]
+ irho: 1/rho [1/m]
+ phi1: entrance edge angle [degres]
+ phi2: exit edge angle [degres]
+ K : gradient = n/Rho
- Output:
- M transfer matrix
+ Output:
+ M transfer matrix
- Return:
- none
+ Return:
+ none
- Global variables:
- none
+ Global variables:
+ none
- Specific functions:
- quadmat, make3by3
- UnitMat, MulRMat, psi
+ Specific functions:
+ quadmat, make3by3
+ UnitMat, MulRMat, psi
- Comments:
- none
+ Comments:
+ none
****************************************************************************/
static void bendmat(Matrix &M, double L, double irho, double phi1, double phi2,
@@ -2236,28 +2207,29 @@ static void bendmat(Matrix &M, double L, double irho, double phi1, double phi2,
if (k == 0.0) {
/* simple vertical dipole magnet */
- /*
- H-plane
- 2 2 2
- px + py 2 x
- H = -------- - h*x*dP + h ---
- 2*(1+dP) 2
-
- 2 2
- dx h h*dP
- --2 + ---- x = -----
- ds 1+dP 1+dP
-
-
- let be u = Lh/sqrt(1+dP) then the transfert matrix becomes:
-
- ( sin(u) 1- cos(u) )
- ( cos(u) -------------- ----------------- )
- ( h sqrt(1+dP) h )
- ( -sin(u)*sqrt(1+dP)*h cos(u) sin(u)*sqrt(1+dP) )
- ( 0 0 1 )
-
- */
+ /* simple vertical dipole magnet */
+ /*
+ H-plane
+ 2 2 2
+ px + py 2 x
+ H = -------- - h*x*dP + h ---
+ 2*(1+dP) 2
+
+ 2 2
+ dx h h*dP
+ --2 + ---- x = -----
+ ds 1+dP 1+dP
+
+
+ let be u = Lh/sqrt(1+dP) then the transfert matrix becomes:
+
+ ( sin(u) 1- cos(u) )
+ ( cos(u) -------------- ----------------- )
+ ( h sqrt(1+dP) h )
+ ( -sin(u)*sqrt(1+dP)*h cos(u) sin(u)*sqrt(1+dP) )
+ ( 0 0 1 )
+
+ */
c = cos(r);
s = sin(r);
make3by3(ah, c, s / (irho * scoef), (1.0 - c) / irho,
@@ -2282,29 +2254,29 @@ static void bendmat(Matrix &M, double L, double irho, double phi1, double phi2,
sk = sqrt(afk);
p = L * sk / scoef;
if (fk < 0.0) {
- /*
- H-plane
- 2 2 2
- px + py 2 x
- H = -------- - h*x*dP + (k+h ) ---
- 2*(1+dP) 2
+ /*
+ H-plane
+ 2 2 2
+ px + py 2 x
+ H = -------- - h*x*dP + (k+h ) ---
+ 2*(1+dP) 2
- 2 2
- dx k+h h*dP
- --2 + ---- x = -----
- ds 1+dP 1+dP
+ 2 2
+ dx k+h h*dP
+ --2 + ---- x = -----
+ ds 1+dP 1+dP
- let be u = Lsqrt(|h*h+b2|)/sqrt(1+dP)
- then the transfert matrix becomes:
+ let be u = Lsqrt(|h*h+b2|)/sqrt(1+dP)
+ then the transfert matrix becomes:
- ( sin(u) 1- cos(u) )
- ( cos(u) -------------- ----------------- )
- ( sk*sqrt(1+dP) |k+h*h|*h )
- ( -sin(u)*sqrt(1+dP)*sk cos(u) h*sin(u)*sqrt(1+dP)/sk)
- ( 0 0 1 )
+ ( sin(u) 1- cos(u) )
+ ( cos(u) -------------- ----------------- )
+ ( sk*sqrt(1+dP) |k+h*h|*h )
+ ( -sin(u)*sqrt(1+dP)*sk cos(u) h*sin(u)*sqrt(1+dP)/sk)
+ ( 0 0 1 )
- */
+ */
c = cos(p);
s = sin(p);
make3by3(ah, c, s / sk / scoef, irho * (1.0 - c) / (coef * afk),
@@ -2331,21 +2303,21 @@ static void bendmat(Matrix &M, double L, double irho, double phi1, double phi2,
}
/* Edge focusing, no effect due to gap between AU and AD */
- /*
- (1 0 0)
- Edge(theta) = (h*tan(theta) 1 0)
- (0 0 1)
-
- (1 0 0)
- Edge(theta) = (-h*tan(theta-psi) 1 0)
- (0 0 1)
-
- 2
- K1*gap*h*(1 + sin phi)
- psi = -----------------------, K1 = 1/2
- cos phi
-
- */
+ /*
+ (1 0 0)
+ Edge(theta) = (h*tan(theta) 1 0)
+ (0 0 1)
+
+ (1 0 0)
+ Edge(theta) = (-h*tan(theta-psi) 1 0)
+ (0 0 1)
+
+ 2
+ K1*gap*h*(1 + sin phi)
+ psi = -----------------------, K1 = 1/2
+ cos phi
+
+ */
if (phi1 != 0.0 || gap > 0.0) {
UnitMat(3L, edge);
edge[1][0] = irho * tan(dtor(phi1));