Skip to content

GitLab

Commit 34b8cafa authored by snuverink_j's avatar snuverink_j
Browse files

add B-field threshold for applying trim coils (trim coils are typically on the magnets only)

parent 786b1516
Hide whitespace changes
Inline Side-by-side
Showing with 99 additions and 77 deletions
src/Algorithms/ParallelCyclotronTracker.cpp
View file @ 34b8cafa
......@@ -573,7 +573,7 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {
} else if(type == std::string("BANDRF")) {
fieldflag = 6;
} else if(type == std::string("SYNCHROCYCLOTRON")) {
fieldflag = 7;
fieldflag = 7;
} else //(type == "RING")
fieldflag = 1;
......
src/Classic/AbsBeamline/Cyclotron.cpp
View file @ 34b8cafa
......@@ -136,6 +136,14 @@ double Cyclotron::getPZinit() const {
return pzinit_m;
}
void Cyclotron::setTrimCoilThreshold(double trimCoilThreshold) {
trimCoilThreshold_m = trimCoilThreshold;
}
double Cyclotron::getTrimCoilThreshold() const {
return trimCoilThreshold_m;
}
void Cyclotron::setSpiralFlag(bool spiral_flag) {
spiral_flag_m = spiral_flag;
}
......@@ -336,7 +344,7 @@ bool Cyclotron::apply(const Vector_t &R, const Vector_t &P, const double &t, Vec
/* if((R[0] > 0) && (R[1] >= 0)) tet = tempv;
else*/
if((R[0] < 0) && (R[1] >= 0)) tet = pi + tempv;
if ((R[0] < 0) && (R[1] >= 0)) tet = pi + tempv;
else if((R[0] < 0) && (R[1] <= 0)) tet = pi + tempv;
else if((R[0] > 0) && (R[1] <= 0)) tet = 2.0 * pi + tempv;
else if((R[0] == 0) && (R[1] > 0)) tet = pi / 2.0;
......@@ -443,112 +451,115 @@ bool Cyclotron::apply(const Vector_t &R, const Vector_t &P, const double &t, Vec
// bt = -( btf - btcub );
bt = - btf;
applyTrimCoil(rad, R[2], &br, &bz);
if (std::abs(br) > trimCoilThreshold_m || std::abs(bz) > trimCoilThreshold_m)
applyTrimCoil(rad, R[2], &br, &bz);
/* Br Btheta -> Bx By */
B[0] = br * cos(tet_rad) - bt * sin(tet_rad);
B[1] = br * sin(tet_rad) + bt * cos(tet_rad);
B[2] = bz;
//*gmsg << "R = " << rad << ", Theta = " << tet << ", B = (" << B[0] << "/" << B[1] << "/" << B[2] << ")" << endl;
//*gmsg << "R = " << rad << ", Theta = " << tet << ", B = (" << B[0] << "/" << B[1] << "/" << B[2] << ")" << endl;
} else {
return true;
}
if(myBFieldType_m == SYNCHRO || myBFieldType_m == BANDRF) {
//The RF field is supposed to be sampled on a cartesian grid
vector<Fieldmap *>::const_iterator fi = RFfields_m.begin();
vector<double>::const_iterator rffi = rffrequ_m.begin();
vector<double>::const_iterator rfphii = rfphi_m.begin();
vector<double>::const_iterator escali = escale_m.begin();
vector<bool>::const_iterator superposei = superpose_m.begin();
vector< vector<double> >::const_iterator rffci;
vector< vector<double> >::const_iterator rfvci;
if(myBFieldType_m != SYNCHRO && myBFieldType_m != BANDRF) {
return false;
}
//The RF field is supposed to be sampled on a cartesian grid
vector<Fieldmap *>::const_iterator fi = RFfields_m.begin();
vector<double>::const_iterator rffi = rffrequ_m.begin();
vector<double>::const_iterator rfphii = rfphi_m.begin();
vector<double>::const_iterator escali = escale_m.begin();
vector<bool>::const_iterator superposei = superpose_m.begin();
vector< vector<double> >::const_iterator rffci;
vector< vector<double> >::const_iterator rfvci;
if(myBFieldType_m == SYNCHRO) {
rffci = rffc_m.begin();
rfvci = rfvc_m.begin();
}
double xBegin(0), xEnd(0), yBegin(0), yEnd(0), zBegin(0), zEnd(0);
int fcount = 0;
for(; fi != RFfields_m.end(); ++fi, ++rffi, ++rfphii, ++escali, ++superposei) {
if(myBFieldType_m == SYNCHRO) {
rffci = rffc_m.begin();
rfvci = rfvc_m.begin();
++rffci, ++rfvci;
}
double xBegin(0), xEnd(0), yBegin(0), yEnd(0), zBegin(0), zEnd(0);
int fcount = 0;
for(; fi != RFfields_m.end(); ++fi, ++rffi, ++rfphii, ++escali, ++superposei) {
if(myBFieldType_m == SYNCHRO) {
++rffci, ++rfvci;
}
(*fi)->getFieldDimensions(xBegin, xEnd, yBegin, yEnd, zBegin, zEnd);
bool SuperPose = *superposei;
if (fcount > 0 && !SuperPose) {
//INFOMSG ("Field maps taken : " << fcount << "Superpose false" << endl);
break;
}
(*fi)->getFieldDimensions(xBegin, xEnd, yBegin, yEnd, zBegin, zEnd);
bool SuperPose = *superposei;
if (fcount > 0 && !SuperPose) {
//INFOMSG ("Field maps taken : " << fcount << "Superpose false" << endl);
break;
}
// Ok, this is a total patch job, but now that the internal cyclotron units are in m, we have to
// change stuff here to match with the input units of mm in the fieldmaps. -DW
const Vector_t temp_R = R * Vector_t(1000.0); //Keep this until we have transitioned fully to m -DW
// Ok, this is a total patch job, but now that the internal cyclotron units are in m, we have to
// change stuff here to match with the input units of mm in the fieldmaps. -DW
const Vector_t temp_R = R * Vector_t(1000.0); //Keep this until we have transitioned fully to m -DW
if ((temp_R(0) >= xBegin && temp_R(0) <= xEnd && temp_R(1) >= yBegin && temp_R(1) <= yEnd && temp_R(2) >= zBegin && temp_R(2) <= zEnd) == false)
continue;
if ((temp_R(0) >= xBegin && temp_R(0) <= xEnd && temp_R(1) >= yBegin && temp_R(1) <= yEnd && temp_R(2) >= zBegin && temp_R(2) <= zEnd) == false)
continue;
Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
if((*fi)->getFieldstrength(temp_R, tmpE, tmpB) == true) // out of bounds
continue;
Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
if((*fi)->getFieldstrength(temp_R, tmpE, tmpB) == true) // out of bounds
continue;
++fcount;
++fcount;
double frequency = (*rffi); // frequency in MHz
double ebscale = (*escali); // E and B field scaling
double frequency = (*rffi); // frequency in MHz
double ebscale = (*escali); // E and B field scaling
if(myBFieldType_m == SYNCHRO) {
double powert = 1;
for(const double fcoef : (*rffci)) {
powert *= (t * 1e-9);
frequency += fcoef * powert; // Add frequency ramp (in MHz/s^n)
}
if(myBFieldType_m == SYNCHRO) {
double powert = 1;
for(const double fcoef : (*rffci)) {
powert *= (t * 1e-9);
frequency += fcoef * powert; // Add frequency ramp (in MHz/s^n)
}
powert = 1;
for(const double vcoef : (*rfvci)) {
powert *= (t * 1e-9);
ebscale += vcoef * powert; // Add frequency ramp (in MHz/s^n)
}
powert = 1;
for(const double vcoef : (*rfvci)) {
powert *= (t * 1e-9);
ebscale += vcoef * powert; // Add frequency ramp (in MHz/s^n)
}
}
double phase = 2.0 * pi * 1.0E-3 * frequency * t + (*rfphii); // f in [MHz], t in [ns]
double phase = 2.0 * pi * 1.0E-3 * frequency * t + (*rfphii); // f in [MHz], t in [ns]
E += ebscale * cos(phase) * tmpE;
B -= ebscale * sin(phase) * tmpB;
E += ebscale * cos(phase) * tmpE;
B -= ebscale * sin(phase) * tmpB;
// INFOMSG("Field " << fcount << " BANDRF E= " << tmpE << " R= " << R << " phase " << phase << endl);
// INFOMSG("Field " << fcount << " BANDRF E= " << tmpE << " R= " << R << " phase " << phase << endl);
if(myBFieldType_m != SYNCHRO)
continue;
if(myBFieldType_m != SYNCHRO)
continue;
// Some phase output -DW
// Some phase output -DW
if (tet >= 90.0 && waiting_for_gap == 1) {
if (tet >= 90.0 && waiting_for_gap == 1) {
double phase_print = 180.0 * phase / pi;
phase_print = fmod(phase_print, 360) - 360.0;
double phase_print = 180.0 * phase / pi;
phase_print = fmod(phase_print, 360) - 360.0;
*gmsg << endl << "Gap 1 phase = " << phase_print << " Deg" << endl;
*gmsg << "Gap 1 E-Field = (" << E[0] << "/" << E[1] << "/" << E[2] << ")" << endl;
*gmsg << "Gap 1 B-Field = (" << B[0] << "/" << B[1] << "/" << B[2] << ")" << endl;
*gmsg << "RF Frequency = " << frequency << " MHz" << endl;
*gmsg << endl << "Gap 1 phase = " << phase_print << " Deg" << endl;
*gmsg << "Gap 1 E-Field = (" << E[0] << "/" << E[1] << "/" << E[2] << ")" << endl;
*gmsg << "Gap 1 B-Field = (" << B[0] << "/" << B[1] << "/" << B[2] << ")" << endl;
*gmsg << "RF Frequency = " << frequency << " MHz" << endl;
waiting_for_gap = 2;
}
else if (tet >= 270.0 && waiting_for_gap == 2) {
waiting_for_gap = 2;
}
else if (tet >= 270.0 && waiting_for_gap == 2) {
double phase_print = 180.0 * phase / pi;
phase_print = fmod(phase_print, 360) - 360.0;
double phase_print = 180.0 * phase / pi;
phase_print = fmod(phase_print, 360) - 360.0;
*gmsg << endl << "Gap 2 phase = " << phase_print << " Deg" << endl;
*gmsg << "Gap 2 E-Field = (" << E[0] << "/" << E[1] << "/" << E[2] << ")" << endl;
*gmsg << "Gap 2 B-Field = (" << B[0] << "/" << B[1] << "/" << B[2] << ")" << endl;
*gmsg << "RF Frequency = " << frequency << " MHz" << endl;
waiting_for_gap = 0;
}
*gmsg << endl << "Gap 2 phase = " << phase_print << " Deg" << endl;
*gmsg << "Gap 2 E-Field = (" << E[0] << "/" << E[1] << "/" << E[2] << ")" << endl;
*gmsg << "Gap 2 B-Field = (" << B[0] << "/" << B[1] << "/" << B[2] << ")" << endl;
*gmsg << "RF Frequency = " << frequency << " MHz" << endl;
waiting_for_gap = 0;
}
}
return false;
......
src/Classic/AbsBeamline/Cyclotron.h
View file @ 34b8cafa
......@@ -179,6 +179,9 @@ public:
void setFMHighE(double e);
virtual double getFMHighE() const;
void setTrimCoilThreshold(double);
virtual double getTrimCoilThreshold() const;
void setSpiralFlag(bool spiral_flag);
virtual bool getSpiralFlag() const;
......@@ -235,6 +238,7 @@ private:
double pzinit_m;
bool spiral_flag_m;
double trimCoilThreshold_m; ///< B-field threshold for applying trim coil
std::string type_m; /* what type of field we use */
double harm_m;
......
src/Elements/OpalCyclotron.cpp
View file @ 34b8cafa
......@@ -70,7 +70,7 @@ OpalCyclotron::OpalCyclotron():
("ESCALE", "Scale factor for the RF field(s)");
itsAttr[SUPERPOSE]= Attributes::makeBoolArray
("SUPERPOSE", "If TRUE, all of the electric field maps are superposed, only used when TYPE = BANDRF");
("SUPERPOSE", "If TRUE, all of the electric field maps are superposed, only used when TYPE = BANDRF");
itsAttr[RFMAPFN] = Attributes::makeStringArray
("RFMAPFN", "Filename(s) for the RF fieldmap(s)");
......@@ -110,7 +110,10 @@ OpalCyclotron::OpalCyclotron():
itsAttr[SPIRAL] = Attributes::makeBool
("SPIRAL","Flag whether or not this is a spiral inflector simulation", false);
itsAttr[TRIMCOILTHRESHOLD] = Attributes::makeReal
("TRIMCOILTHRESHOLD","Minimum magnetic field [T] for which trim coils are applied", 0.0);
itsAttr[TRIMCOIL] = Attributes::makeStringArray
("TRIMCOIL", "List of trim coils");
......@@ -138,6 +141,7 @@ OpalCyclotron::OpalCyclotron():
registerRealAttribute("MAXR");
registerRealAttribute("FMLOWE");
registerRealAttribute("FMHIGHE");
registerRealAttribute("TRIMCOILTHRESHOLD");
registerOwnership();
......@@ -192,6 +196,7 @@ void OpalCyclotron::update() {
double fmHighE = Attributes::getReal(itsAttr[FMHIGHE]);
bool spiral_flag = Attributes::getBool(itsAttr[SPIRAL]);
double trimCoilThreshold = Attributes::getReal(itsAttr[TRIMCOILTHRESHOLD]);
cycl->setFieldMapFN(fmapfm);
cycl->setSymmetry(symmetry);
......@@ -216,6 +221,7 @@ void OpalCyclotron::update() {
cycl->setFMHighE(fmHighE);
cycl->setSpiralFlag(spiral_flag);
cycl->setTrimCoilThreshold(trimCoilThreshold);
std::vector<std::string> fm_str = Attributes::getStringArray(itsAttr[RFMAPFN]);
std::vector<std::string> rffcfn_str = Attributes::getStringArray(itsAttr[RFFCFN]);
......
src/Elements/OpalCyclotron.h
View file @ 34b8cafa
......@@ -59,6 +59,7 @@ public:
FMLOWE, // minimal energy of the field map
FMHIGHE, // maximal energy of the field map
SPIRAL, // flag whether or not this is a spiral inflector simulation
TRIMCOILTHRESHOLD, // minimum B-field for which trim coils are applied
TRIMCOIL, // list of trim coils
SIZE
};
......
Markdown is supported
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