Skip to content
Snippets Groups Projects
Commit 4f16249a authored by nadolski's avatar nadolski
Browse files

Format restaure from 1.6 hoping there is no side effect !!!

parent 60874daa
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
tracy/tracy/src/t2elem.cc
+ 191
− 219
View file @ 4f16249a
......@@ -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));
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

The place to be for your source code