Commit 751a9978 authored by ext-rogers_c's avatar ext-rogers_c
Browse files

Some code cleanup and testing; identified bug in routines for non-straight dipole

parent 410a1c3e
......@@ -65,7 +65,6 @@ MultipoleT::MultipoleT(const MultipoleT &right):
verticalApert_m(right.verticalApert_m),
horizApert_m(right.horizApert_m),
dummy() {
RefPartBunch_m = right.RefPartBunch_m;
}
......@@ -88,17 +87,16 @@ bool MultipoleT::apply(const Vector_t &R, const Vector_t &P,
/** If magnet is not straight go to coords along the magnet */
if(angle_m != 0.0) {
Vector_t X = R_prime;
R_prime = transformCoords(X);
R_prime = transformCoords(X);
}
if (insideAperture(R_prime)) {
B[0] = getBx(R_prime);
B[1] = getBz(R_prime);
B[2] = getBs(R_prime);
return false;
}
else {
B[1] = getBz(R_prime);
B[2] = getBs(R_prime);
return false;
} else {
for(int i = 0; i < 3; i++) {
B[i] = 0.0;
B[i] = 0.0;
}
return true;
}
......@@ -155,47 +153,45 @@ Vector_t MultipoleT::transformCoords(const Vector_t &R) {
// if radius not variable
if (not variableRadius_m) {
double radius = length_m / angle_m;
// Transform from Cartesian to Frenet-Seret along the magnet
double alpha = atan(R[2] / (R[0] + radius ));
if (alpha != 0.0) {
X[0] = R[2] / sin(alpha) - radius;
X[1] = R[1];
X[2] = radius * alpha;
}
else {
X = R;
}
}
// if radius is variable need to transform coordinates at each
// point along the trajectory
else {
// Transform from Cartesian to Frenet-Seret along the magnet
double alpha = atan(R[2] / (R[0] + radius ));
if (alpha != 0.0) {
X[0] = R[2] / sin(alpha) - radius;
X[1] = R[1];
X[2] = radius * alpha;
} else {
X = R;
}
} else {
// if radius is variable need to transform coordinates at each
// point along the trajectory
double deltaAlpha, S = 0.0, localRadius;
double stepSize = varStep_m; // mm -> has a big effect on tracking time
if (abs(R[2]) <= stepSize) {
return R; // no transformation around origin
}
Y = R;
Vector_t temp;
while (abs(Y[2]) > stepSize and getRadius(S) != -1) {
localRadius = getRadius(S);
deltaAlpha = stepSize / localRadius;
if (R[2] < 0) {
deltaAlpha *= - 1.; // rotate in the other direction
double stepSize = varStep_m; // mm -> has a big effect on tracking time
if (abs(R[2]) <= stepSize) {
return R; // no transformation around origin
}
Y = R;
Vector_t temp;
while (abs(Y[2]) > stepSize and getRadius(S) != -1) {
localRadius = getRadius(S);
deltaAlpha = stepSize / localRadius;
if (R[2] < 0) {
deltaAlpha *= - 1.; // rotate in the other direction
}
temp = Y;
// translation
temp[0] += localRadius * (1 - cos(deltaAlpha));
temp[2] -= localRadius * sin(deltaAlpha);
// + rotation along the ideal trajectory
Y[2] = temp[2] * cos(deltaAlpha) - temp[0] * sin(deltaAlpha);
Y[0] = temp[2] * sin(deltaAlpha) + temp[0] * cos(deltaAlpha);
S += localRadius * deltaAlpha;
// until we reach the actual point from Cartesian coordinates
// translation
temp[0] += localRadius * (1 - cos(deltaAlpha));
temp[2] -= localRadius * sin(deltaAlpha);
// + rotation along the ideal trajectory
Y[2] = temp[2] * cos(deltaAlpha) - temp[0] * sin(deltaAlpha);
Y[0] = temp[2] * sin(deltaAlpha) + temp[0] * cos(deltaAlpha);
S += localRadius * deltaAlpha;
// until we reach the actual point from Cartesian coordinates
}
X[0] = Y[0];
X[1] = Y[1];
X[2] = S;
}
X[0] = Y[0];
X[1] = Y[1];
X[2] = S;
}
return X;
}
......@@ -205,15 +201,17 @@ void MultipoleT::setTransProfile(unsigned int n, double dTn) {
transProfile_m.resize(n+1, 0.0);
}
transProfile_m[n] = dTn;
}
}
bool MultipoleT::setFringeField(double s0, double lambda_l, double lambda_r) {
fringeField_l.Tanh::setLambda(lambda_l);
fringeField_r.Tanh::setLambda(lambda_r);
fringeField_l.Tanh::setX0(s0);
fringeField_l.Tanh::setTanhDiffIndices(2 * maxOrder_m + 1);
fringeField_r.Tanh::setLambda(lambda_r);
fringeField_r.Tanh::setX0(s0);
fringeField_r.Tanh::setTanhDiffIndices(2 * maxOrder_m + 1);
return true;
}
......@@ -225,26 +223,25 @@ double MultipoleT::getBz(const Vector_t &R) {
if (angle_m == 0.0) {
// Straight geometry -> use corresponding field expansion
for(unsigned int n = 0; n <= maxOrder_m; n++) {
double f_n = 0.0;
for(unsigned int i = 0; i <= n; i++) {
f_n += gsl_sf_choose(n, i) * getTransDeriv(2 * i, R[0]) *
getFringeDeriv(2 * n - 2 * i, R[2]);
double f_n = 0.0;
for(unsigned int i = 0; i <= n; i++) {
f_n += gsl_sf_choose(n, i) * getTransDeriv(2 * i, R[0]) *
getFringeDeriv(2 * n - 2 * i, R[2]);
}
f_n *= gsl_sf_pow_int(-1.0, n);
Bz += gsl_sf_pow_int(R[1], 2 * n) / gsl_sf_fact(2 * n) * f_n;
}
f_n *= gsl_sf_pow_int(-1.0, n);
Bz += gsl_sf_pow_int(R[1], 2 * n) / gsl_sf_fact(2 * n) * f_n;
}
}
else {
} else {
if (variableRadius_m == true and getFringeDeriv(0, R[2]) < 1.0e-12) {
// return 0 if end of fringe field is reached
// this is to avoid functions being called at infinite radius
return 0.0;
}
}
// Curved geometry -> use corresponding field expansion
for(unsigned int n = 0; n <= maxOrder_m; n++) {
double f_n = getFn(n, R[0], R[2]);
Bz += gsl_sf_pow_int(R[1], 2 * n) / gsl_sf_fact(2 * n) * f_n;
}
Bz += gsl_sf_pow_int(R[1], 2 * n) / gsl_sf_fact(2 * n) * f_n;
}
}
return Bz;
}
......@@ -257,28 +254,27 @@ double MultipoleT::getBx(const Vector_t &R) {
if (angle_m == 0.0) {
// Straight geometry -> use corresponding field expansion
for(unsigned int n = 0; n <= maxOrder_m; n++) {
double f_n = 0.0;
for(unsigned int i = 0; i <= n; i++) {
f_n += gsl_sf_choose(n, i) * getTransDeriv(2 * i + 1, R[0]) *
getFringeDeriv(2 * n - 2 * i, R[2]);
double f_n = 0.0;
for(unsigned int i = 0; i <= n; i++) {
f_n += gsl_sf_choose(n, i) * getTransDeriv(2 * i + 1, R[0]) *
getFringeDeriv(2 * n - 2 * i, R[2]);
}
f_n *= gsl_sf_pow_int(-1.0, n);
Bx += gsl_sf_pow_int(R[1], 2 * n + 1) /
gsl_sf_fact(2 * n + 1) * f_n;
}
f_n *= gsl_sf_pow_int(-1.0, n);
Bx += gsl_sf_pow_int(R[1], 2 * n + 1) /
gsl_sf_fact(2 * n + 1) * f_n;
}
}
else {
} else {
if (variableRadius_m == true and getFringeDeriv(0, R[2]) < 1.0e-12) {
// return 0 if end of fringe field is reached
// this is to avoid functions being called at infinite radius
return 0.0;
}
}
// Curved geometry -> use corresponding field expansion
for(unsigned int n = 0; n <= maxOrder_m; n++) {
double partialX_fn = getFnDerivX(n, R[0], R[2]);
Bx += gsl_sf_pow_int(R[1], 2 * n + 1) / gsl_sf_fact(2 * n + 1);
Bx *= partialX_fn;
}
Bx += gsl_sf_pow_int(R[1], 2 * n + 1) / gsl_sf_fact(2 * n + 1);
Bx *= partialX_fn;
}
}
return Bx;
}
......@@ -291,29 +287,28 @@ double MultipoleT::getBs(const Vector_t &R) {
if (angle_m == 0.0) {
// Straight geometry -> use corresponding field expansion
for(unsigned int n = 0; n <= maxOrder_m; n++) {
double f_n = 0.0;
for(unsigned int i = 0; i <= n; i++) {
f_n += gsl_sf_choose(n, i) * getTransDeriv(2 * i, R[0]) *
getFringeDeriv(2 * n - 2 * i + 1, R[2]);
double f_n = 0.0;
for(unsigned int i = 0; i <= n; i++) {
f_n += gsl_sf_choose(n, i) * getTransDeriv(2 * i, R[0]) *
getFringeDeriv(2 * n - 2 * i + 1, R[2]);
}
f_n *= gsl_sf_pow_int(-1.0, n);
Bs += gsl_sf_pow_int(R[1], 2 * n + 1) /
gsl_sf_fact(2 * n + 1) * f_n;
}
f_n *= gsl_sf_pow_int(-1.0, n);
Bs += gsl_sf_pow_int(R[1], 2 * n + 1) /
gsl_sf_fact(2 * n + 1) * f_n;
}
}
else {
} else {
if (variableRadius_m == true and getFringeDeriv(0, R[2]) < 1.0e-12) {
// return 0 if end of fringe field is reached
// this is to avoid functions being called at infinite radius
return 0.0;
}
}
// Curved geometry -> use corresponding field expansion
for(unsigned int n = 0; n <= maxOrder_m; n++) {
double partialS_fn = getFnDerivS(n, R[0], R[2]);
Bs += gsl_sf_pow_int(R[1], 2 * n + 1) / gsl_sf_fact(2 * n + 1);
Bs *= partialS_fn;
}
Bs /= getScaleFactor(R[0], R[2]);
double partialS_fn = getFnDerivS(n, R[0], R[2]);
Bs += gsl_sf_pow_int(R[1], 2 * n + 1) / gsl_sf_fact(2 * n + 1);
Bs *= partialS_fn;
}
Bs /= getScaleFactor(R[0], R[2]);
}
return Bs;
}
......@@ -336,23 +331,22 @@ double MultipoleT::getTransDeriv(unsigned int n, double x) {
std::vector<double> temp = transProfile_m;
if (n <= transMaxOrder_m) {
if (n != 0) {
for(unsigned int i = 1; i <= n; i++) {
for(unsigned int j = 0; j <= transMaxOrder_m; j++) {
if (j <= transMaxOrder_m - i) {
// move terms to the left and multiply by power
temp[j] = temp[j + 1] * (j + 1);
for(unsigned int i = 1; i <= n; i++) {
for(unsigned int j = 0; j <= transMaxOrder_m; j++) {
if (j <= transMaxOrder_m - i) {
// move terms to the left and multiply by power
temp[j] = temp[j + 1] * (j + 1);
} else {
// put 0 at the end for missing higher powers
temp[j] = 0.0;
}
}
else {
// put 0 at the end for missing higher powers
temp[j] = 0.0;
}
}
}
}
// Now use the vector to calculate value of the function
// Now use the vector to calculate value of the function
for(unsigned int k = 0; k <= transMaxOrder_m; k++) {
func += temp[k] * gsl_sf_pow_int(x, k);
}
}
}
return func;
}
......@@ -401,8 +395,7 @@ double MultipoleT::getRadius(double s) {
// move to current position on central axis
if (getFringeDeriv(0, s) != 0.0) {
return propCoeff / getFringeDeriv(0, s);
}
else {
} else {
return -1; // return -1 if radius is infinite
}
}
......@@ -433,12 +426,11 @@ double MultipoleT::getScaleFactor(double x, double s) {
if (not variableRadius_m) {
double rho = length_m / angle_m;
return 1 + x / rho;
}
else {
} else {
double temp = 0.0;
temp += gsl_sf_pow_int(getRadiusFirstDeriv(s), 2.0);
temp += gsl_sf_pow_int(1 + x / getRadius(s), 2.0);
return gsl_sf_pow_int(temp, 0.5);
temp += gsl_sf_pow_int(getRadiusFirstDeriv(s), 2.0);
temp += gsl_sf_pow_int(1 + x / getRadius(s), 2.0);
return gsl_sf_pow_int(temp, 0.5);
}
}
......@@ -546,24 +538,23 @@ double MultipoleT::getFn(unsigned int n, double x, double s) {
if (not variableRadius_m) {
double rho = length_m / angle_m;
double h_s = 1 + x / rho;
double temp = 0.0;
temp += getFnDerivX(n - 1, x, s) / (h_s * rho);
temp += getFnSecDerivX(n - 1, x, s);
temp += getFnSecDerivS(n - 1, x, s) / gsl_sf_pow_int(h_s, 2.0);
temp *= -1.0;
return temp;
}
// If radius is variable use corresponding recursion
else {
double temp = 0.0;
double h_s = getScaleFactor(x, s);
temp += getFnDerivX(n - 1, x, s) / (h_s * rho);
temp += getFnSecDerivX(n - 1, x, s);
temp += getFnSecDerivS(n - 1, x, s) / gsl_sf_pow_int(h_s, 2.0);
temp *= -1.0;
return temp;
} else {
// If radius is variable use corresponding recursion
double temp = 0.0;
double h_s = getScaleFactor(x, s);
temp += getScaleFactorDerivX(x, s) * getFnDerivX(n - 1, x, s);
temp -= getScaleFactorDerivS(x, s) * getFnDerivS(n - 1, x, s) /
gsl_sf_pow_int(h_s, 2.0);
temp += getFnSecDerivX(n - 1, x, s) * h_s;
temp += getFnSecDerivS(n - 1, x, s) / h_s;
temp *= -1.0 / h_s;
return temp;
temp -= getScaleFactorDerivS(x, s) * getFnDerivS(n - 1, x, s) /
gsl_sf_pow_int(h_s, 2.0);
temp += getFnSecDerivX(n - 1, x, s) * h_s;
temp += getFnSecDerivS(n - 1, x, s) / h_s;
temp *= -1.0 / h_s;
return temp;
}
}
......
......@@ -37,7 +37,7 @@
* ---------------------------------------------------------------------
*
* Class category: AbsBeamline \n
* $Author: titus
* $Author: Titus Dascalu, Chris Rogers
*
* ---------------------------------------------------------------------
*
......
......@@ -11,17 +11,17 @@ vector< vector<double> > partialsDerivB(const Vector_t &R,const Vector_t B, doub
double t = 0 ;
Vector_t P, E;
for(int i = 0; i < 3; i++)
{
// B at the previous and next grid points R_prev, R_next
Vector_t R_prev = R, R_next = R;
R_prev[i] -= stepSize;
R_next[i] += stepSize;
Vector_t B_prev, B_next;
dummyField->apply(R_prev, P, t, E, B_prev);
dummyField->apply(R_next, P, t, E, B_next);
for(int j = 0; j < 3; j++)
allPartials[i][j] = (B_next[j] - B_prev[j]) / (2 * stepSize);
}
{
// B at the previous and next grid points R_prev, R_next
Vector_t R_prev = R, R_next = R;
R_prev[i] -= stepSize;
R_next[i] += stepSize;
Vector_t B_prev, B_next;
dummyField->apply(R_prev, P, t, E, B_prev);
dummyField->apply(R_next, P, t, E, B_next);
for(int j = 0; j < 3; j++)
allPartials[i][j] = (B_next[j] - B_prev[j]) / (2 * stepSize);
}
return allPartials;
}
......@@ -32,22 +32,22 @@ vector< vector<double> > partialsDerivB_5(const Vector_t &R,const Vector_t B, do
double t = 0 ;
Vector_t P, E;
for(int i = 0; i < 3; i++)
{
// B at the previous and next grid points R_prev, R_next
Vector_t R_pprev = R, R_prev = R, R_next = R, R_nnext = R;
R_pprev(i) -= 2 * stepSize;
R_nnext(i) += 2 * stepSize;
R_prev(i) -= stepSize;
R_next(i) += stepSize;
Vector_t B_prev, B_next, B_pprev, B_nnext;
dummyField->apply(R_prev, P, t, E, B_prev);
dummyField->apply(R_next, P, t, E, B_next);
dummyField->apply(R_pprev, P, t, E, B_pprev);
dummyField->apply(R_nnext, P, t, E, B_nnext);
for(int j = 0; j < 3; j++)
allPartials[i][j] = (B_pprev[j] - 8 * B_prev[j] + 8 * B_next[j] - B_nnext[j]) / (12 * stepSize);
}
return allPartials;
{
// B at the previous and next grid points R_prev, R_next
Vector_t R_pprev = R, R_prev = R, R_next = R, R_nnext = R;
R_pprev(i) -= 2 * stepSize;
R_nnext(i) += 2 * stepSize;
R_prev(i) -= stepSize;
R_next(i) += stepSize;
Vector_t B_prev, B_next, B_pprev, B_nnext;
dummyField->apply(R_prev, P, t, E, B_prev);
dummyField->apply(R_next, P, t, E, B_next);
dummyField->apply(R_pprev, P, t, E, B_pprev);
dummyField->apply(R_nnext, P, t, E, B_nnext);
for(int j = 0; j < 3; j++)
allPartials[i][j] = (B_pprev[j] - 8 * B_prev[j] + 8 * B_next[j] - B_nnext[j]) / (12 * stepSize);
}
return allPartials;
}
double calcDivB(Vector_t &R, Vector_t B, double stepSize, MultipoleT* dummyField )
......@@ -71,7 +71,7 @@ vector<double> calcCurlB(Vector_t &R, Vector_t B, double stepSize, MultipoleT* d
return curl;
}
TEST(MultipoleT, Field)
TEST(MultipoleTTest, Field)
{
MultipoleT* myMagnet = new MultipoleT("Quadrupole");
double centralField = 5;
......@@ -82,29 +82,31 @@ TEST(MultipoleT, Field)
myMagnet->setTransMaxOrder(1);
myMagnet->setDipoleConstant(0.0);
myMagnet->setTransProfile(1, 100.0);
cout<<"test: "<<myMagnet->getTransProfile(1)<<' '<<myMagnet->getTransMaxOrder()<<endl;
cout << "test: " << myMagnet->getTransProfile(1) << ' '
<< myMagnet->getTransMaxOrder() << endl;
//highest power in the field is z ^ (2 * maxOrder + 1)
// !!! should be less than max_index / 2 !!!
myMagnet->setMaxOrder(3);
ofstream fout("3D_quad_");
Vector_t R(0., 0., 0.), P(3), E(3);
double t = 0., x, z, s;
for(x = - 0.02; x <= 0.02 ; x += 0.005)
for(z = 0.0; z <= 0.02 ; z += 0.005)
for(s = -10; s <= 10 ; s += 1)
{
R[0] = x;
R[1] = z;
R[2] = s;
Vector_t B(0., 0., 0.);
myMagnet->apply(R, P, t, E, B);
fout<<x<<' '<<z<<' '<<s<<' '<<B[0]<<' '<<B[1]<<' '<<B[2]<<endl;
}
for(x = - 0.02; x <= 0.02 ; x += 0.005) {
for(z = 0.0; z <= 0.02 ; z += 0.005) {
for(s = -10; s <= 10 ; s += 1) {
R[0] = x;
R[1] = z;
R[2] = s;
Vector_t B(0., 0., 0.);
myMagnet->apply(R, P, t, E, B);
fout << x << ' ' << z << ' ' << s << ' '
<< B[0] << ' ' << B[1] << ' ' << B[2] << endl;
}
}
}
fout.close();
}
TEST(MultipoleT, Maxwell)
{
TEST(MultipoleTTest, Maxwell) {
MultipoleT* myMagnet = new MultipoleT("Quadrupole");
double centralField = 5;
double fringeLength = 0.5;
......@@ -121,75 +123,37 @@ TEST(MultipoleT, Maxwell)
//ofstream fout("Quad_CurlB_off");
Vector_t R(0., 0., 0.), P(3), E(3);
double t = 0., x, z, s, stepSize= 1e-6;
for(x = 0.0; x <= 0.0 ; x += 0.001)
for(z = 0.0; z <= 0.02 ; z += 0.001)
for(s = -10; s <= 10 ; s += 0.1)
{
R[0] = x;
R[1] = z;
R[2] = s;
Vector_t B(0., 0., 0.);
myMagnet->apply(R, P, t, E, B);
double div = 0.;
div = calcDivB(R, B, stepSize, myMagnet);
EXPECT_NEAR(div, 0.0, 0.01);
//fout<<x<<' '<<z<<' '<<s<<' '<<B[0]<<' '<<B[1]<<' '<<B[2]<<' '<<div<<endl;
vector<double> curl;
curl = calcCurlB(R, B, stepSize, myMagnet);
EXPECT_NEAR(curl[0], 0.0, 1e-4);
EXPECT_NEAR(curl[1], 0.0, 1e-4);
EXPECT_NEAR(curl[2], 0.0, 1e-4);
//fout<<x<<' '<<z<<' '<<s<<' '<<B[0]<<' '<<B[1]<<' '<<B[2]<<' '<<
//sqrt(pow(curl[0], 2)+pow(curl[1],2)+pow(curl[2],2))<<endl;
}
//fout.close();
for(x = 0.0; x <= 0.0 ; x += 0.001) {
for(z = 0.0; z <= 0.02 ; z += 0.001) {
for(s = -10; s <= 10 ; s += 0.1) {
R[0] = x;
R[1] = z;
R[2] = s;
Vector_t B(0., 0., 0.);
myMagnet->apply(R, P, t, E, B);
double div = 0.;
div = calcDivB(R, B, stepSize, myMagnet);
EXPECT_NEAR(div, 0.0, 0.01);
vector<double> curl;
curl = calcCurlB(R, B, stepSize, myMagnet);
EXPECT_NEAR(curl[0], 0.0, 1e-4);
EXPECT_NEAR(curl[1], 0.0, 1e-4);
EXPECT_NEAR(curl[2], 0.0, 1e-4);
}
}
}
}
TEST(MultipoleT, Convergence_radius)
{
MultipoleT* myMagnet = new MultipoleT("Dipole");
double centralField = 1;
// the highest differential in the fringe field
double max_index = 10;
//Set the magnet
myMagnet->setTransMaxOrder(0);
myMagnet->setDipoleConstant(1.0);
/**for(int order = 1; order <= 10; order ++)
myMagnet->setTransProfile(order, 1.0);*/
//highest power in the field is z ^ (2 * maxOrder + 1)
// !!! should be less than max_index / 2 !!!
myMagnet->setMaxOrder(3);
ofstream fout("Convergence9");
Vector_t R(0., 0., 0.), P(3), E(3);
double t = 0., x, z, s;
double fringeLength = 0.5;
myMagnet->setFringeField(centralField, fringeLength, max_index);
for(x = 0.0; x <= 0.0 ; x += 0.001)
for(z = 0.0; z <= 0.7 ; z += 0.01)
for(s = -2; s <= 2 ; s += 0.001)
{
R[0] = x;
R[1] = z;
R[2] = s;
Vector_t B(0., 0., 0.);
myMagnet->apply(R, P, t, E, B);
fout<<x<<' '<<z<<' '<<s<<' '<<B[0]<<' '<<B[1]<<' '<<B[2]<<' '<<fringeLength<<endl;
}
fout.close();
}
TEST(MultipoleT, Curved_magnet)
{
TEST(MultipoleTTest, CurvedMagnet) {
MultipoleT* myMagnet = new MultipoleT("Combined function");
myMagnet->setLength(4.0);
myMagnet->setBendAngle(1.57);
myMagnet->setBendAngle(0.0); // BUG small, non-zero bend angle ruins the convergence
myMagnet->setAperture(0.4, 0.4);
myMagnet->setFringeField(2, 0.5, 0.5);
myMagnet->setVarRadius();
myMagnet->setVarStep(0.1);
myMagnet->setTransMaxOrder(1);