Commit 5d6c9d35 authored by Steve Russell's avatar Steve Russell
Browse files

Fixes to SBend and RBend. User guide updates to come (about 1/2 done).

parent 67d65ee6
......@@ -43,6 +43,8 @@ RBend::RBend():
designEnergy_m(0.0),
designRadius_m(0.0),
fieldAmplitude_m(0.0),
bX_m(0.0),
bY_m(0.0),
entranceAngle_m(0.0),
exitAngle_m(0.0),
gradient_m(0.0),
......@@ -95,6 +97,8 @@ RBend::RBend(const RBend &right):
designEnergy_m(right.designEnergy_m),
designRadius_m(right.designRadius_m),
fieldAmplitude_m(right.fieldAmplitude_m),
bX_m(right.bX_m),
bY_m(right.bY_m),
entranceAngle_m(right.entranceAngle_m),
exitAngle_m(right.exitAngle_m),
gradient_m(right.gradient_m),
......@@ -149,6 +153,8 @@ RBend::RBend(const std::string &name):
designEnergy_m(0.0),
designRadius_m(0.0),
fieldAmplitude_m(0.0),
bX_m(0.0),
bY_m(0.0),
entranceAngle_m(0.0),
exitAngle_m(0.0),
gradient_m(0.0),
......@@ -478,8 +484,9 @@ void RBend::SetEntranceAngle(double entranceAngle) {
entranceAngle_m = entranceAngle;
}
void RBend::SetFieldAmplitude(double fieldAmplitude) {
fieldAmplitude_m = fieldAmplitude;
void RBend::SetFieldAmplitude(double k0, double k0s) {
bY_m = k0;
bX_m = k0s;
}
void RBend::SetFieldMapFN(std::string fileName) {
......@@ -502,13 +509,6 @@ void RBend::SetRotationAboutZ(double rotation) {
Orientation_m(2) = rotation;
}
void RBend::SetRotationAboutZ(double k0, double k0s) {
Orientation_m(2) = atan2(k0s, k0);
fieldAmplitude_m = sqrt(pow(k0, 2.0) + pow(k0s, 2.0));
}
void RBend::AdjustFringeFields(double ratio) {
double delta = std::abs(entranceParameter1_m - entranceParameter2_m);
......@@ -1070,16 +1070,16 @@ bool RBend::FindIdealBendParameters(double bendLength) {
return true;
} else {
} else if(bX_m == 0.0) {
// Negative angle is a positive bend rotated 180 degrees.
if((refCharge > 0.0 && fieldAmplitude_m < 0.0)
|| (refCharge < 0.0 && fieldAmplitude_m > 0.0)) {
if((refCharge > 0.0 && bY_m < 0.0)
|| (refCharge < 0.0 && bY_m > 0.0)) {
gradient_m *= -1.0;
Orientation_m(2) += Physics::pi;
}
fieldAmplitude_m = std::abs(fieldAmplitude_m);
fieldAmplitude_m = refCharge * std::abs(bY_m / refCharge);
designRadius_m = std::abs(refBetaGamma * refMass / (Physics::c * fieldAmplitude_m));
double angle = asin(bendLength / designRadius_m - sin(entranceAngle_m));
angle_m = angle + entranceAngle_m;
......@@ -1087,6 +1087,23 @@ bool RBend::FindIdealBendParameters(double bendLength) {
return false;
} else {
Orientation_m(2) += atan2(bX_m, bY_m);
if(refCharge < 0.0) {
gradient_m *= -1.0;
Orientation_m(2) -= Physics::pi;
}
fieldAmplitude_m = refCharge
* std::abs(sqrt(pow(bY_m, 2.0) + pow(bX_m, 2.0))
/ refCharge);
designRadius_m = std::abs(refBetaGamma * refMass / (Physics::c * fieldAmplitude_m));
double angle = asin(bendLength / designRadius_m - sin(entranceAngle_m));
angle_m = angle + entranceAngle_m;
exitAngle_m = angle_m - entranceAngle_m;
return false;
}
}
......@@ -1228,9 +1245,9 @@ void RBend::Print(Inform &msg, double bendAngleX, double bendAngleY) {
<< " m (in s coordinates)"
<< endl;
msg << endl
<< "Bend Magnet Properties"
<< "Ideal (Hard Edge) Bend Magnet Properties"
<< endl
<< "======================"
<< "========================================"
<< endl << endl;
msg << "Bend angle magnitude: "
<< angle_m
......@@ -1262,6 +1279,10 @@ void RBend::Print(Inform &msg, double bendAngleX, double bendAngleY) {
<< designRadius_m
<< " m"
<< endl;
msg << "Bend design energy: "
<< designEnergy_m
<< " eV"
<< endl;
msg << "Rotation about x axis: "
<< Orientation_m(1)
<< " rad ("
......@@ -1281,9 +1302,9 @@ void RBend::Print(Inform &msg, double bendAngleX, double bendAngleY) {
<< " degrees)"
<< endl;
msg << endl
<< "Reference Trajectory Properties"
<< "Reference Trajectory Properties Through Bend Magnet with Fringe Fields"
<< endl
<< "==============================="
<< "======================================================================"
<< endl << endl;
msg << "Reference particle is bent: "
<< bendAngleX
......@@ -1599,6 +1620,15 @@ bool RBend::TreatAsDrift(Inform &msg) {
<< endl;
designRadius_m = 0.0;
return true;
} else if(angle_m == 0.0 &&
pow(bX_m, 2.0) + pow(bY_m, 2.0) == 0.0) {
msg << "Warning bend angle/strength is zero. Treating as drift."
<< endl;
designRadius_m = 0.0;
return true;
} else
return false;
}
......@@ -215,7 +215,7 @@ public:
void SetBeta(double beta);
void SetDesignEnergy(double energy);
void SetEntranceAngle(double entranceAngle);
void SetFieldAmplitude(double fieldAmplitude);
void SetFieldAmplitude(double k0, double k0s);
/*
* Set the name of the field map.
......@@ -232,7 +232,6 @@ public:
/// Set rotation about z axis in bend frame.
void SetRotationAboutZ(double rotation);
void SetRotationAboutZ(double k0, double k0s);
private:
......@@ -303,11 +302,12 @@ private:
/// (Not currently used.)
double angle_m; /// Bend angle for reference particle with bend
/// design energy (radians).
double aperture_m; /// Aperture of magnet in non-bend (horizontal)
/// plane.
double aperture_m; /// Aperture of magnet in non-bend plane.
double designEnergy_m; /// Bend design energy (eV).
double designRadius_m; /// Bend design radius (m).
double fieldAmplitude_m; /// Amplitude of magnet field (T).
double bX_m; /// Amplitude of Bx field (T).
double bY_m; /// Amplitude of By field (T).
double entranceAngle_m; /// Angle between incoming reference trajectory
/// and the entrance face of the magnet (radians).
double exitAngle_m; /// Angle between outgoing reference trajectory
......
......@@ -40,6 +40,8 @@ SBend::SBend():
designEnergy_m(0.0),
designRadius_m(0.0),
fieldAmplitude_m(0.0),
bX_m(0.0),
bY_m(0.0),
angleGreaterThanPi_m(false),
entranceAngle_m(0.0),
exitAngle_m(0.0),
......@@ -90,6 +92,8 @@ SBend::SBend(const SBend &right):
designEnergy_m(right.designEnergy_m),
designRadius_m(right.designRadius_m),
fieldAmplitude_m(right.fieldAmplitude_m),
bX_m(right.bX_m),
bY_m(right.bY_m),
angleGreaterThanPi_m(right.angleGreaterThanPi_m),
entranceAngle_m(right.entranceAngle_m),
exitAngle_m(right.exitAngle_m),
......@@ -145,6 +149,8 @@ SBend::SBend(const std::string &name):
designEnergy_m(0.0),
designRadius_m(0.0),
fieldAmplitude_m(0.0),
bX_m(0.0),
bY_m(0.0),
angleGreaterThanPi_m(false),
entranceAngle_m(0.0),
exitAngle_m(0.0),
......@@ -443,8 +449,9 @@ void SBend::setExitAngle(double exitAngle) {
exitAngle_m = exitAngle;
}
void SBend::SetFieldAmplitude(double fieldAmplitude) {
fieldAmplitude_m = fieldAmplitude;
void SBend::SetFieldAmplitude(double k0, double k0s) {
bY_m = k0;
bX_m = k0s;
}
void SBend::SetFieldMapFN(std::string fileName) {
......@@ -467,13 +474,6 @@ void SBend::SetRotationAboutZ(double rotation) {
Orientation_m(2) = rotation;
}
void SBend::SetRotationAboutZ(double k0, double k0s) {
Orientation_m(2) = atan2(k0s, k0);
fieldAmplitude_m = sqrt(pow(k0, 2.0) + pow(k0s, 2.0));
}
void SBend::AdjustFringeFields(double ratio) {
double delta = std::abs(entranceParameter1_m - entranceParameter2_m);
......@@ -1034,16 +1034,16 @@ bool SBend::FindIdealBendParameters(double chordLength) {
/ (Physics::c * designRadius_m);
return true;
} else {
} else if(bX_m == 0.0) {
// Negative angle is a positive bend rotated 180 degrees.
if((refCharge > 0.0 && fieldAmplitude_m < 0.0)
|| (refCharge < 0.0 && fieldAmplitude_m > 0.0)) {
if((refCharge > 0.0 && bY_m < 0.0)
|| (refCharge < 0.0 && bY_m > 0.0)) {
gradient_m *= -1.0;
Orientation_m(2) += Physics::pi;
}
fieldAmplitude_m = std::abs(fieldAmplitude_m);
fieldAmplitude_m = refCharge * std::abs(bY_m / refCharge);
designRadius_m = std::abs(refBetaGamma * refMass / (Physics::c * fieldAmplitude_m));
double bendAngle = 2.0 * std::asin(chordLength / (2.0 * designRadius_m));
......@@ -1053,6 +1053,26 @@ bool SBend::FindIdealBendParameters(double chordLength) {
angle_m = bendAngle;
return false;
} else {
Orientation_m(2) += atan2(bX_m, bY_m);
if(refCharge < 0.0) {
gradient_m *= -1.0;
Orientation_m(2) -= Physics::pi;
}
fieldAmplitude_m = refCharge
* std::abs(sqrt(pow(bY_m, 2.0) + pow(bX_m, 2.0))
/ refCharge);
designRadius_m = std::abs(refBetaGamma * refMass / (Physics::c * fieldAmplitude_m));
double bendAngle = 2.0 * std::asin(chordLength / (2.0 * designRadius_m));
if(angleGreaterThanPi_m)
bendAngle = 2.0 * Physics::pi - bendAngle;
angle_m = bendAngle;
return false;
}
}
......@@ -1194,9 +1214,9 @@ void SBend::Print(Inform &msg, double bendAngleX, double bendAngleY) {
<< " m (in s coordinates)"
<< endl;
msg << endl
<< "Bend Magnet Properties"
<< "Ideal (Hard Edge) Bend Magnet Properties"
<< endl
<< "======================"
<< "========================================"
<< endl << endl;
msg << "Bend angle magnitude: "
<< angle_m
......@@ -1228,6 +1248,10 @@ void SBend::Print(Inform &msg, double bendAngleX, double bendAngleY) {
<< designRadius_m
<< " m"
<< endl;
msg << "Bend design energy: "
<< designEnergy_m
<< " eV"
<< endl;
msg << "Rotation about x axis: "
<< Orientation_m(1)
<< " rad ("
......@@ -1247,9 +1271,9 @@ void SBend::Print(Inform &msg, double bendAngleX, double bendAngleY) {
<< " degrees)"
<< endl;
msg << endl
<< "Reference Trajectory Properties"
<< "Reference Trajectory Properties Through Bend Magnet with Fringe Fields"
<< endl
<< "==============================="
<< "======================================================================"
<< endl << endl;
msg << "Reference particle is bent: "
<< bendAngleX
......@@ -1484,7 +1508,7 @@ bool SBend::SetupBendGeometry(Inform &msg, double &startField, double &endField)
if(!FindChordLength(msg, chordLength, chordLengthFromMap))
return false;
if(TreatAsDrift(msg)) {
if(TreatAsDrift(msg, chordLength)) {
startField_m = startField;
endField_m = startField + chordLength;
return true;
......@@ -1558,12 +1582,39 @@ void SBend::SetupPusher(PartBunch *bunch) {
}
bool SBend::TreatAsDrift(Inform &msg) {
bool SBend::TreatAsDrift(Inform &msg, double chordLength) {
if(designEnergy_m <= 0.0) {
msg << "Warning: bend design energy is zero. Treating as drift."
<< endl;
designRadius_m = 0.0;
return true;
} else if(angle_m == 0.0) {
double refMass = RefPartBunch_m->getM();
double refGamma = designEnergy_m / refMass + 1.0;
double refBetaGamma = sqrt(pow(refGamma, 2.0) - 1.0);
double amplitude = std::abs(fieldAmplitude_m);
double radius = std::abs(refBetaGamma * refMass / (Physics::c * amplitude));
double sinArgument = chordLength / (2.0 * radius);
if(std::abs(sinArgument) > 1.0) {
msg << "Warning: bend strength and defined reference trajectory "
<< "chord length are not consistent. Treating bend as drift."
<< endl;
designRadius_m = 0.0;
return true;
} else
return false;
} else if(angle_m == 0.0 &&
pow(bY_m, 2.0) + pow(bX_m, 2.0) == 0.0) {
msg << "Warning bend angle/strength is zero. Treating as drift."
<< endl;
designRadius_m = 0.0;
return true;
} else
return false;
}
......
......@@ -206,7 +206,7 @@ public:
void SetDesignEnergy(double energy);
void SetEntranceAngle(double entranceAngle);
void setExitAngle(double exitAngle);
void SetFieldAmplitude(double fieldAmplitude);
void SetFieldAmplitude(double k0, double k0s);
/*
* Set the name of the field map.
......@@ -223,7 +223,6 @@ public:
/// Set rotation about z axis in bend frame.
void SetRotationAboutZ(double rotation);
void SetRotationAboutZ(double k0, double k0s);
private:
......@@ -284,7 +283,7 @@ private:
bool SetupDefaultFieldMap(Inform &msg);
void SetFieldBoundaries(double startField, double endField);
void SetupPusher(PartBunch *bunch);
bool TreatAsDrift(Inform &msg);
bool TreatAsDrift(Inform &msg, double chordlength);
BorisPusher pusher_m; /// Pusher used to integrate reference particle
/// through the bend.
......@@ -294,11 +293,12 @@ private:
/// (Not currently used.)
double angle_m; /// Bend angle for reference particle with bend
/// design energy (radians).
double aperture_m; /// Aperture of magnet in non-bend (horizontal)
/// plane.
double aperture_m; /// Aperture of magnet in non-bend plane.
double designEnergy_m; /// Bend design energy (eV).
double designRadius_m; /// Bend design radius (m).
double fieldAmplitude_m; /// Amplitude of magnet field (T).
double bX_m; /// Amplitude of Bx field (T).
double bY_m; /// Amplitude of By field (T).
bool angleGreaterThanPi_m; /// Set to true if bend angle is greater than
/// 180 degrees.
double entranceAngle_m; /// Angle between incoming reference trajectory
......
......@@ -67,8 +67,8 @@ OpalBend::OpalBend(const char *name, const char *help):
("FMAPFN", "Filename for the fieldmap");
itsAttr[GAP] = Attributes::makeReal
("GAP", "Full gap height of the magnet (m)", 0.0);
itsAttr[APERTURE] = Attributes::makeReal
("APERTURE", "Non-bend plane magnet aperture (m)", 0.0);
itsAttr[HAPERT] = Attributes::makeReal
("HAPERT", "Non-bend plane magnet aperture (m)", 0.0);
itsAttr[ROTATION] = Attributes::makeReal
("ROTATION", "Magnet rotation about z axis in degrees");
itsAttr[ALPHA] = Attributes::makeReal
......@@ -84,9 +84,9 @@ OpalBend::OpalBend(const char *name, const char *help):
"Set to true if bend angle is greater than 180 degrees",
false);
itsAttr[DX] = Attributes::makeReal
("DX", "Misalignment in x direction",0.0);
("DX", "Misalignment in x direction", 0.0);
itsAttr[DY] = Attributes::makeReal
("DY", "Misalignment in y direction",0.0);
("DY", "Misalignment in y direction", 0.0);
registerRealAttribute("ANGLE");
registerRealAttribute("K0L");
......@@ -107,7 +107,7 @@ OpalBend::OpalBend(const char *name, const char *help):
registerRealAttribute("STEPSIZE");
registerStringAttribute("FMAPFN");
registerRealAttribute("GAP");
registerRealAttribute("APERTURE");
registerRealAttribute("HAPERT");
registerRealAttribute("ROTATION");
registerRealAttribute("ALPHA");
registerRealAttribute("BETA");
......
......@@ -44,7 +44,7 @@ public:
SLICES, STEPSIZE, // Parameters used to determine slicing.
FMAPFN, // File name containing on-axis field.
GAP, // Full gap of magnet.
APERTURE, // Horizontal (non-bend plane) magnet aperture.
HAPERT, // Horizontal aperture of magnet.
ROTATION, // Magnet rotation about z axis.
ALPHA, // The edge angle 1
BETA, // The edge angle 2
......@@ -52,7 +52,7 @@ public:
EXITANGLE, // the relative angle between the entry and the exit face
GREATERTHANPI, // Boolean flag set to true if bend angle is greater
// than 180 degrees.
DX, // Misalignment: translation in x direction
DX, // Misalignment: translation in x direction
DY, // Misalignment: translation in y direction
SIZE // Total number of attributes.
};
......
......@@ -150,23 +150,26 @@ void OpalRBend::update() {
field.setSkewComponent(4, factor * Attributes::getReal(itsAttr[K3S]) / 6.0);
bend->setField(field);
// Set field amplitude or bend angle and the magnet rotation about the z axis.
// Set field amplitude or bend angle.
if(itsAttr[ANGLE])
bend->SetBendAngle(Attributes::getReal(itsAttr[ANGLE]));
else
bend->SetFieldAmplitude(sqrt(k0 * k0 + k0s * k0s));
bend->SetFieldAmplitude(k0, k0s);
if(itsAttr[ROTATION])
bend->SetRotationAboutZ(Attributes::getReal(itsAttr[ROTATION]));
else if(itsAttr[K0] || itsAttr[K0S])
bend->SetRotationAboutZ(k0, k0s);
else
bend->SetRotationAboutZ(0.0);
if(itsAttr[FMAPFN])
bend->SetFieldMapFN(Attributes::getString(itsAttr[FMAPFN]));
else if(bend->getName() != "RBEND")
ERRORMSG(bend->getName() << ": No filename for a field map given" << endl);
else if(bend->getName() != "RBEND") {
ERRORMSG(bend->getName() << ": No filename for a field map given. "
"Will assume the default map "
"\"1DPROFILE1-DEFAULT\"."
<< endl);
bend->SetFieldMapFN("1DPROFILE1-DEFAULT");
}
if(itsAttr[E1])
bend->SetEntranceAngle(Attributes::getReal(itsAttr[E1]));
......@@ -180,9 +183,9 @@ void OpalRBend::update() {
else
bend->SetBeta(0.0);
// Convert from MeV to eV.
// Energy in eV.
if(itsAttr[DESIGNENERGY]) {
bend->SetDesignEnergy(Attributes::getReal(itsAttr[DESIGNENERGY]) * 1e6);
bend->SetDesignEnergy(Attributes::getReal(itsAttr[DESIGNENERGY]));
}
if(itsAttr[GAP])
......@@ -190,8 +193,8 @@ void OpalRBend::update() {
else
bend->SetFullGap(0.0);
if(itsAttr[APERTURE])
bend->SetAperture(Attributes::getReal(itsAttr[APERTURE]));
if(itsAttr[HAPERT])
bend->SetAperture(Attributes::getReal(itsAttr[HAPERT]));
else
bend->SetAperture(0.0);
......@@ -211,6 +214,10 @@ void OpalRBend::update() {
else
bend->SetK1(0.0);
bend->setMisalignment(Attributes::getReal(itsAttr[DX]),
Attributes::getReal(itsAttr[DY]),
0.0);
std::vector<double> apert = getApert();
double apert_major = -1., apert_minor = -1.;
if(apert.size() > 0) {
......
......@@ -155,11 +155,11 @@ void OpalSBend::update() {
field.setSkewComponent(4, factor * Attributes::getReal(itsAttr[K3S]) / 6.0);
bend->setField(field);
// Set field amplitude or bend angle and the magnet rotation about the z axis.
// Set field amplitude or bend angle.
if(itsAttr[ANGLE])
bend->SetBendAngle(Attributes::getReal(itsAttr[ANGLE]));
else
bend->SetFieldAmplitude(sqrt(k0 * k0 + k0s * k0s));
bend->SetFieldAmplitude(k0, k0s);
if(itsAttr[GREATERTHANPI])
bend->SetAngleGreaterThanPiFlag(Attributes::getBool(itsAttr[GREATERTHANPI]));
......@@ -168,15 +168,18 @@ void OpalSBend::update() {
if(itsAttr[ROTATION])
bend->SetRotationAboutZ(Attributes::getReal(itsAttr[ROTATION]));
else if(itsAttr[K0] || itsAttr[K0S])
bend->SetRotationAboutZ(k0, k0s);
else
bend->SetRotationAboutZ(0.0);
if(itsAttr[FMAPFN])
bend->SetFieldMapFN(Attributes::getString(itsAttr[FMAPFN]));
else if(bend->getName() != "SBEND")
ERRORMSG(bend->getName() << ": No filename for a field map given" << endl);
else if(bend->getName() != "SBEND") {
ERRORMSG(bend->getName() << ": No filename for a field map given. "
"Will assume the default map "
"\"1DPROFILE1-DEFAULT\"."
<< endl);
bend->SetFieldMapFN("1DPROFILE1-DEFAULT");
}
if(itsAttr[E1])
bend->SetEntranceAngle(Attributes::getReal(itsAttr[E1]));
......@@ -190,9 +193,9 @@ void OpalSBend::update() {
else
bend->SetBeta(0.0);
// Convert from MeV to eV.
// Units are eV.
if(itsAttr[DESIGNENERGY]) {
bend->SetDesignEnergy(Attributes::getReal(itsAttr[DESIGNENERGY]) * 1e6);
bend->SetDesignEnergy(Attributes::getReal(itsAttr[DESIGNENERGY]));
}
if(itsAttr[GAP])
......@@ -200,8 +203,8 @@ void OpalSBend::update() {
else
bend->SetFullGap(0.0);
if(itsAttr[APERTURE])
bend->SetAperture(Attributes::getReal(itsAttr[APERTURE]));
if(itsAttr[HAPERT])
bend->SetAperture(Attributes::getReal(itsAttr[HAPERT]));
else
bend->SetAperture(0.0);
......@@ -228,6 +231,11 @@ void OpalSBend::update() {
bend->SetK1(Attributes::getReal(itsAttr[K1]));
else
bend->SetK1(0.0);
bend->setMisalignment(Attributes::getReal(itsAttr[DX]),
Attributes::getReal(itsAttr[DY]),
0.0);
// if (itsAttr[L])
// bend->setL(Attributes::getReal(itsAttr[L]));
// else
......
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