From 791825846e0f3317b5b7b0df0d7c498b4bf78537 Mon Sep 17 00:00:00 2001
From: nadolski <nadolski@9a6e40ed-f3a0-4838-9b4a-bf418f78e88d>
Date: Sun, 30 Jun 2013 05:12:34 +0000
Subject: [PATCH] Cleaning the code
---
tracy/tracy/src/t2cell.cc | 3 +-
tracy/tracy/src/t2elem.cc | 275 +++++++++++++++++++++++++++++---------
tracy/tracy/src/t2lat.cc | 23 +---
3 files changed, 220 insertions(+), 81 deletions(-)
diff --git a/tracy/tracy/src/t2cell.cc b/tracy/tracy/src/t2cell.cc
index f23f391..1f0be8e 100644
--- a/tracy/tracy/src/t2cell.cc
+++ b/tracy/tracy/src/t2cell.cc
@@ -259,9 +259,8 @@ void Cell_Pass(const long i0, const long i1, ss_vect<T> &x, long &lastpos)
{
lastpos = i1;
for (i = i0; i <= i1; i++)
- // for (i = i0+1L; i <= i1; i++)
{
- Elem_Pass(i, x);
+ Elem_Pass(i, x);
if (!CheckAmpl(x, i))
{ lastpos = i; break; }
}
diff --git a/tracy/tracy/src/t2elem.cc b/tracy/tracy/src/t2elem.cc
index 4af1255..49495d7 100644
--- a/tracy/tracy/src/t2elem.cc
+++ b/tracy/tracy/src/t2elem.cc
@@ -32,7 +32,7 @@ double **C_;
****************************************************************************/
template<typename T>
void GtoL(ss_vect<T> &X, Vector2 &S, Vector2 &R,
- const double c0, const double c1, const double s1)
+ const double c0, const double c1, const double s1)
{
ss_vect<T> x1;
@@ -60,7 +60,7 @@ void GtoL(ss_vect<T> &X, Vector2 &S, Vector2 &R,
****************************************************************************/
template<typename T>
void LtoG(ss_vect<T> &X, Vector2 &S, Vector2 &R,
- double c0, double c1, double s1)
+ double c0, double c1, double s1)
{
ss_vect<T> x1;
@@ -122,8 +122,8 @@ 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)
+ H(x,y,l,px,py,delta) = -------------
+ 2(1+delta)
Otherwise, using exact Hamiltonian
@@ -194,10 +194,9 @@ void zero_mat(const int n, double** A) {
****************************************************************************/
void identity_mat(const int n, double** A) {
- int i, j;
- for (i = 1; i <= n; i++)
- for (j = 1; j <= n; j++)
+ for (int i = 1; i <= n; i++)
+ for (int j = 1; j <= n; j++)
A[i][j] = (i == j) ? 1.0 : 0.0;
}
@@ -210,7 +209,7 @@ void identity_mat(const int n, double** A) {
****************************************************************************/
double det_mat(const int n, double **A) {
- int i, *indx;
+ int *indx;
double **U, d;
indx = ivector(1, n);
@@ -219,7 +218,7 @@ double det_mat(const int n, double **A) {
dmcopy(A, n, n, U);
dludcmp(U, n, indx, &d);
- for (i = 1; i <= n; i++)
+ for (int i = 1; i <= n; i++)
d *= U[i][i];
free_dmatrix(U, 1, n, 1, n);
@@ -237,11 +236,10 @@ double det_mat(const int n, double **A) {
****************************************************************************/
double trace_mat(const int n, double **A) {
- int i;
double d;
d = 0.0;
- for (i = 1; i <= n; i++)
+ for (int i = 1; i <= n; i++)
d += A[i][i];
return d;
@@ -652,7 +650,7 @@ static double get_psi(double irho, double phi, double gap) {
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
@@ -1699,7 +1697,7 @@ void FieldMap_Pass(CellType &Cell, ss_vect<T> &X) {
}
/********************************************************************
- void Insertion_Pass(CellType &Cell, ss_vect<T> &x)
+ void Insertion_Pass(CellType &Cell, ss_vect<T> &X)
Purpose:
Track vector x through an insertion
@@ -1730,85 +1728,238 @@ void FieldMap_Pass(CellType &Cell, ss_vect<T> &X) {
*******************************************************************/
template<typename T>
-void Insertion_Pass(CellType &Cell, ss_vect<T> &x) {
+void Insertion_Pass(CellType &Cell, ss_vect<T> &X) {
elemtype *elemp;
double LN = 0.0;
+ double L1 = 0.0, L2 = 0.0, K1 = 0.0, K2 = 0.0;
T tx1, tz1; /* thetax and thetaz retrieved from
interpolation routine First order kick*/
T tx2, tz2; /* thetax and thetaz retrieved from
interpolation routine Second order Kick */
T d;
+
double alpha0 = 0.0; // 1/ brh0
double alpha02 = 0.0; // alpha square
int Nslice = 0;
int i = 0;
bool outoftable = false;
- elemp = &Cell.Elem;
- Nslice = elemp->ID->PN;
- alpha0 = c0 / globval.Energy * 1E-9 * elemp->ID->scaling1;
- alpha02 = (c0 / globval.Energy * 1E-9) * (c0 / globval.Energy * 1E-9)
- * elemp->ID->scaling2;
+ elemp = &Cell.Elem;
+ Nslice = elemp->ID->PN;
+ alpha0 = c0 / globval.Energy * 1E-9 * elemp->ID->scaling1;
+ alpha02 = sqr(c0 / globval.Energy * 1E-9) * elemp->ID->scaling2;
- // /* Global -> Local */
- // GtoL(X, Cell->dS, Cell->dT, 0.0, 0.0, 0.0);
+ /* Global -> Local */
+ //GtoL(X, elemp->ID->dS, elemp->ID->dT, 0.0, 0.0, 0.0);
// Nslice*2*1/2 drifts, Nslice kicks
LN = elemp->PL / Nslice;
- Drift(LN/2.0, x);
+
+
+ if ((elemp->ID->firstorder && elemp->ID->scaling1 == 0.0) &
+ (elemp->ID->secondorder || elemp->ID->scaling2 == 0.0)){
+ Drift(elemp->PL, X); // only a drift if kick is zero
+ if (trace) fprintf(stdout, "Insertion_Pass: used a drift since scaling factors are zeros for elem %s\n",
+ elemp->PName);
+ }
+ else{ // Integrator scheme
+ switch (elemp->ID->Pmethod) {
+ case Meth_Second:
+ for (i = 1; i <= Nslice; i++) {
+ Drift(LN/2.0, X);
+ // fprintf(stdout, "Second order L= %f, N= %d LN=%f S1=%f S2 = %f\n",
+ //elemp->PL, Nslice, LN, elemp->ID->scaling1, elemp->ID->scaling1);
+
+ // First order kick map
+ if (elemp->ID->firstorder && elemp->ID->scaling1 != 0.0) {
+ if (!elemp->ID->linear) {
+ //cout << "InsertionPass: Spline\n";
+ SplineInterpolation2(X[x_], X[y_], tx1, tz1, Cell, outoftable,
+ 1);
+ } else{
+ LinearInterpolation2(X[x_], X[y_], tx1, tz1, Cell, outoftable,
+ 1);
+ }
+ if (outoftable) {
+ X[x_] = NAN;
+ return;
+ }
+
+ d = alpha0 / Nslice;
+ X[px_] += d * tx1;
+ X[py_] += d * tz1;
+ }
+
+ // Second order kick map
+ if (elemp->ID->secondorder && elemp->ID->scaling2 != 0.0) {
+ if (!elemp->ID->linear) {
+ //cout << "InsertionPass: Spline\n";
+ SplineInterpolation2(X[x_], X[y_], tx2, tz2, Cell, outoftable, 2);
+ } else {
+ //cout << "InsertionPass: Linear\n";
+ LinearInterpolation2(X[x_], X[y_], tx2, tz2, Cell, outoftable,
+ 2);
+ }
+ if (outoftable) {
+ X[x_] = NAN;
+ return;
+ }
+ d = alpha02 / Nslice / (1.0 + X[delta_]);
+ X[px_] += d * tx2;
+ X[py_] += d * tz2;
+ }
+ if (i != Nslice)
+ Drift(LN, X);
+ else
+ Drift(LN/2.0, X);
+ } // for
+ break;
+ case Meth_Fourth: /* 4-th order integrator */
+ L1 = c_1 * LN;
+ L2 = c_2 * LN;
+ K1 = d_1;
+ K2 = d_2;
+
+// Drift(L1, X);
+ for (i = 1; i <= Nslice; i++) {
+ Drift(L1, X);
+ // First order kick map
+ if (elemp->ID->firstorder && elemp->ID->scaling1 != 0.0) {
+ if (!elemp->ID->linear) {
+ //cout << "InsertionPass: Spline\n";
+ SplineInterpolation2(X[x_], X[y_], tx1, tz1, Cell, outoftable,
+ 1);
+ } else{
+ LinearInterpolation2(X[x_], X[y_], tx1, tz1, Cell, outoftable,
+ 1);
+ }
+ if (outoftable) {
+ X[x_] = NAN;
+ return;
+ }
- for (i = 1; i <= Nslice; i++) {
- // First order kick map
- if (elemp->ID->firstorder) {
+ d = alpha0 / Nslice;
+ X[px_] += K1*d * tx1;
+ X[py_] += K1*d * tz1;
+ }
+ // Second order kick map
+ if (elemp->ID->secondorder && elemp->ID->scaling2 != 0.0) {
if (!elemp->ID->linear) {
//cout << "InsertionPass: Spline\n";
- SplineInterpolation2(x[x_], x[y_], tx1, tz1, Cell, outoftable,
+ SplineInterpolation2(X[x_], X[y_], tx2, tz2, Cell, outoftable,
+ 2);
+ } else {
+ //cout << "InsertionPass: Linear\n";
+ LinearInterpolation2(X[x_], X[y_], tx2, tz2, Cell, outoftable,
+ 2);
+ }
+ if (outoftable) {
+ X[x_] = NAN;
+ return;
+ }
+ d = alpha02 / Nslice / (1.0 + X[delta_]);
+ X[px_] += K1*d * tx2;
+ X[py_] += K1*d * tz2;
+ }
+
+ Drift(L2, X);
+
+ // First order kick map
+ if (elemp->ID->firstorder && elemp->ID->scaling1 != 0.0) {
+ if (!elemp->ID->linear) {
+ //cout << "InsertionPass: Spline\n";
+ SplineInterpolation2(X[x_], X[y_], tx1, tz1, Cell, outoftable,
1);
} else{
- LinearInterpolation2(x[x_], x[y_], tx1, tz1, Cell, outoftable,
+ LinearInterpolation2(X[x_], X[y_], tx1, tz1, Cell, outoftable,
1);
}
if (outoftable) {
- x[x_] = NAN;
+ X[x_] = NAN;
return;
}
d = alpha0 / Nslice;
- x[px_] += d * tx1;
- x[py_] += d * tz1;
+ X[px_] += K2*d * tx1;
+ X[py_] += K2*d * tz1;
}
- // Second order kick map
- if (elemp->ID->secondorder) {
+ // Second order kick map
+ if (elemp->ID->secondorder && elemp->ID->scaling2 != 0.0) {
if (!elemp->ID->linear) {
//cout << "InsertionPass: Spline\n";
- SplineInterpolation2(x[x_], x[y_], tx2, tz2, Cell, outoftable,
+ SplineInterpolation2(X[x_], X[y_], tx2, tz2, Cell, outoftable,
2);
} else {
//cout << "InsertionPass: Linear\n";
- LinearInterpolation2(x[x_], x[y_], tx2, tz2, Cell, outoftable,
+ LinearInterpolation2(X[x_], X[y_], tx2, tz2, Cell, outoftable,
2);
}
if (outoftable) {
- x[x_] = NAN;
+ X[x_] = NAN;
return;
}
+ d = alpha02 / Nslice / (1.0 + X[delta_]);
+ X[px_] += K2*d * tx2;
+ X[py_] += K2*d * tz2;
+ }
- d = alpha02 / Nslice / (1.0 + x[delta_]);
- x[px_] += d * tx2;
- x[py_] += d * tz2;
+ Drift(L2, X);
+
+ // First order kick map
+ if (elemp->ID->firstorder && elemp->ID->scaling1 != 0.0) {
+ if (!elemp->ID->linear) {
+ //cout << "InsertionPass: Spline\n";
+ SplineInterpolation2(X[x_], X[y_], tx1, tz1, Cell, outoftable,
+ 1);
+ } else{
+ LinearInterpolation2(X[x_], X[y_], tx1, tz1, Cell, outoftable,
+ 1);
+ }
+ if (outoftable) {
+ X[x_] = NAN;
+ return;
+ }
+
+ d = alpha0 / Nslice;
+ X[px_] += K1*d * tx1;
+ X[py_] += K1*d * tz1;
}
- if (i != Nslice)
- Drift(LN, x);
- else
- Drift(LN/2.0, x);
- }
- //CopyVec(6L, x, Cell->BeamPos);
+ // Second order kick map
+ if (elemp->ID->secondorder && elemp->ID->scaling2 != 0.0){
+ if (!elemp->ID->linear) {
+ //cout << "InsertionPass: Spline\n";
+ SplineInterpolation2(X[x_], X[y_], tx2, tz2, Cell, outoftable,
+ 2);
+ } else {
+ //cout << "InsertionPass: Linear\n";
+ LinearInterpolation2(X[x_], X[y_], tx2, tz2, Cell, outoftable,
+ 2);
+ }
+ if (outoftable) {
+ X[x_] = NAN;
+ return;
+ }
+ d = alpha02 / Nslice / (1.0 + X[delta_]);
+ X[px_] += K1*d * tx2;
+ X[py_] += K1*d * tz2;
+ }
+
+// if (i != Nslice)
+// Drift(L1*2.0, X);
+// else
+ Drift(L1, X);
+ break;
+ } // for
- /* Local -> Global */
- //LtoG(X, Cell->dS, Cell->dT, 0.0, 0.0, 0.0);
+ } // switch
+ } // if
+ //CopyVec(6L, X, Cell->BeamPos);
+
+ /* Local -> Global */
+ //LtoG(X, Cell->dS, Cell->dT, 0.0, 0.0, 0.0);
}
template<typename T>
@@ -1833,7 +1984,7 @@ void sol_pass(const elemtype &elem, ss_vect<T> &x) {
x[px_] -= dpx;
x[py_] -= AyoBrho;
- x[y_] += 0.5 * hd * x[py_];
+ x[y_] += 0.5 * hd * x[py_];
x[ct_] += sqr(0.5) * hd * sqr(x[py_]) / (1.0 + x[delta_]);
AyoBrho = elem.Sol->BoBrho * x[x_] / 2.0;
@@ -2129,15 +2280,15 @@ void Drift_SetMatrix(int Fnum1, int Knum1) {
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);
+ 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);
+ 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) {
@@ -2678,9 +2829,8 @@ void Insertion_SetMatrix(int Fnum1, int Knum1) {
elemp = &cellp->Elem;
ID = elemp->ID;
Nslice = ID->PN;
- alpha0 = c0 / globval.Energy * 1E-9 * ID->scaling1;
- alpha02 = (c0 / globval.Energy * 1E-9) * c0 / globval.Energy * 1E-9 / (1.0
- + globval.dPparticle) * ID->scaling2;
+ alpha0 = c0 / globval.Energy * 1E-9 * ID->scaling1;
+ alpha02 = sqr(c0 / globval.Energy * 1E-9) / (1.0+ globval.dPparticle) * ID->scaling2;
UnitMat(6L, ID->D55);
UnitMat(6L, ID->K55);
@@ -2699,10 +2849,10 @@ void Insertion_SetMatrix(int Fnum1, int Knum1) {
/* quadrupole Kick matrix linearized around closed orbit */
if (!ID->linear) {
// SplineInterpDeriv3(cellp->BeamPos[0], cellp->BeamPos[2],
- // &DTHXDX, &DTHXDZ, &DTHZDX, &DTHZDZ, cellp);
+ // &DTHXDX, &DTHXDZ, &DTHZDX, &DTHZDZ, cellp);
} else {
// LinearInterpDeriv2(cellp->BeamPos[0], cellp->BeamPos[2],
- // &DTHXDX, &DTHXDZ, &DTHZDX, &DTHZDZ, cellp, 1);
+ // &DTHXDX, &DTHXDZ, &DTHZDX, &DTHZDZ, cellp, 1);
}
ID->K55[1][0] = ID->K55[1][0] + alpha0 * DTHXDX / Nslice;
ID->K55[1][2] = ID->K55[1][2] + alpha0 * DTHXDZ / Nslice;
@@ -2715,10 +2865,10 @@ void Insertion_SetMatrix(int Fnum1, int Knum1) {
/* quadrupole Kick matrix linearized around closed orbit */
if (!ID->linear) {
// SplineInterpDeriv3(cellp->BeamPos[0], cellp->BeamPos[2],
- // &DTHXDX, &DTHXDZ, &DTHZDX, &DTHZDZ, cellp);
+ // &DTHXDX, &DTHXDZ, &DTHZDX, &DTHZDZ, cellp);
} else {
// LinearInterpDeriv2(cellp->BeamPos[0], cellp->BeamPos[2],
- // &DTHXDX, &DTHXDZ, &DTHZDX, &DTHZDZ, cellp, 2);
+ // &DTHXDX, &DTHXDZ, &DTHZDX, &DTHZDZ, cellp, 2);
}
ID->K55[1][0] = ID->K55[1][0] + alpha02 * DTHXDX / Nslice;
ID->K55[1][2] = ID->K55[1][2] + alpha02 * DTHXDZ / Nslice;
@@ -2738,8 +2888,8 @@ void Insertion_SetMatrix(int Fnum1, int Knum1) {
// {
// L = elemp->PL/(1.0 + globval.dPparticle); /* L = L/(1 + dP) */
// make4by5(ID->KD55,
- // 1.0, L, 0.0, 0.0, 1.0, 0.0,
- // 1.0, L, 0.0, 0.0, 1.0, 0.0);
+ // 1.0, L, 0.0, 0.0, 1.0, 0.0,
+ // 1.0, L, 0.0, 0.0, 1.0, 0.0);
// }
}
@@ -2783,6 +2933,7 @@ void Insertion_Pass_M(CellType &Cell, Vector &xref, Matrix &M) {
// {
MulLMat(5, elemp->ID->KD55, M); /* M<-KD55*M */
LinTrans(5, elemp->ID->KD55, xref);
+ cout << "Insertion_Pass_M Hello\n";
// }
// else
// {
@@ -4161,7 +4312,7 @@ void Mpole_DefPBsys(int Fnum1, int Knum1, int Order, double PBsys) {
}
/*********************************************************************
void Mpole_SetdS(int Fnum1, int Knum1)
-
+
Purpose:
Set misalignment error to the element with Fnum and Knum
@@ -4187,7 +4338,7 @@ void Mpole_SetdS(int Fnum1, int Knum1) {
/*********************************************************************
void Mpole_SetdT(int Fnum1, int Knum1)
-
+
Purpose:
Set rotation error to the element with Fnum and Knum
diff --git a/tracy/tracy/src/t2lat.cc b/tracy/tracy/src/t2lat.cc
index 9dec889..8175b9f 100644
--- a/tracy/tracy/src/t2lat.cc
+++ b/tracy/tracy/src/t2lat.cc
@@ -7,9 +7,9 @@ L. Nadolski SOLEIL 2002 Link to NAFF, Radia field maps
J. Bengtsson NSLS-II, BNL 2004 -
*/
-/* Current revision $Revision: 1.20 $
+/* Current revision $Revision: 1.21 $
On branch $Name: not supported by cvs2svn $
- Latest change $Date: 2012-01-27 16:47:39 $ by $Author: zhang $
+ Latest change $Date: 2013-06-30 05:12:34 $ by $Author: nadolski $
*/
@@ -3237,7 +3237,7 @@ static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
QKxV = 0e0;
QKS = 0e0;
k1 = 1; // number of slices of the lattice element
- k2 = 3; // 1 linear interpolation, 3 spline interpolation
+ k2 = 2; // 1 linear interpolation, 3 spline interpolation
dt = 0e0;
scaling1 = 1.0; // scaling factor
scaling2 = 1.0; // scaling factor
@@ -3253,8 +3253,7 @@ static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
getest__(P_expset(SET, 1L << ((long) eql)), "<=> expected", &V);
//read the parameters setting from the lattice
- switch (sym1) {
-
+ switch (sym1) {
case nsym: /* Read number of slices */
k1 = abs((long) floor(EVAL_(&V)));
GetSym__(&V);
@@ -3291,7 +3290,7 @@ static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
GetSym__(&V);
break;
- case mthsym: // method for interpolation: 1 means linear 3 spline
+ case mthsym: // method for interpolation: 1 means linear >1 means spline
k2 = (long) floor(EVAL_(&V));
if (k2 != Meth_Linear)
globval.MatMeth = false;
@@ -3353,7 +3352,7 @@ static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
exit_(-1);
}
- if (k2 == 3 | k2 == 2) { // cubic interpolation
+ if (k2 > 1) { // cubic interpolation
WITH5->linear = false;
} else { // linear interpolation
WITH5->linear = true;
@@ -3380,16 +3379,6 @@ static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
Matrices4Spline(WITH5,2);}
}
- // to put somewhere
- // /** Free memory **/
- // free(tab1);
- // free(tab2);
- //
- // free_matrix(tx,1,nz,1,nx);
- // free_matrix(tz,1,nz,1,nx);
- // free_matrix(f2x,1,nz,1,nx);
- // free_matrix(f2z,1,nz,1,nx);
-
} else {
printf("Elem_nFamMax exceeded: %ld(%ld)\n", globval.Elem_nFam,
(long) Elem_nFamMax);
--
GitLab