diff --git a/src/Algorithms/MultiBunchHandler.cpp b/src/Algorithms/MultiBunchHandler.cpp
index 45aed13e4274efff01066eebec9900a0714622f8..6c482a455ce7068d423aa2162d99950420d783d7 100644
--- a/src/Algorithms/MultiBunchHandler.cpp
+++ b/src/Algorithms/MultiBunchHandler.cpp
@@ -271,9 +271,7 @@ short MultiBunchHandler::injectBunch(PartBunchBase<double, 3>* beam,
*gmsg << "* MBM: Checking for automatically injecting new bunch ..." << endl;
- //beam->R *= Vector_t(0.001); // mm --> m
beam->calcBeamParameters();
- //beam->R *= Vector_t(1000.0); // m --> mm
Vector_t Rmean = beam->get_centroid(); // m
@@ -418,7 +416,6 @@ void MultiBunchHandler::setRadiusTurns(const double& radius) {
} else {
*gmsg << "Initial radial position = ";
}
- // New OPAL 2.0: Init in m -DW
*gmsg << radiusThisTurn_m << " m" << endl;
}
@@ -513,7 +510,7 @@ bool MultiBunchHandler::calcBunchBeamParameters(PartBunchBase<double, 3>* beam,
double invN = 1.0 / double(bunchTotalNum);
binfo.ekin = local[0] * invN;
- binfo.time = beam->getT() * Units::s2ns;
+ binfo.time = beam->getT();
binfo.nParticles = bunchTotalNum;
for (unsigned int i = 0; i < dim; ++i) {
diff --git a/src/Algorithms/MultiBunchHandler.h b/src/Algorithms/MultiBunchHandler.h
index 783b281971e65dfb257258224177e885f33deb3a..563af851d9763d54d9628081e622c1b852c371ce 100644
--- a/src/Algorithms/MultiBunchHandler.h
+++ b/src/Algorithms/MultiBunchHandler.h
@@ -42,10 +42,10 @@ public:
, radius(0.0)
{ };
- double time; // ns
+ double time; // s
double pathlength; // m
double azimuth; // deg
- double radius; // mm
+ double radius; // m
};
struct beaminfo_t {
diff --git a/src/Algorithms/ParallelCyclotronTracker.cpp b/src/Algorithms/ParallelCyclotronTracker.cpp
index 54f214a1d7253798fc6cdec062d663b30ef8a63e..253773b85225819c5e1d9d64d816c9cc1fc7dd94 100644
--- a/src/Algorithms/ParallelCyclotronTracker.cpp
+++ b/src/Algorithms/ParallelCyclotronTracker.cpp
@@ -26,7 +26,6 @@
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#include "Algorithms/ParallelCyclotronTracker.h"
-#include "Algorithms/BoostMatrix.h"
#include "AbsBeamline/CCollimator.h"
#include "AbsBeamline/Corrector.h"
@@ -61,6 +60,7 @@
#include "AbstractObjects/Element.h"
#include "AbstractObjects/OpalData.h"
+#include "Algorithms/BoostMatrix.h"
#include "Algorithms/Ctunes.h"
#include "Algorithms/PolynomialTimeDependence.h"
@@ -92,13 +92,13 @@
#include <limits>
#include <numeric>
-constexpr double c_mmtns = Physics::c * Units::m2mm / Units::s2ns;
-
Vector_t const ParallelCyclotronTracker::xaxis = Vector_t(1.0, 0.0, 0.0);
Vector_t const ParallelCyclotronTracker::yaxis = Vector_t(0.0, 1.0, 0.0);
Vector_t const ParallelCyclotronTracker::zaxis = Vector_t(0.0, 0.0, 1.0);
extern Inform *gmsg;
+extern Inform *gmsgALL;
+
/**
* Constructor ParallelCyclotronTracker
@@ -199,9 +199,8 @@ void ParallelCyclotronTracker::bgf_main_collision_test() {
itsBunch_m->Q[i], itsBunch_m->M[i]),
std::make_pair(turnnumber_m, itsBunch_m->bunchNum[i]));
itsBunch_m->Bin[i] = -1;
- Inform gmsgALL("OPAL", INFORM_ALL_NODES);
- gmsgALL << level4 << "* Particle " << itsBunch_m->ID[i]
- << " lost on boundary geometry" << endl;
+ *gmsgALL << level4 << "* Particle " << itsBunch_m->ID[i]
+ << " lost on boundary geometry" << endl;
}
}
}
@@ -241,7 +240,7 @@ void ParallelCyclotronTracker::computePathLengthUpdate(std::vector<double>& dl,
{
// the last element in dotP is the dot-product over all particles
std::vector<double> dotP(dl.size());
- if ( Options::psDumpFrame == DumpFrame::BUNCH_MEAN || isMultiBunch()) {
+ if ( Options::psDumpFrame == DumpFrame::BUNCH_MEAN || isMultiBunch() ) {
for (unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) {
dotP[itsBunch_m->bunchNum[i]] += dot(itsBunch_m->P[i], itsBunch_m->P[i]);
@@ -269,8 +268,7 @@ void ParallelCyclotronTracker::computePathLengthUpdate(std::vector<double>& dl,
for (size_t i = 0; i < dotP.size(); ++i) {
double const gamma = std::sqrt(1.0 + dotP[i]);
double const beta = std::sqrt(dotP[i]) / gamma;
- dl[i] = c_mmtns * dt * Units::mm2m * beta;
-
+ dl[i] = Physics::c * dt * beta;
}
}
@@ -349,14 +347,12 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron& cycl) {
if (!OpalData::getInstance()->inRestartRun()) {
// Get reference values from cyclotron element
- // For now, these are still stored in mm. should be the only ones. -DW
referenceR = cycl_m->getRinit();
referenceTheta = cycl_m->getPHIinit();
referenceZ = cycl_m->getZinit();
referencePr = cycl_m->getPRinit();
referencePz = cycl_m->getPZinit();
-
- referencePtot = itsReference.getGamma() * itsReference.getBeta();
+ referencePtot = itsReference.getGamma() * itsReference.getBeta();
// Calculate reference azimuthal (tangential) momentum from total-, z- and radial momentum:
float insqrt = referencePtot * referencePtot - \
@@ -398,29 +394,32 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron& cycl) {
}
// Adjust some of the reference variables from the h5 file
- referencePhi *= Units::deg2rad;
- referencePsi *= Units::deg2rad;
- referencePtot = bega;
+ referenceR *= Units::mm2m;
+ referenceZ *= Units::mm2m;
+ referenceTheta *= Units::deg2rad;
+ referencePhi *= Units::deg2rad;
+ referencePsi *= Units::deg2rad;
+ referencePtot = bega;
}
cycl_m->checkInitialReferenceParticle(referenceR, referenceTheta, referenceZ);
- sinRefTheta_m = std::sin(referenceTheta * Units::deg2rad);
- cosRefTheta_m = std::cos(referenceTheta * Units::deg2rad);
+ sinRefTheta_m = std::sin(referenceTheta);
+ cosRefTheta_m = std::cos(referenceTheta);
*gmsg << endl;
*gmsg << "* Bunch global starting position:" << endl;
- *gmsg << "* RINIT = " << referenceR << " [mm]" << endl;
- *gmsg << "* PHIINIT = " << referenceTheta << " [deg]" << endl;
- *gmsg << "* ZINIT = " << referenceZ << " [mm]" << endl;
+ *gmsg << "* RINIT = " << referenceR << " [m]" << endl;
+ *gmsg << "* PHIINIT = " << referenceTheta * Units::rad2deg << " [deg]" << endl;
+ *gmsg << "* ZINIT = " << referenceZ << " [m]" << endl;
*gmsg << endl;
*gmsg << "* Bunch global starting momenta:" << endl;
*gmsg << "* Initial gamma = " << itsReference.getGamma() << endl;
*gmsg << "* Initial beta = " << itsReference.getBeta() << endl;
*gmsg << "* Reference total momentum = " << referencePtot << " [beta gamma]" << endl;
- *gmsg << "* Reference azimuthal momentum (Pt) = " << referencePt << " [beta gamma]" << endl;
- *gmsg << "* Reference radial momentum (Pr) = " << referencePr << " [beta gamma]" << endl;
- *gmsg << "* Reference axial momentum (Pz) = " << referencePz << " [beta gamma]" << endl;
+ *gmsg << "* Reference azimuthal momentum (Pt) = " << referencePt << " [beta gamma]" << endl;
+ *gmsg << "* Reference radial momentum (Pr) = " << referencePr << " [beta gamma]" << endl;
+ *gmsg << "* Reference axial momentum (Pz) = " << referencePz << " [beta gamma]" << endl;
*gmsg << endl;
double sym = cycl_m->getSymmetry();
@@ -438,15 +437,15 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron& cycl) {
double rmin = cycl_m->getMinR();
double rmax = cycl_m->getMaxR();
- *gmsg << "* Radial aperture = " << rmin << " ... " << rmax<<" [m] "<< endl;
+ *gmsg << "* Radial aperture = " << rmin << " ... " << rmax << " [m]" << endl;
double zmin = cycl_m->getMinZ();
double zmax = cycl_m->getMaxZ();
- *gmsg << "* Vertical aperture = " << zmin << " ... " << zmax<<" [m]"<< endl;
+ *gmsg << "* Vertical aperture = " << zmin << " ... " << zmax << " [m]" << endl;
double h = cycl_m->getCyclHarm();
*gmsg << "* Number of trimcoils = " << cycl_m->getNumberOfTrimcoils() << endl;
- *gmsg << "* Harmonic number h = " << h << " " << endl;
+ *gmsg << "* Harmonic number h = " << h << endl;
std::vector<double> rffrequ = cycl_m->getRfFrequ();
*gmsg << "* RF frequency = " << Util::doubleVectorToString(rffrequ) << " [MHz]" << endl;
@@ -470,8 +469,8 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron& cycl) {
cycl_m->initialise(itsBunch_m, cycl_m->getBScale());
double BcParameter[8] = {};
- BcParameter[0] = Units::mm2m * cycl_m->getRmin();
- BcParameter[1] = Units::mm2m * cycl_m->getRmax();
+ BcParameter[0] = cycl_m->getRmin();
+ BcParameter[1] = cycl_m->getRmax();
// Store inner radius and outer radius of cyclotron field map in the list
buildupFieldList(BcParameter, ElementType::CYCLOTRON, cycl_m);
@@ -809,9 +808,9 @@ void ParallelCyclotronTracker::visitRFCavity(const RFCavity& as) {
double BcParameter[8] = {};
- BcParameter[0] = Units::mm2m * rmin;
- BcParameter[1] = Units::mm2m * rmax;
- BcParameter[2] = Units::mm2m * pdis;
+ BcParameter[0] = rmin;
+ BcParameter[1] = rmax;
+ BcParameter[2] = pdis;
BcParameter[3] = angle;
buildupFieldList(BcParameter, ElementType::RFCAVITY, elptr);
@@ -837,9 +836,9 @@ void ParallelCyclotronTracker::visitRing(const Ring& ring) {
referenceTheta = opalRing_m->getBeamPhiInit();
beamInitialRotation = opalRing_m->getBeamThetaInit();
- if (referenceTheta <= -180.0 || referenceTheta > 180.0) {
+ if (!Util::angleBetweenAngles(referenceTheta, -Physics::pi, Physics::pi)) {
throw OpalException("Error in ParallelCyclotronTracker::visitRing",
- "PHIINIT is out of [-180, 180)");
+ "BEAM_PHIINIT is out of [-180, 180)");
}
referenceZ = 0.0;
@@ -851,15 +850,15 @@ void ParallelCyclotronTracker::visitRing(const Ring& ring) {
if (referencePtot < 0.0)
referencePt *= -1.0;
- sinRefTheta_m = std::sin(referenceTheta * Units::deg2rad);
- cosRefTheta_m = std::cos(referenceTheta * Units::deg2rad);
+ sinRefTheta_m = std::sin(referenceTheta);
+ cosRefTheta_m = std::cos(referenceTheta);
double BcParameter[8] = {}; // zero initialise array
buildupFieldList(BcParameter, ElementType::RING, opalRing_m);
// Finally print some diagnostic
- *gmsg << "* Initial beam radius = " << referenceR << " [mm] " << endl;
+ *gmsg << "* Initial beam radius = " << referenceR << " [m] " << endl;
*gmsg << "* Initial gamma = " << itsReference.getGamma() << endl;
*gmsg << "* Initial beta = " << itsReference.getBeta() << endl;
*gmsg << "* Total reference momentum = " << referencePtot << " [beta gamma]" << endl;
@@ -1191,9 +1190,9 @@ void ParallelCyclotronTracker::execute() {
lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(bgf_m->getOpalName(),!Options::asciidump));
// External field arrays for dumping
- for (int k = 0; k < 2; k++)
+ for (int k = 0; k < 2; k++) {
FDext_m[k] = Vector_t(0.0, 0.0, 0.0);
-
+ }
extE_m = Vector_t(0.0, 0.0, 0.0);
extB_m = Vector_t(0.0, 0.0, 0.0);
DumpFields::writeFields((((*FieldDimensions.begin())->second).second));
@@ -1243,8 +1242,8 @@ void ParallelCyclotronTracker::MtsTracker() {
/*
* variable unit meaning
* ------------------------------------------------
- * t [ns] time
- * dt [ns] time step
+ * t [s] time
+ * dt [s] time step
* oldReferenceTheta [rad] azimuth angle
* itsBunch_m->R [m] particle position
*
@@ -1257,7 +1256,7 @@ void ParallelCyclotronTracker::MtsTracker() {
*gmsg << "* MTS: Number of substeps per step is " << numSubsteps << endl;
double const dt_inner = dt / double(numSubsteps);
- *gmsg << "* MTS: The inner time step is therefore " << dt_inner << " [ns]" << endl;
+ *gmsg << "* MTS: The inner time step is therefore " << dt_inner * Units::s2ns << " [ns]" << endl;
// int SteptoLastInj = itsBunch_m->getSteptoLastInj();
@@ -1265,8 +1264,9 @@ void ParallelCyclotronTracker::MtsTracker() {
*gmsg << "* ---------------------- Start tracking ---------------------------------- *" << endl;
- if ( itsBunch_m->hasFieldSolver() )
+ if (itsBunch_m->hasFieldSolver()) {
computeSpaceChargeFields_m();
+ }
for (; (step_m < maxSteps_m) && (itsBunch_m->getTotalNum()>0); step_m++) {
@@ -1284,9 +1284,9 @@ void ParallelCyclotronTracker::MtsTracker() {
}
// Substeps for external field integration
- for (int n = 0; n < numSubsteps; ++n)
+ for (int n = 0; n < numSubsteps; ++n) {
borisExternalFields(dt_inner);
-
+ }
// bunch injection
// TODO: Where is correct location for this piece of code? Beginning/end of step? Before field solve?
injectBunch(flagTransition);
@@ -1309,9 +1309,9 @@ void ParallelCyclotronTracker::MtsTracker() {
}
// Second half kick from space charge force
- if (itsBunch_m->hasFieldSolver())
+ if (itsBunch_m->hasFieldSolver()) {
kick(0.5 * dt);
-
+ }
// recalculate bingamma and reset the BinID for each particles according to its current gamma
if (isMultiBunch() && (step_m % Options::rebinFreq == 0)) {
mbHandler_m->updateParticleBins(itsBunch_m);
@@ -1319,7 +1319,6 @@ void ParallelCyclotronTracker::MtsTracker() {
// dump some data after one push in single particle tracking
if ( mode_m == TrackingMode::SINGLE ) {
-
unsigned int i = 0;
double temp_meanTheta = calculateAngle2(itsBunch_m->R[i](0),
@@ -1354,14 +1353,12 @@ void ParallelCyclotronTracker::MtsTracker() {
// printing + updating bunch parameters + updating t
update_m(t, dt, finishedTurn);
-
}
// Some post-integration stuff
*gmsg << endl;
*gmsg << "* ---------------------- DONE TRACKING PARTICLES ------------------------- *" << endl;
-
//FIXME
dvector_t Ttime = dvector_t();
dvector_t Tdeltr = dvector_t();
@@ -1376,8 +1373,8 @@ void ParallelCyclotronTracker::GenericTracker() {
/*
* variable unit meaning
* ------------------------------------------------
- * t [ns] time
- * dt [ns] time step
+ * t [s] time
+ * dt [s] time step
* oldReferenceTheta [rad] azimuth angle
* itsBunch_m->R [m] particle position
*
@@ -1480,7 +1477,6 @@ bool ParallelCyclotronTracker::getFieldsAtPoint(const double& t, const size_t& P
*
* @return
*/
-
bool ParallelCyclotronTracker::checkGapCross(Vector_t Rold, Vector_t Rnew,
RFCavity* rfcavity, double& Dold) {
bool flag = false;
@@ -1976,21 +1972,19 @@ inline void ParallelCyclotronTracker::getQuaternionTwoVectors(Vector_t u, Vector
bool ParallelCyclotronTracker::push(double h) {
- /* h [ns]
- * dt1 [ns]
- * dt2 [ns]
- */
- IpplTimings::startTimer(IntegrationTimer_m);
- h *= Units::ns2s;
+ IpplTimings::startTimer(IntegrationTimer_m);
bool flagNeedUpdate = false;
for (unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) {
+
Vector_t const oldR = itsBunch_m->R[i];
- double const gamma = std::sqrt(1.0 + dot(itsBunch_m->P[i], itsBunch_m->P[i]));
+ double const gamma = Util::getGamma(itsBunch_m->P[i]);
double const c_gamma = Physics::c / gamma;
Vector_t const v = itsBunch_m->P[i] * c_gamma;
+
itsBunch_m->R[i] += h * v;
+
for (const auto & ccd : cavCrossDatas_m) {
double const distNew = (itsBunch_m->R[i][0] * ccd.sinAzimuth - itsBunch_m->R[i][1] * ccd.cosAzimuth) - ccd.perpenDistance;
bool tagCrossing = false;
@@ -2000,16 +1994,14 @@ bool ParallelCyclotronTracker::push(double h) {
if (distOld > 0.0) tagCrossing = true;
}
if (tagCrossing) {
- double const dt1 = distOld / std::sqrt(dot(v, v));
+ double const dt1 = distOld / euclidean_norm(v);
double const dt2 = h - dt1;
- // Retrack particle from the old postion to cavity gap point
+ // Re-track particle from the old position to cavity gap point
itsBunch_m->R[i] = oldR + dt1 * v;
// Momentum kick
- //itsBunch_m->R[i] *= 1000.0; // RFkick uses [itsBunch_m->R] == mm
- RFkick(ccd.cavity, itsBunch_m->getT() * Units::s2ns, dt1, i);
- //itsBunch_m->R[i] *= 0.001;
+ RFkick(ccd.cavity, itsBunch_m->getT(), dt1, i);
itsBunch_m->R[i] += dt2 * itsBunch_m->P[i] * c_gamma;
}
@@ -2034,7 +2026,7 @@ bool ParallelCyclotronTracker::kick(double h) {
pusher.kick(itsBunch_m->R[i], itsBunch_m->P[i],
itsBunch_m->Ef[i], itsBunch_m->Bf[i],
- h * Units::ns2s, M, q);
+ h, M, q);
flagNeedUpdate |= (itsBunch_m->Bin[i] < 0);
}
IpplTimings::stopTimer(IntegrationTimer_m);
@@ -2043,7 +2035,6 @@ bool ParallelCyclotronTracker::kick(double h) {
void ParallelCyclotronTracker::borisExternalFields(double h) {
- // h in [ns]
// push particles for first half step
bool flagNeedUpdate = push(0.5 * h);
@@ -2055,7 +2046,7 @@ void ParallelCyclotronTracker::borisExternalFields(double h) {
itsBunch_m->Ef[i] = Vector_t(0.0, 0.0, 0.0);
itsBunch_m->Bf[i] = Vector_t(0.0, 0.0, 0.0);
- computeExternalFields_m(i, itsBunch_m->getT() * Units::s2ns,
+ computeExternalFields_m(i, itsBunch_m->getT(),
itsBunch_m->Ef[i], itsBunch_m->Bf[i]);
}
IpplTimings::stopTimer(IntegrationTimer_m);
@@ -2068,7 +2059,7 @@ void ParallelCyclotronTracker::borisExternalFields(double h) {
// apply the plugin elements: probe, collimator, stripper, septum
flagNeedUpdate |= applyPluginElements(h);
- // destroy particles if they are marked as Bin=-1 in the plugin elements or out of global apeture
+ // destroy particles if they are marked as Bin=-1 in the plugin elements or out of global aperture
deleteParticle(flagNeedUpdate);
}
@@ -2204,15 +2195,13 @@ bool ParallelCyclotronTracker::deleteParticle(bool flagNeedUpdate) {
double const phi = calculateAngle(meanP(0), meanP(1)) - 0.5 * Physics::pi;
// Bunch (local) elevation at meanR w.r.t. xy plane
- double const psi = 0.5 * Physics::pi - std::acos(meanP(2) / std::sqrt(dot(meanP, meanP)));
+ double const psi = 0.5 * Physics::pi - std::acos(meanP(2) / euclidean_norm(meanP));
// For statistics data, the bunch is transformed into a local coordinate system
// with meanP in direction of y-axis -DW
globalToLocal(itsBunch_m->R, phi, psi, meanR);
globalToLocal(itsBunch_m->P, phi, psi, Vector_t(0.0)); // P should be rotated, but not shifted
- //itsBunch_m->R *= Vector_t(0.001); // mm --> m
-
// Now destroy particles and update pertinent parameters in local frame
// Note that update() will be called within boundp() -DW
itsBunch_m->boundp();
@@ -2220,8 +2209,6 @@ bool ParallelCyclotronTracker::deleteParticle(bool flagNeedUpdate) {
itsBunch_m->calcBeamParameters();
- //itsBunch_m->R *= Vector_t(1000.0); // m --> mm
-
localToGlobal(itsBunch_m->R, phi, psi, meanR);
localToGlobal(itsBunch_m->P, phi, psi, Vector_t(0.0));
@@ -2242,7 +2229,6 @@ void ParallelCyclotronTracker::initTrackOrbitFile() {
outfTrackOrbit_m.precision(8);
if (myNode_m == 0) {
-
if (OpalData::getInstance()->inRestartRun()) {
outfTrackOrbit_m.open(f.c_str(), std::ios::app);
outfTrackOrbit_m << "# Restart at integration step " << itsBunch_m->getLocalTrackStep() << std::endl;
@@ -2259,7 +2245,7 @@ void ParallelCyclotronTracker::initDistInGlobalFrame() {
if (!OpalData::getInstance()->inRestartRun()) {
// Start a new run (no restart)
- double const initialReferenceTheta = referenceTheta * Units::deg2rad;
+ double const initialReferenceTheta = referenceTheta;
// TODO: Replace with TracerParticle
// Force the initial phase space values of the particle with ID = 0 to zero,
@@ -2274,11 +2260,9 @@ void ParallelCyclotronTracker::initDistInGlobalFrame() {
}
// Initialize global R
- //itsBunch_m->R *= Vector_t(1000.0); // m --> mm
-
- Vector_t const initMeanR = Vector_t(Units::mm2m * referenceR * cosRefTheta_m,
- Units::mm2m * referenceR * sinRefTheta_m,
- Units::mm2m * referenceZ);
+ Vector_t const initMeanR = Vector_t(referenceR * cosRefTheta_m,
+ referenceR * sinRefTheta_m,
+ referenceZ);
localToGlobal(itsBunch_m->R, initialReferenceTheta, initMeanR);
@@ -2293,7 +2277,7 @@ void ParallelCyclotronTracker::initDistInGlobalFrame() {
// Out of the three coordinates of meanR (R, Theta, Z) only the angle
// changes the momentum vector...
- double angle = initialReferenceTheta+beamInitialRotation*Units::deg2rad;
+ double angle = initialReferenceTheta + beamInitialRotation;
localToGlobal(itsBunch_m->P, angle);
DistributionType distType = itsBunch_m->getDistType();
@@ -2322,11 +2306,10 @@ void ParallelCyclotronTracker::initDistInGlobalFrame() {
if ((Options::psDumpFrame != DumpFrame::GLOBAL)) {
*gmsg << "* Restart in the local frame" << endl;
- //itsBunch_m->R *= Vector_t(1000.0); // m --> mm
- Vector_t const initMeanR = Vector_t(Units::mm2m * referenceR * cosRefTheta_m,
- Units::mm2m * referenceR * sinRefTheta_m,
- Units::mm2m * referenceZ);
+ Vector_t const initMeanR = Vector_t(referenceR * cosRefTheta_m,
+ referenceR * sinRefTheta_m,
+ referenceZ);
// Do the transformations
localToGlobal(itsBunch_m->R, referencePhi, referencePsi, initMeanR);
@@ -2342,7 +2325,6 @@ void ParallelCyclotronTracker::initDistInGlobalFrame() {
*gmsg << "* Restart in the global frame" << endl;
pathLength_m = itsBunch_m->get_sPos();
- //itsBunch_m->R *= Vector_t(1000.0); // m --> mm
}
}
@@ -2357,9 +2339,9 @@ void ParallelCyclotronTracker::initDistInGlobalFrame() {
double const phi = calculateAngle(meanP(0), meanP(1)) - 0.5 * Physics::pi;
// Bunch (local) elevation at meanR w.r.t. xy plane
- double const psi = 0.5 * Physics::pi - std::acos(meanP(2) / std::sqrt(dot(meanP, meanP)));
+ double const psi = 0.5 * Physics::pi - std::acos(meanP(2) / euclidean_norm(meanP));
- double radius = std::sqrt(meanR[0] * meanR[0] + meanR[1] * meanR[1]); // [m]
+ double radius = std::sqrt(meanR[0] * meanR[0] + meanR[1] * meanR[1]);
if ( isMultiBunch() ) {
mbHandler_m->setRadiusTurns(radius);
@@ -2369,8 +2351,6 @@ void ParallelCyclotronTracker::initDistInGlobalFrame() {
globalToLocal(itsBunch_m->R, phi, psi, meanR);
globalToLocal(itsBunch_m->P, phi, psi); // P should be rotated, but not shifted
- //itsBunch_m->R *= Vector_t(0.001); // mm --> m
-
itsBunch_m->boundp();
checkNumPart(std::string("\n* Before repartition: "));
@@ -2382,8 +2362,6 @@ void ParallelCyclotronTracker::initDistInGlobalFrame() {
*gmsg << endl << "* *********************** Bunch information in local frame: ************************";
*gmsg << *itsBunch_m << endl;
- //itsBunch_m->R *= Vector_t(1000.0); // m --> mm
-
localToGlobal(itsBunch_m->R, phi, psi, meanR);
localToGlobal(itsBunch_m->P, phi, psi);
@@ -2397,10 +2375,7 @@ void ParallelCyclotronTracker::initDistInGlobalFrame() {
step_m += 1;
}
- // Print out the Bunch information at beginning of the run. Because the bunch information
- // displays in units of m we have to change back and forth one more time.
- //itsBunch_m->R *= Vector_t(0.001); // mm --> m
-
+ // Print out the Bunch information at beginning of the run.
itsBunch_m->calcBeamParameters();
// multi-bunch simulation only
@@ -2408,8 +2383,6 @@ void ParallelCyclotronTracker::initDistInGlobalFrame() {
*gmsg << endl << "* *********************** Bunch information in global frame: ***********************";
*gmsg << *itsBunch_m << endl;
-
- //itsBunch_m->R *= Vector_t(1000.0); // m --> mm
}
void ParallelCyclotronTracker::setTimeStep() {
@@ -2601,7 +2574,7 @@ void ParallelCyclotronTracker::bunchDumpStatData(){
phi = calculateAngle(meanP(0), meanP(1)) - 0.5 * Physics::pi;
// Bunch (local) elevation at meanR w.r.t. xy plane
- psi = 0.5 * Physics::pi - std::acos(meanP(2) / std::sqrt(dot(meanP, meanP)));
+ psi = 0.5 * Physics::pi - std::acos(meanP(2) / euclidean_norm(meanP));
// Rotate so Pmean is in positive y direction. No shift, so that normalized emittance and
// unnormalized emittance as well as centroids are calculated correctly in
@@ -2622,7 +2595,7 @@ void ParallelCyclotronTracker::bunchDumpStatData(){
}
// --------------------------------- Get some Values ---------------------------------------- //
- double const temp_t = itsBunch_m->getT() * Units::s2ns;
+ double const temp_t = itsBunch_m->getT();
Vector_t meanR;
Vector_t meanP;
if (Options::psDumpFrame == DumpFrame::BUNCH_MEAN) {
@@ -2657,7 +2630,7 @@ void ParallelCyclotronTracker::bunchDumpStatData(){
phi = calculateAngle(meanP(0), meanP(1)) - 0.5 * Physics::pi;
// Bunch (local) elevation at meanR w.r.t. xy plane
- psi = 0.5 * Physics::pi - std::acos(meanP(2) / std::sqrt(dot(meanP, meanP)));
+ psi = 0.5 * Physics::pi - std::acos(meanP(2) / euclidean_norm(meanP));
// Rotate so Pmean is in positive y direction. No shift, so that normalized emittance and
// unnormalized emittance as well as centroids are calculated correctly in
@@ -2668,16 +2641,12 @@ void ParallelCyclotronTracker::bunchDumpStatData(){
globalToLocal(itsBunch_m->P, phi, psi);
}
- //itsBunch_m->R *= Vector_t(0.001); // mm -> m
-
FDext_m[0] = extB_m * Units::kG2T;
FDext_m[1] = extE_m; // kV/mm? -DW
// Save the stat file
itsDataSink->dumpSDDS(itsBunch_m, FDext_m, azimuth_m);
- //itsBunch_m->R *= Vector_t(1000.0); // m -> mm
-
// If we are in local mode, transform back after saving
if (Options::psDumpFrame != DumpFrame::GLOBAL) {
localToGlobal(itsBunch_m->R, phi, psi);
@@ -2698,7 +2667,7 @@ void ParallelCyclotronTracker::bunchDumpPhaseSpaceData() {
IpplTimings::startTimer(DumpTimer_m);
// --------------------------------- Get some Values ---------------------------------------- //
- double const temp_t = itsBunch_m->getT() * Units::s2ns;
+ double const temp_t = itsBunch_m->getT();
Vector_t meanR;
Vector_t meanP;
@@ -2712,7 +2681,7 @@ void ParallelCyclotronTracker::bunchDumpPhaseSpaceData() {
meanP = itsBunch_m->P[0];
}
- double const betagamma_temp = std::sqrt(dot(meanP, meanP));
+ double const betagamma_temp = euclidean_norm(meanP);
double const E = itsBunch_m->get_meanKineticEnergy();
@@ -2723,13 +2692,13 @@ void ParallelCyclotronTracker::bunchDumpPhaseSpaceData() {
double const phi = calculateAngle(meanP(0), meanP(1)) - 0.5 * Physics::pi;
// Bunch (local) elevation at meanR w.r.t. xy plane
- double const psi = 0.5 * Physics::pi - std::acos(meanP(2) / std::sqrt(dot(meanP, meanP)));
+ double const psi = 0.5 * Physics::pi - std::acos(meanP(2) / betagamma_temp);
// ---------------- Re-calculate reference values in format of input values ----------------- //
// Position:
// New OPAL 2.0: Init in m (back to mm just for output) -DW
referenceR = computeRadius(meanR); // includes m --> mm conversion
- referenceTheta = theta / Units::deg2rad;
+ referenceTheta = Units::rad2deg * theta;
referenceZ = Units::m2mm * meanR(2);
// Momentum in Theta-hat, R-hat, Z-hat coordinates at position meanR:
@@ -2763,7 +2732,6 @@ void ParallelCyclotronTracker::bunchDumpPhaseSpaceData() {
// The bunch is transformed into a local coordinate system with meanP in direction of y-axis
globalToLocal(itsBunch_m->R, phi, psi, meanR);
- //globalToLocal(itsBunch_m->R, phi, psi, meanR * Vector_t(0.001));
globalToLocal(itsBunch_m->P, phi, psi); // P should only be rotated
globalToLocal(extB_m, phi, psi);
@@ -2781,16 +2749,15 @@ void ParallelCyclotronTracker::bunchDumpPhaseSpaceData() {
referenceR,
referenceTheta,
referenceZ,
- phi / Units::deg2rad, // P_mean azimuth
+ phi * Units::rad2deg, // P_mean azimuth
// at ref. R/Th/Z
- psi / Units::deg2rad, // P_mean elevation
+ psi * Units::rad2deg, // P_mean elevation
// at ref. R/Th/Z
- dumpLocalFrame); // Flag localFrame
+ dumpLocalFrame); // Flag localFrame
if (dumpLocalFrame == true) {
// Return to global frame
localToGlobal(itsBunch_m->R, phi, psi, meanR);
- //localToGlobal(itsBunch_m->R, phi, psi, meanR * Vector_t(0.001));
localToGlobal(itsBunch_m->P, phi, psi);
}
@@ -2806,15 +2773,15 @@ void ParallelCyclotronTracker::bunchDumpPhaseSpaceData() {
}
// Print dump information on screen
- *gmsg << "* T = " << temp_t << " ns"
+ *gmsg << "* T = " << temp_t << " s"
<< ", Live Particles: " << itsBunch_m->getTotalNum() << endl;
*gmsg << "* E = " << E << " MeV"
<< ", beta * gamma = " << betagamma_temp << endl;
*gmsg << "* Bunch position: R = " << referenceR << " mm"
<< ", Theta = " << referenceTheta << " Deg"
<< ", Z = " << referenceZ << " mm" << endl;
- *gmsg << "* Local Azimuth = " << phi / Units::deg2rad << " Deg"
- << ", Local Elevation = " << psi / Units::deg2rad << " Deg" << endl;
+ *gmsg << "* Local Azimuth = " << phi * Units::rad2deg << " Deg"
+ << ", Local Elevation = " << psi * Units::rad2deg << " Deg" << endl;
IpplTimings::stopTimer(DumpTimer_m);
}
@@ -2868,8 +2835,8 @@ void ParallelCyclotronTracker::update_m(double& t, const double& dt,
std::tuple<double, double, double> ParallelCyclotronTracker::initializeTracking_m() {
// Read in some control parameters
- setup_m.scSolveFreq = (spiral_flag) ? 1 : Options::scSolveFreq;
- setup_m.stepsPerTurn = itsBunch_m->getStepsPerTurn();
+ setup_m.scSolveFreq = (spiral_flag) ? 1 : Options::scSolveFreq;
+ setup_m.stepsPerTurn = itsBunch_m->getStepsPerTurn();
// Define 3 special azimuthal angles where dump particle's six parameters
// at each turn into 3 ASCII files. only important in single particle tracking
@@ -2879,10 +2846,10 @@ std::tuple<double, double, double> ParallelCyclotronTracker::initializeTracking_
azimuth_angle_m[2] = 45.0 * Units::deg2rad;
double harm = getHarmonicNumber();
- double dt = itsBunch_m->getdT() * Units::s2ns * harm;
- double t = itsBunch_m->getT() * Units::s2ns;
+ double dt = itsBunch_m->getdT() * harm;
+ double t = itsBunch_m->getT();
- double oldReferenceTheta = referenceTheta * Units::deg2rad; // init here, reset each step
+ double oldReferenceTheta = referenceTheta; // init here, reset each step
setup_m.deltaTheta = Physics::pi / (setup_m.stepsPerTurn); // half of the average angle per step
//int stepToLastInj = itsBunch_m->getSteptoLastInj(); // TODO: Do we need this? -DW
@@ -2917,12 +2884,12 @@ std::tuple<double, double, double> ParallelCyclotronTracker::initializeTracking_
}
// --- Output to user --- //
- *gmsg << "* Beginning of this run is at t = " << t << " [ns]" << endl;
- *gmsg << "* The time step is set to dt = " << dt << " [ns]" << endl;
+ *gmsg << "* Beginning of this run is at t = " << t * Units::s2ns << " [ns]" << endl;
+ *gmsg << "* The time step is set to dt = " << dt * Units::s2ns << " [ns]" << endl;
if ( isMultiBunch() ) {
*gmsg << "* MBM: Time interval between neighbour bunches is set to "
- << setup_m.stepsPerTurn * dt << "[ns]" << endl;
+ << setup_m.stepsPerTurn * dt * Units::s2ns << "[ns]" << endl;
*gmsg << "* MBM: The particles energy bin reset frequency is set to "
<< Options::rebinFreq << endl;
}
@@ -3088,7 +3055,7 @@ void ParallelCyclotronTracker::seoMode_m(double& t, const double dt, bool& /*fin
// store dx and dz for future tune calculation if higher precision needed, reduce freqSample.
if (step_m % Options::sptDumpFreq == 0) {
- Ttime.push_back(t * Units::ns2s);
+ Ttime.push_back(t);
Tdeltz.push_back(z_tuning[1]);
Tdeltr.push_back(r_tuning[1] - r_tuning[0]);
TturnNumber.push_back(turnnumber_m);
@@ -3139,17 +3106,18 @@ void ParallelCyclotronTracker::singleMode_m(double& t, const double dt,
oldReferenceTheta = temp_meanTheta;
// used for gap crossing checking
- Vector_t Rold = itsBunch_m->R[i]; // [x,y,z] (mm)
+ Vector_t Rold = itsBunch_m->R[i]; // [x,y,z] (m)
Vector_t Pold = itsBunch_m->P[i]; // [px,py,pz] (beta*gamma)
// integrate for one step in the lab Cartesian frame (absolute value).
itsStepper_mp->advance(itsBunch_m, i, t, dt);
// If gap crossing happens, do momenta kicking (not if gap crossing just happened)
- if (itsBunch_m->cavityGapCrossed[i] == true)
+ if (itsBunch_m->cavityGapCrossed[i] == true) {
itsBunch_m->cavityGapCrossed[i] = false;
- else
+ } else {
gapCrossKick_m(i, t, dt, Rold, Pold);
+ }
IpplTimings::stopTimer(IntegrationTimer_m);
// Destroy particles if they are marked as Bin = -1 in the plugin elements
@@ -3174,7 +3142,7 @@ void ParallelCyclotronTracker::bunchMode_m(double& t, const double dt, bool& fin
// oldReferenceTheta = calculateAngle2(meanR(0), meanR(1));
// Calculate SC field before each time step and keep constant during integration.
- // Space Charge effects are included only when total macropaticles number is NOT LESS THAN 1000.
+ // Space Charge effects are included only when total macroparticles number is NOT LESS THAN 1000.
if (itsBunch_m->hasFieldSolver()) {
if (step_m % setup_m.scSolveFreq == 0) {
@@ -3213,7 +3181,7 @@ void ParallelCyclotronTracker::bunchMode_m(double& t, const double dt, bool& fin
for (size_t i = 0; i < itsBunch_m->getLocalNum(); i++) {
// used for gap crossing checking
- Vector_t Rold = itsBunch_m->R[i]; // [x,y,z] (mm)
+ Vector_t Rold = itsBunch_m->R[i]; // [x,y,z] (m)
Vector_t Pold = itsBunch_m->P[i]; // [px,py,pz] (beta*gamma)
// Integrate for one step in the lab Cartesian frame (absolute value).
@@ -3275,11 +3243,8 @@ void ParallelCyclotronTracker::gapCrossKick_m(size_t i, double t,
if ( tag_crossing ) {
itsBunch_m->cavityGapCrossed[i] = true;
- double oldMomentum2 = dot(Pold, Pold);
- double oldBetgam = std::sqrt(oldMomentum2);
- double oldGamma = std::sqrt(1.0 + oldMomentum2);
- double oldBeta = oldBetgam / oldGamma;
- double dt1 = DistOld / (oldBeta * Physics::c / Units::s2ns); // ns, c in [m/ns]
+ double oldBeta = euclidean_norm(Util::getBeta(Pold));
+ double dt1 = DistOld / (oldBeta * Physics::c);
double dt2 = dt - dt1;
// retrack particle from the old postion to cavity gap point
@@ -3315,7 +3280,7 @@ void ParallelCyclotronTracker::dumpAzimuthAngles_m(const double& t,
( temp_meanTheta >= azimuth_angle_m[i] - setup_m.deltaTheta))
{
*(outfTheta_m[i]) << "#Turn number = " << turnnumber_m
- << ", Time = " << t << " [ns]" << std::endl
+ << ", Time = " << t * Units::s2ns << " [ns]" << std::endl
<< " " << std::hypot(R(0), R(1))
<< " " << P(0) * std::cos(temp_meanTheta) + P(1) * std::sin(temp_meanTheta)
<< " " << temp_meanTheta * Units::rad2deg
@@ -3338,7 +3303,7 @@ void ParallelCyclotronTracker::dumpThetaEachTurn_m(const double& t,
*gmsg << "* SPT: Finished turn " << turnnumber_m - 1 << endl;
*(outfTheta_m[3]) << "#Turn number = " << turnnumber_m
- << ", Time = " << t << " [ns]" << std::endl
+ << ", Time = " << t * Units::s2ns << " [ns]" << std::endl
<< " " << std::sqrt(R(0) * R(0) + R(1) * R(1))
<< " " << P(0) * std::cos(temp_meanTheta) +
P(1) * std::sin(temp_meanTheta)
@@ -3362,32 +3327,25 @@ void ParallelCyclotronTracker::computeSpaceChargeFields_m() {
// If we are doing a spiral inflector simulation and are using the SAAMG solver
// we don't rotate or translate the bunch and gamma is 1.0 (non-relativistic).
- //itsBunch_m->R *= Vector_t(0.001); // mm --> m
-
itsBunch_m->setGlobalMeanR(Vector_t(0.0, 0.0, 0.0));
itsBunch_m->setGlobalToLocalQuaternion(Quaternion_t(1.0, 0.0, 0.0, 0.0));
itsBunch_m->computeSelfFields_cycl(1.0);
- //itsBunch_m->R *= Vector_t(1000.0); // m --> mm
-
} else {
- Vector_t const meanR = calcMeanR(); // (m)
+ Vector_t const meanR = calcMeanR();
PreviousMeanP = calcMeanP();
// Since calcMeanP takes into account all particles of all bins (TODO: Check this! -DW)
// Using the quaternion method with PreviousMeanP and yaxis should give the correct result
-
Quaternion_t quaternionToYAxis;
getQuaternionTwoVectors(PreviousMeanP, yaxis, quaternionToYAxis);
globalToLocal(itsBunch_m->R, quaternionToYAxis, meanR);
- //itsBunch_m->R *= Vector_t(0.001); // mm --> m
-
if ((step_m + 1) % Options::boundpDestroyFreq == 0) {
itsBunch_m->boundp_destroyCycl();
} else {
@@ -3404,9 +3362,7 @@ void ParallelCyclotronTracker::computeSpaceChargeFields_m() {
// Calculate space charge field for each energy bin
for (int b = 0; b < itsBunch_m->getLastemittedBin(); b++) {
-
itsBunch_m->setBinCharge(b, itsBunch_m->getChargePerParticle());
- //itsBunch_m->setGlobalMeanR(0.001 * meanR);
itsBunch_m->setGlobalMeanR(meanR);
itsBunch_m->setGlobalToLocalQuaternion(quaternionToYAxis);
itsBunch_m->computeSelfFields_cycl(b);
@@ -3416,21 +3372,15 @@ void ParallelCyclotronTracker::computeSpaceChargeFields_m() {
} else {
// --- Single bunch mode --- //
-
- double temp_meangamma = std::sqrt(1.0 + dot(PreviousMeanP, PreviousMeanP));
+ double temp_meangamma = Util::getGamma(PreviousMeanP);
repartition();
- //itsBunch_m->setGlobalMeanR(0.001 * meanR);
itsBunch_m->setGlobalMeanR(meanR);
itsBunch_m->setGlobalToLocalQuaternion(quaternionToYAxis);
-
itsBunch_m->computeSelfFields_cycl(temp_meangamma);
}
- //scale coordinates back
- //itsBunch_m->R *= Vector_t(1000.0); // m --> mm
-
// Transform coordinates back to global
localToGlobal(itsBunch_m->R, quaternionToYAxis, meanR);
@@ -3523,7 +3473,7 @@ void ParallelCyclotronTracker::saveInjectValues() {
MultiBunchHandler::injection_t& inj = mbHandler_m->getInjectionValues();
- inj.time = itsBunch_m->getT() * Units::s2ns;
+ inj.time = itsBunch_m->getT();
inj.pathlength = itsBunch_m->get_sPos();
inj.azimuth = azimuth_m;
inj.radius = radius;
@@ -3552,10 +3502,9 @@ void ParallelCyclotronTracker::updatePathLength(const double& dt) {
void ParallelCyclotronTracker::updateTime(const double& dt) {
- // t is in seconds
double t = itsBunch_m->getT();
- itsBunch_m->setT(t + dt * Units::ns2s);
+ itsBunch_m->setT(t + dt);
if ( isMultiBunch() ) {
mbHandler_m->updateTime(dt);
diff --git a/src/Algorithms/ParallelCyclotronTracker.h b/src/Algorithms/ParallelCyclotronTracker.h
index fc27ffdd1762d3cf5db6ab2fae33e21f139dd39d..a612ec2d433456ea18e5822db3c8d1c3bc478017 100644
--- a/src/Algorithms/ParallelCyclotronTracker.h
+++ b/src/Algorithms/ParallelCyclotronTracker.h
@@ -479,7 +479,7 @@ private:
void update_m(double& t, const double& dt, const bool& finishedTurn);
/*!
- * @returns the time t [ns], time step dt [ns] and the azimuth angle [rad]
+ * @returns the time t [s], time step dt [s] and the azimuth angle [rad]
*/
std::tuple<double, double, double> initializeTracking_m();
@@ -530,12 +530,12 @@ private:
bool hasMultiBunch() const;
/*
- * @param dt time step in ns
+ * @param dt time step in s
*/
void updatePathLength(const double& dt);
/*
- * @param dt time step in ns
+ * @param dt time step in s
*/
void updateTime(const double& dt);
diff --git a/src/Classic/AbsBeamline/Cyclotron.cpp b/src/Classic/AbsBeamline/Cyclotron.cpp
index 4838bcbbdae824b73d80ac31bc5e19f1e1aa476d..01207240ff4928df54fddb6015f108fc8d683305 100644
--- a/src/Classic/AbsBeamline/Cyclotron.cpp
+++ b/src/Classic/AbsBeamline/Cyclotron.cpp
@@ -162,7 +162,7 @@ double Cyclotron::getPZinit() const {
}
void Cyclotron::setTrimCoilThreshold(double trimCoilThreshold) {
- trimCoilThreshold_m = Units::T2kG * trimCoilThreshold;
+ trimCoilThreshold_m = trimCoilThreshold;
}
double Cyclotron::getTrimCoilThreshold() const {
@@ -308,15 +308,11 @@ double Cyclotron::getRmax() const {
}
void Cyclotron::setMinR(double r) {
- // DW: This is to let the user keep using mm in the input file for now
- // while switching internally to m
- minr_m = Units::mm2m * r;
+ minr_m = r;
}
void Cyclotron::setMaxR(double r) {
- // DW: This is to let the user keep using mm in the input file for now
- // while switching internally to m
- maxr_m = Units::mm2m * r;
+ maxr_m = r;
}
double Cyclotron::getMinR() const {
@@ -338,9 +334,7 @@ double Cyclotron::getMaxR() const {
}
void Cyclotron::setMinZ(double z) {
- // DW: This is to let the user keep using mm in the input file for now
- // while switching internally to m
- minz_m = Units::mm2m * z;
+ minz_m = z;
}
double Cyclotron::getMinZ() const {
@@ -348,9 +342,7 @@ double Cyclotron::getMinZ() const {
}
void Cyclotron::setMaxZ(double z) {
- // DW: This is to let the user keep using mm in the input file for now
- // while switching internally to m
- maxz_m = Units::mm2m * z;
+ maxz_m = z;
}
double Cyclotron::getMaxZ() const {
@@ -409,10 +401,6 @@ Cyclotron::BFieldType Cyclotron::getBFieldType() const {
void Cyclotron::checkInitialReferenceParticle(double refR,
double refTheta,
double refZ) {
- refZ *= Units::mm2m;
- refR *= Units::mm2m;
- refTheta *= Units::deg2rad;
-
if ( (refZ < minz_m || refZ > maxz_m) ||
(refR < minr_m || refR > maxr_m) ) {
throw GeneralClassicException(
@@ -437,23 +425,23 @@ bool Cyclotron::apply(const size_t& id, const double& t, Vector_t& E, Vector_t&
if (zpos > maxz_m || zpos < minz_m || rpos > maxr_m || rpos < minr_m) {
flagNeedUpdate = true;
*gmsgALL << level4 << getName() << ": Particle " << id
- << " out of the global aperture of cyclotron!" << endl;
+ << " out of the global aperture of cyclotron!" << endl;
*gmsgALL << level4 << getName()
- << ": Coords: "<< RefPartBunch_m->R[id] << " m" << endl;
+ << ": Coords: "<< RefPartBunch_m->R[id] << " m" << endl;
} else {
flagNeedUpdate = apply(RefPartBunch_m->R[id], RefPartBunch_m->P[id], t, E, B);
if (flagNeedUpdate) {
*gmsgALL << level4 << getName() << ": Particle "<< id
- << " out of the field map boundary!" << endl;
+ << " out of the field map boundary!" << endl;
*gmsgALL << level4 << getName()
- << ": Coords: "<< RefPartBunch_m->R[id] << " m" << endl;
+ << ": Coords: "<< RefPartBunch_m->R[id] << " m" << endl;
}
}
if (flagNeedUpdate) {
lossDs_m->addParticle(OpalParticle(id, RefPartBunch_m->R[id], RefPartBunch_m->P[id],
- t*Units::ns2s, RefPartBunch_m->Q[id], RefPartBunch_m->M[id]),
+ t, RefPartBunch_m->Q[id], RefPartBunch_m->M[id]),
std::make_pair(0, RefPartBunch_m->bunchNum[id]));
RefPartBunch_m->Bin[id] = -1;
}
@@ -465,7 +453,7 @@ bool Cyclotron::apply(const Vector_t& R, const Vector_t& /*P*/,
const double& t, Vector_t& E, Vector_t& B) {
double tet = 0.0;
- if (std::fabs(R[0]) < 1.0E-10) {
+ if (std::abs(R[0]) < 1.0E-10) {
if (R[1] >= 0.0) {
tet = Physics::pi / 2.0;
} else {
@@ -480,19 +468,15 @@ bool Cyclotron::apply(const Vector_t& R, const Vector_t& /*P*/,
tet = Physics::two_pi + std::atan(R[1] / R[0]);
}
}
- double tet_rad = tet;
-
- // the actual angle of particle in degree
- tet *= Units::rad2deg;
// Necessary for gap phase output -DW
- if (0 <= tet && tet <= 45) waitingGap_m = 1;
+ if ( 0 <= tet && tet <= (Physics::pi / 4) ) waitingGap_m = 1;
// dB_{z}/dr, dB_{z}/dtheta, B_{z}
double brint = 0.0, btint = 0.0, bzint = 0.0;
- const double rad = std::hypot(R[0],R[1]);
- if ( this->interpolate(rad, tet_rad, brint, btint, bzint) ) {
+ const double rad = std::hypot(R[0],R[1]);
+ if ( this->interpolate(rad, tet, brint, btint, bzint) ) {
/* Br */
double br = - brint * R[2];
@@ -503,11 +487,11 @@ bool Cyclotron::apply(const Vector_t& R, const Vector_t& /*P*/,
/* Bz */
double bz = - bzint;
- this->applyTrimCoil(rad, R[2], tet_rad, br, bz);
+ this->applyTrimCoil(rad, R[2], tet, br, bz);
/* Br Btheta -> Bx By */
- B[0] = br * std::cos(tet_rad) - bt * std::sin(tet_rad);
- B[1] = br * std::sin(tet_rad) + bt * std::cos(tet_rad);
+ B[0] = br * std::cos(tet) - bt * std::sin(tet);
+ B[1] = br * std::sin(tet) + bt * std::cos(tet);
B[2] = bz;
} else {
return true;
@@ -518,10 +502,10 @@ bool Cyclotron::apply(const Vector_t& R, const Vector_t& /*P*/,
}
//The RF field is supposed to be sampled on a cartesian grid
- std::vector<Fieldmap>::const_iterator fi = RFfields_m.begin();
- std::vector<double>::const_iterator rffi = rffrequ_m.begin();
- std::vector<double>::const_iterator rfphii = rfphi_m.begin();
- std::vector<double>::const_iterator escali = escale_m.begin();
+ std::vector<Fieldmap>::const_iterator fi = RFfields_m.begin();
+ std::vector<double>::const_iterator rffi = rffrequ_m.begin();
+ std::vector<double>::const_iterator rfphii = rfphi_m.begin();
+ std::vector<double>::const_iterator escali = escale_m.begin();
std::vector<bool>::const_iterator superposei = superpose_m.begin();
std::vector<std::vector<double>>::const_iterator rffci = rffc_m.begin();
std::vector<std::vector<double>>::const_iterator rfvci = rfvc_m.begin();
@@ -553,17 +537,17 @@ bool Cyclotron::apply(const Vector_t& R, const Vector_t& /*P*/,
if (fieldType_m == BFieldType::SYNCHRO) {
double powert = 1;
for (const double fcoef : (*rffci)) {
- powert *= (t * Units::ns2s);
+ powert *= t;
frequency += fcoef * powert; // Add frequency ramp [MHz/s^n]
}
powert = 1;
for (const double vcoef : (*rfvci)) {
- powert *= (t * Units::ns2s);
+ powert *= t;
ebscale += vcoef * powert; // Add frequency ramp [MHz/s^n]
}
}
- double phase = Physics::two_pi * Units::MHz2Hz * frequency * t * Units::ns2s + (*rfphii);
+ double phase = Physics::two_pi * Units::MHz2Hz * frequency * t + (*rfphii);
E += ebscale * std::cos(phase) * tmpE;
B -= ebscale * std::sin(phase) * tmpB;
@@ -1146,7 +1130,6 @@ void Cyclotron::getFieldFromFile_FFA(const double& /*scaleFactor*/) {
for (int i = 0; i < num_of_header_lines; ++i)
file_to_read.ignore(max_num_of_char_in_a_line, '\n');
- // TEMP for OPAL 2.0 changing this to m -DW
while(!file_to_read.eof()) {
double r, th, x, y, bz;
file_to_read >> r >> th >> x >> y >> bz;
diff --git a/src/Classic/AbsBeamline/PluginElement.cpp b/src/Classic/AbsBeamline/PluginElement.cpp
index 6b620876bf6d19d1d76a19f4539c281b2ca41d28..c9e7e89837b3ef4b6737a0f49a1cf755157fd79e 100644
--- a/src/Classic/AbsBeamline/PluginElement.cpp
+++ b/src/Classic/AbsBeamline/PluginElement.cpp
@@ -116,11 +116,11 @@ void PluginElement::setDimensions(double xstart, double xend, double ystart, dou
void PluginElement::setGeom(const double dist) {
double slope;
- if (xend_m == xstart_m)
+ if (xend_m == xstart_m) {
slope = 1.0e12;
- else
+ } else {
slope = (yend_m - ystart_m) / (xend_m - xstart_m);
-
+ }
double coeff2 = std::sqrt(1 + slope * slope);
double coeff1 = slope / coeff2;
double halfdist = dist / 2.0;
@@ -143,9 +143,7 @@ void PluginElement::setGeom(const double dist) {
}
void PluginElement::changeWidth(PartBunchBase<double, 3> *bunch, int i, const double tstep, const double tangle) {
-
- constexpr double c_mtns = Physics::c / Units::s2ns; // m/s --> m/ns
- double lstep = euclidean_norm(bunch->P[i]) / Util::getGamma(bunch->P[i]) * c_mtns * tstep; // [m]
+ double lstep = euclidean_norm(bunch->P[i]) / Util::getGamma(bunch->P[i]) * Physics::c * tstep;
double sWidth = lstep / std::sqrt( 1 + 1/tangle/tangle );
setGeom(sWidth);
}
diff --git a/src/Classic/AbsBeamline/RFCavity.cpp b/src/Classic/AbsBeamline/RFCavity.cpp
index ba3e0ff7817c9478e99e8dd9c11098312604662c..e1a15560d0caedcc2a64325fe0ef9020b36cb74b 100644
--- a/src/Classic/AbsBeamline/RFCavity.cpp
+++ b/src/Classic/AbsBeamline/RFCavity.cpp
@@ -132,6 +132,7 @@ bool RFCavity::apply(const Vector_t& R,
if (R(2) >= startField_m &&
R(2) < startField_m + getElementLength()) {
+
Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
bool outOfBounds = fieldmap_m->getFieldstrength(R, tmpE, tmpB);
@@ -139,7 +140,6 @@ bool RFCavity::apply(const Vector_t& R,
E += (scale_m + scaleError_m) * std::cos(frequency_m * t + phase_m + phaseError_m) * tmpE;
B -= (scale_m + scaleError_m) * std::sin(frequency_m * t + phase_m + phaseError_m) * tmpB;
-
}
return false;
}
@@ -399,7 +399,8 @@ void RFCavity::getMomentaKick(const double normalRadius,
Voltage *= Ufactor;
// rad/s, ns --> rad
- double nphase = (frequency * (t + dtCorrt) * Units::ns2s) - phi0_m;
+
+ double nphase = (frequency * (t + dtCorrt)) - phi0_m;
double dgam = Voltage * std::cos(nphase) / (restMass);
double tempdegree = std::fmod(nphase * Units::rad2deg, 360.0);
@@ -412,7 +413,7 @@ void RFCavity::getMomentaKick(const double normalRadius,
double pr = momentum[0] * cosAngle_m + momentum[1] * sinAngle_m;
double ptheta = std::sqrt(newmomentum2 - std::pow(pr, 2));
double px = pr * cosAngle_m - ptheta * sinAngle_m ; // x
- double py = pr * sinAngle_m + ptheta * cosAngle_m; // y
+ double py = pr * sinAngle_m + ptheta * cosAngle_m; // y
double rotate = -derivate * (scale_m * Units::MVpm2Vpm) / (rmax_m - rmin_m) * std::sin(nphase) / (frequency * Physics::two_pi) / (betgam * restMass / Physics::c / chargenumber); // radian
@@ -422,12 +423,10 @@ void RFCavity::getMomentaKick(const double normalRadius,
if (PID == 0) {
Inform m("OPAL", *gmsg, Ippl::myNode());
-
m << "* Cavity " << getName() << " Phase= " << tempdegree << " [deg] transit time factor= " << Ufactor
- << " dE= " << dgam *restMass * Units::eV2MeV << " [MeV]"
- << " E_kin= " << (gamma - 1.0)*restMass * Units::eV2MeV << " [MeV] Time dep freq = " << frequencyTD_m->getValue(t) << endl;
+ << " dE= " << dgam * restMass * Units::eV2MeV << " [MeV]"
+ << " E_kin= " << (gamma - 1.0) * restMass * Units::eV2MeV << " [MeV] Time dep freq = " << frequencyTD_m->getValue(t) << endl;
}
-
}
/* cubic spline subrutine */
diff --git a/src/Classic/AbsBeamline/Ring.cpp b/src/Classic/AbsBeamline/Ring.cpp
index 8fbd3bbaa948393f360b11ea1a0a12707fb180fd..1b54b3b5b6255dfb3a99ec2ba4dc0410c85e9c75 100644
--- a/src/Classic/AbsBeamline/Ring.cpp
+++ b/src/Classic/AbsBeamline/Ring.cpp
@@ -1,43 +1,34 @@
-/*
- * Copyright (c) 2012-2023, Chris Rogers
- * All rights reserved.
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions are met:
- * 1. Redistributions of source code must retain the above copyright notice,
- * this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright notice,
- * this list of conditions and the following disclaimer in the documentation
- * and/or other materials provided with the distribution.
- * 3. Neither the name of STFC nor the names of its contributors may be used to
- * endorse or promote products derived from this software without specific
- * prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
- * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
- * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
- * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
- * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
- * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
- * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
- * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
- * POSSIBILITY OF SUCH DAMAGE.
- */
-
+//
+// Class Ring
+// Defines the abstract interface for a ring type geometry.
+//
+// Copyright (c) 2012 - 2023, Chris Rogers, RAL-UKRI, England, UK
+// All rights reserved
+//
+// This file is part of OPAL.
+//
+// OPAL is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// You should have received a copy of the GNU General Public License
+// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
+//
#include "AbsBeamline/Ring.h"
-#include <sstream>
-#include <limits>
-#include <cmath>
-
-#include "Utility/Inform.h" // ippl
+#include "Utility/Inform.h"
-#include "Fields/EMField.h"
#include "AbsBeamline/BeamlineVisitor.h"
+#include "Algorithms/PartBunchBase.h"
#include "BeamlineGeometry/Euclid3D.h"
+#include "Fields/EMField.h"
+#include "Physics/Units.h"
#include "Structure/LossDataSink.h"
-#include "Algorithms/PartBunchBase.h"
+
+#include <cmath>
+#include <limits>
+#include <sstream>
// fairly generous tolerance here... maybe too generous? Maybe should be
// user parameter?
@@ -101,16 +92,18 @@ Ring::~Ring() {
delete section_list_m[i];
}
-bool Ring::apply(const size_t &id, const double &t, Vector_t &E,
- Vector_t &B) {
- bool flagNeedUpdate =
- apply(refPartBunch_m->R[id], refPartBunch_m->P[id], t, E, B);
- if(flagNeedUpdate) {
+bool Ring::apply(const size_t &id, const double &t,
+ Vector_t &E, Vector_t &B) {
+
+ bool flagNeedUpdate = apply(refPartBunch_m->R[id], refPartBunch_m->P[id], t, E, B);
+
+ if (flagNeedUpdate) {
*gmsgALL << level4 << getName() << ": particle " << id
<< " at " << refPartBunch_m->R[id]
<< " m out of the field map boundary" << endl;
+
lossDS_m->addParticle(OpalParticle(id,
- refPartBunch_m->R[id] * Vector_t(1000.0), refPartBunch_m->P[id],
+ refPartBunch_m->R[id] , refPartBunch_m->P[id],
t,
refPartBunch_m->Q[id], refPartBunch_m->M[id]));
@@ -138,7 +131,7 @@ bool Ring::apply(const Vector_t &R, const Vector_t &/*P*/,
Vector_t B_temp(0.0, 0.0, 0.0);
Vector_t E_temp(0.0, 0.0, 0.0);
// Super-TEMP! cyclotron tracker now uses m internally, have to change to mm here to match old field limits -DW
- outOfBounds &= sections[i]->getFieldValue(R * Vector_t(1000.0), refPartBunch_m->get_centroid() * Vector_t(1000.0), t, E_temp, B_temp);
+ outOfBounds &= sections[i]->getFieldValue(R * Vector_t(Units::m2mm), refPartBunch_m->get_centroid() * Units::m2mm, t, E_temp, B_temp);
B += (scale_m * B_temp);
E += (scale_m * E_temp);
}
diff --git a/src/Classic/AbsBeamline/Ring.h b/src/Classic/AbsBeamline/Ring.h
index 6329a70577f3cd00c5267e0f80caf386e0d55675..9951c481ddc6520d8d0f5d8afabeef5f8a1046a6 100644
--- a/src/Classic/AbsBeamline/Ring.h
+++ b/src/Classic/AbsBeamline/Ring.h
@@ -1,30 +1,20 @@
-/*
- * Copyright (c) 2012-2014, Chris Rogers
- * All rights reserved.
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions are met:
- * 1. Redistributions of source code must retain the above copyright notice,
- * this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright notice,
- * this list of conditions and the following disclaimer in the documentation
- * and/or other materials provided with the distribution.
- * 3. Neither the name of STFC nor the names of its contributors may be used to
- * endorse or promote products derived from this software without specific
- * prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
- * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
- * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
- * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
- * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
- * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
- * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
- * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
- * POSSIBILITY OF SUCH DAMAGE.
- */
-
+//
+// Class Ring
+// Defines the abstract interface for a ring type geometry.
+//
+// Copyright (c) 2012 - 2023, Chris Rogers, RAL-UKRI, England, UK
+// All rights reserved
+//
+// This file is part of OPAL.
+//
+// OPAL is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// You should have received a copy of the GNU General Public License
+// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
+//
#ifndef RING_H
#define RING_H
diff --git a/src/Distribution/ClosedOrbitFinder.h b/src/Distribution/ClosedOrbitFinder.h
index dddec5457dca63c9f151a320956d47505c73e29d..5f900a15d171f9efedc1e79e554c9d4b512c605e 100644
--- a/src/Distribution/ClosedOrbitFinder.h
+++ b/src/Distribution/ClosedOrbitFinder.h
@@ -37,15 +37,13 @@
#include <utility>
#include <vector>
-#include "Utilities/Options.h"
+#include "AbsBeamline/Cyclotron.h"
+#include "AbstractObjects/OpalData.h"
#include "Utilities/OpalException.h"
+#include "Utilities/Options.h"
#include "Physics/Physics.h"
#include "Physics/Units.h"
-#include "AbstractObjects/OpalData.h"
-
-#include "AbsBeamline/Cyclotron.h"
-
// include headers for integration
#include <boost/numeric/odeint/integrate/integrate_n_steps.hpp>
#include <boost/filesystem.hpp>
@@ -53,205 +51,205 @@
extern Inform *gmsg;
template<typename Value_type, typename Size_type, class Stepper>
-class ClosedOrbitFinder
-{
- public:
- /// Type of variables
- typedef Value_type value_type;
- /// Type for specifying sizes
- typedef Size_type size_type;
- /// Type of container for storing quantities (path length, orbit, etc.)
- typedef std::vector<value_type> container_type;
- /// Type for holding state of ODE values
- typedef std::vector<value_type> state_type;
-
- typedef std::function<void(const state_type&, state_type&, const double)> function_t;
-
- /// Sets the initial values for the integration and calls findOrbit().
- /*!
- * @param E0 is the potential energy (particle energy at rest) [MeV].
- * @param q is the particle charge [e]
- * @param N specifies the number of splits (2pi/N), i.e number of integration steps
- * @param cycl is the cyclotron element
- * @param domain is a boolean (default: true). If "true" the closed orbit is computed over a single sector,
- * otherwise over 2*pi.
- * @param Nsectors is an int (default: 1). Number of sectors that the field map is averaged over
- * in order to avoid first harmonic. Only valid if domain is false
- */
- ClosedOrbitFinder(value_type E0, value_type q, size_type N,
- Cyclotron* cycl, bool domain = true, int Nsectors = 1);
+class ClosedOrbitFinder {
+
+public:
+ /// Type of variables
+ typedef Value_type value_type;
+ /// Type for specifying sizes
+ typedef Size_type size_type;
+ /// Type of container for storing quantities (path length, orbit, etc.)
+ typedef std::vector<value_type> container_type;
+ /// Type for holding state of ODE values
+ typedef std::vector<value_type> state_type;
+
+ typedef std::function<void(const state_type&, state_type&, const double)> function_t;
+
+ /// Sets the initial values for the integration and calls findOrbit().
+ /*!
+ * @param E0 is the potential energy (particle energy at rest) [MeV].
+ * @param q is the particle charge [e]
+ * @param N specifies the number of splits (2pi/N), i.e number of integration steps
+ * @param cycl is the cyclotron element
+ * @param domain is a boolean (default: true). If "true" the closed orbit is computed over a single sector,
+ * otherwise over 2*pi.
+ * @param Nsectors is an int (default: 1). Number of sectors that the field map is averaged over
+ * in order to avoid first harmonic. Only valid if domain is false
+ */
+ ClosedOrbitFinder(value_type E0, value_type q, size_type N,
+ Cyclotron* cycl, bool domain = true, int Nsectors = 1);
- /// Returns the inverse bending radius (size of container N+1)
- container_type getInverseBendingRadius(const value_type& angle = 0);
+ /// Returns the inverse bending radius (size of container N+1)
+ container_type getInverseBendingRadius(const value_type& angle = 0);
- /// Returns the step lengths of the path (size of container N+1)
- container_type getPathLength(const value_type& angle = 0);
+ /// Returns the step lengths of the path (size of container N+1)
+ container_type getPathLength(const value_type& angle = 0);
- /// Returns the field index (size of container N+1)
- container_type getFieldIndex(const value_type& angle = 0);
+ /// Returns the field index (size of container N+1)
+ container_type getFieldIndex(const value_type& angle = 0);
- /// Returns the radial and vertical tunes (in that order)
- std::pair<value_type,value_type> getTunes();
+ /// Returns the radial and vertical tunes (in that order)
+ std::pair<value_type,value_type> getTunes();
- /// Returns the closed orbit (size of container N+1) starting at specific angle (only makes sense when computing
- /// the closed orbit for a whole turn) (default value: 0°).
- /// Attention: It computes the starting index of the array. If it's not an integer it just cuts the floating point
- /// part, i.e. it takes the next starting index below. There's no interpolation of the radius.
- /*!
- * @param angle is the start angle for the output. Has to be within [0°,360°[ (default: 0°).
- */
- container_type getOrbit(value_type angle = 0);
-
- /// Returns the momentum of the orbit (size of container N+1)starting at specific angle (only makes sense when
- /// computing the closed orbit for a whole turn) (default value: 0°), \f$ \left[ p_{r} \right] = \si{m}\f$.
- /// Attention: It computes the starting index of the array. If it's not an integer it just cuts the floating point
- /// part, i.e. it takes the next starting index below. There's no interpolation of the momentum.
- /*!
- * @param angle is the start angle for the output. Has to be within [0°,360°[ (default: 0°).
- * @returns the momentum in \f$ \beta * \gamma \f$ units
- */
- container_type getMomentum(value_type angle = 0);
+ /// Returns the closed orbit (size of container N+1) starting at specific angle (only makes sense when computing
+ /// the closed orbit for a whole turn) (default value: 0°).
+ /// Attention: It computes the starting index of the array. If it's not an integer it just cuts the floating point
+ /// part, i.e. it takes the next starting index below. There's no interpolation of the radius.
+ /*!
+ * @param angle is the start angle for the output. Has to be within [0°,360°[ (default: 0°).
+ */
+ container_type getOrbit(value_type angle = 0);
+
+ /// Returns the momentum of the orbit (size of container N+1)starting at specific angle (only makes sense when
+ /// computing the closed orbit for a whole turn) (default value: 0°), \f$ \left[ p_{r} \right] = \si{m}\f$.
+ /// Attention: It computes the starting index of the array. If it's not an integer it just cuts the floating point
+ /// part, i.e. it takes the next starting index below. There's no interpolation of the momentum.
+ /*!
+ * @param angle is the start angle for the output. Has to be within [0°,360°[ (default: 0°).
+ * @returns the momentum in \f$ \beta * \gamma \f$ units
+ */
+ container_type getMomentum(value_type angle = 0);
- /// Returns the average orbit radius
- value_type getAverageRadius();
+ /// Returns the average orbit radius
+ value_type getAverageRadius();
- /// Returns the frequency error
- value_type getFrequencyError();
+ /// Returns the frequency error
+ value_type getFrequencyError();
- /// Returns true if a closed orbit could be found
- bool isConverged();
+ /// Returns true if a closed orbit could be found
+ bool isConverged();
- /// Computes the closed orbit
- /*!
- * @param accuracy specifies the accuracy of the closed orbit
- * @param maxit is the maximal number of iterations done for finding the closed orbit
- * @param ekin energy for which to find closed orbit (in tune mode: upper limit of range)
- * @param dE step increase [MeV]
- * @param rguess initial radius guess in [mm]
- * @param isTuneMode comptute tunes of all energies in one sweep
- */
- bool findOrbit(value_type accuracy, size_type maxit,
- value_type ekin,
- value_type dE = 1.0, value_type rguess = -1.0,
- bool isTuneMode = false);
-
- /// Fills in the values of h_m, ds_m, fidx_m.
- void computeOrbitProperties(const value_type& E);
-
- private:
- /// This function is called by the function getTunes().
- /*! Transfer matrix Y = [y11, y12; y21, y22] (see Gordon paper for more details).
- * @param y are the positions (elements y11 and y12 of Y)
- * @param py2 is the momentum of the second solution (element y22 of Y)
- * @param ncross is the number of sign changes (\#crossings of zero-line)
- */
- value_type computeTune(const std::array<value_type,2>&, value_type, size_type);
+ /// Computes the closed orbit
+ /*!
+ * @param accuracy specifies the accuracy of the closed orbit
+ * @param maxit is the maximal number of iterations done for finding the closed orbit
+ * @param ekin energy for which to find closed orbit (in tune mode: upper limit of range)
+ * @param dE step increase [MeV]
+ * @param rguess initial radius guess in [m]
+ * @param isTuneMode comptute tunes of all energies in one sweep
+ */
+ bool findOrbit(value_type accuracy, size_type maxit,
+ value_type ekin,
+ value_type dE = 1.0, value_type rguess = -1.0,
+ bool isTuneMode = false);
+
+ /// Fills in the values of h_m, ds_m, fidx_m.
+ void computeOrbitProperties(const value_type& E);
+
+private:
+ /// This function is called by the function getTunes().
+ /*! Transfer matrix Y = [y11, y12; y21, y22] (see Gordon paper for more details).
+ * @param y are the positions (elements y11 and y12 of Y)
+ * @param py2 is the momentum of the second solution (element y22 of Y)
+ * @param ncross is the number of sign changes (\#crossings of zero-line)
+ */
+ value_type computeTune(const std::array<value_type,2>&, value_type, size_type);
- // Compute closed orbit for given energy
- bool findOrbitOfEnergy_m(const value_type&, container_type&, value_type&,
- const value_type&, size_type);
+ // Compute closed orbit for given energy
+ bool findOrbitOfEnergy_m(const value_type&, container_type&, value_type&,
+ const value_type&, size_type);
- /// This function computes nzcross_ which is used to compute the tune in z-direction and the frequency error
+ /// This function computes nzcross_ which is used to compute the tune in z-direction and the frequency error
// void computeVerticalOscillations();
- /// This function rotates the calculated closed orbit finder properties to the initial angle
- container_type rotate(value_type angle, const container_type& orbitProperty);
-
- /// Stores current position in horizontal direction for the solutions of the ODE with different initial values
- std::array<value_type,2> x_m; // x_m = [x1, x2]
- /// Stores current momenta in horizontal direction for the solutions of the ODE with different initial values
- std::array<value_type,2> px_m; // px_m = [px1, px2]
- /// Stores current position in vertical direction for the solutions of the ODE with different initial values
- std::array<value_type,2> z_m; // z_m = [z1, z2]
- /// Stores current momenta in vertical direction for the solutions of the ODE with different initial values
- std::array<value_type,2> pz_m; // pz_m = [pz1, pz2]
-
- /// Stores the inverse bending radius
- container_type h_m;
- /// Stores the step length
- container_type ds_m;
- /// Stores the radial orbit (size: N_m+1)
- container_type r_m;
- /// Stores the vertical oribt (size: N_m+1)
- container_type vz_m;
- /// Stores the radial momentum
- container_type pr_m;
- /// Stores the vertical momentum
- container_type vpz_m;
- /// Stores the field index
- container_type fidx_m;
-
- /// Counts the number of zero-line crossings in horizontal direction (used for computing horizontal tune)
- size_type nxcross_m;
- /// Counts the number of zero-line crossings in vertical direction (used for computing vertical tune)
- size_type nzcross_m; //#crossings of zero-line in x- and z-direction
-
- /// Is the rest mass [MeV / c**2]
- value_type E0_m;
-
- // Is the particle charge [e]
- value_type q_m;
-
- /// Is the nominal orbital frequency
- /* (see paper of Dr. C. Baumgarten: "Transverse-Longitudinal
- * Coupling by Space Charge in Cyclotrons" (2012), formula (1))
- */
- value_type wo_m;
- /// Number of integration steps
- size_type N_m;
- /// Is the angle step size
- value_type dtheta_m;
-
- /// Is the average radius
- value_type ravg_m;
-
- /// Is the phase
- value_type phase_m;
-
- /**
- * Stores the last orbit value (since we have to return to the beginning to check the convergence in the
- * findOrbit() function. This last value is then deleted from the array but is stored in lastOrbitVal_m to
- * compute the vertical oscillations)
- */
- /* value_type lastOrbitVal_m; */
+ /// This function rotates the calculated closed orbit finder properties to the initial angle
+ container_type rotate(value_type angle, const container_type& orbitProperty);
+
+ /// Stores current position in horizontal direction for the solutions of the ODE with different initial values
+ std::array<value_type,2> x_m; // x_m = [x1, x2]
+ /// Stores current momenta in horizontal direction for the solutions of the ODE with different initial values
+ std::array<value_type,2> px_m; // px_m = [px1, px2]
+ /// Stores current position in vertical direction for the solutions of the ODE with different initial values
+ std::array<value_type,2> z_m; // z_m = [z1, z2]
+ /// Stores current momenta in vertical direction for the solutions of the ODE with different initial values
+ std::array<value_type,2> pz_m; // pz_m = [pz1, pz2]
+
+ /// Stores the inverse bending radius
+ container_type h_m;
+ /// Stores the step length
+ container_type ds_m;
+ /// Stores the radial orbit (size: N_m+1)
+ container_type r_m;
+ /// Stores the vertical oribt (size: N_m+1)
+ container_type vz_m;
+ /// Stores the radial momentum
+ container_type pr_m;
+ /// Stores the vertical momentum
+ container_type vpz_m;
+ /// Stores the field index
+ container_type fidx_m;
+
+ /// Counts the number of zero-line crossings in horizontal direction (used for computing horizontal tune)
+ size_type nxcross_m;
+ /// Counts the number of zero-line crossings in vertical direction (used for computing vertical tune)
+ size_type nzcross_m; //#crossings of zero-line in x- and z-direction
+
+ /// Is the rest mass [MeV / c**2]
+ value_type E0_m;
+
+ // Is the particle charge [e]
+ value_type q_m;
+
+ /// Is the nominal orbital frequency
+ /* (see paper of Dr. C. Baumgarten: "Transverse-Longitudinal
+ * Coupling by Space Charge in Cyclotrons" (2012), formula (1))
+ */
+ value_type wo_m;
+ /// Number of integration steps
+ size_type N_m;
+ /// Is the angle step size
+ value_type dtheta_m;
+
+ /// Is the average radius
+ value_type ravg_m;
+
+ /// Is the phase
+ value_type phase_m;
+
+ /**
+ * Stores the last orbit value (since we have to return to the beginning to check the convergence in the
+ * findOrbit() function. This last value is then deleted from the array but is stored in lastOrbitVal_m to
+ * compute the vertical oscillations)
+ */
+ /* value_type lastOrbitVal_m; */
- /**
- * Stores the last momentum value (since we have to return to the beginning to check the convergence in the
- * findOrbit() function. This last value is then deleted from the array but is stored in lastMomentumVal_m to
- * compute the vertical oscillations)
- */
- /* value_type lastMomentumVal_m; */
+ /**
+ * Stores the last momentum value (since we have to return to the beginning to check the convergence in the
+ * findOrbit() function. This last value is then deleted from the array but is stored in lastMomentumVal_m to
+ * compute the vertical oscillations)
+ */
+ /* value_type lastMomentumVal_m; */
- /**
- * Boolean which is true by default. "true": orbit integration over one sector only, "false": integration
- * over 2*pi
- */
- bool domain_m;
- /**
- * Number of sectors to average the field map over
- * in order to avoid first harmonic. Only valid if domain is false
- */
- int nSectors_m;
+ /**
+ * Boolean which is true by default. "true": orbit integration over one sector only, "false": integration
+ * over 2*pi
+ */
+ bool domain_m;
+ /**
+ * Number of sectors to average the field map over
+ * in order to avoid first harmonic. Only valid if domain is false
+ */
+ int nSectors_m;
- /// Defines the stepper for integration of the ODE's
- Stepper stepper_m;
+ /// Defines the stepper for integration of the ODE's
+ Stepper stepper_m;
- /*!
- * This quantity is defined in the paper "Transverse-Longitudinal Coupling by Space Charge in Cyclotrons"
- * of Dr. Christian Baumgarten (2012)
- * The lambda function takes the orbital frequency \f$ \omega_{o} \f$ (also defined in paper) as input argument.
- */
- std::function<double(double)> acon_m = [](double wo) { return Physics::c / wo; };
+ /*!
+ * This quantity is defined in the paper "Transverse-Longitudinal Coupling by Space Charge in Cyclotrons"
+ * of Dr. Christian Baumgarten (2012)
+ * The lambda function takes the orbital frequency \f$ \omega_{o} \f$ (also defined in paper) as input argument.
+ */
+ std::function<double(double)> acon_m = [](double wo) { return Physics::c / wo; };
- /// Cyclotron unit \f$ \left[T\right] \f$ (Tesla)
- /*!
- * The lambda function takes the orbital frequency \f$ \omega_{o} \f$ as input argument.
- */
- std::function<double(double, double)> bcon_m = [this](double e0, double wo) {
- return e0 * 1.0e7 / (q_m * Physics::c * Physics::c / wo);
- };
+ /// Cyclotron unit \f$ \left[T\right] \f$ (Tesla)
+ /*!
+ * The lambda function takes the orbital frequency \f$ \omega_{o} \f$ as input argument.
+ */
+ std::function<double(double, double)> bcon_m = [this](double e0, double wo) {
+ return e0 * 1.0e7 / (q_m * Physics::c * Physics::c / wo);
+ };
- Cyclotron* cycl_m;
+ Cyclotron* cycl_m;
};
// -----------------------------------------------------------------------------------------------------------------------
@@ -313,28 +311,31 @@ template<typename Value_type, typename Size_type, class Stepper>
inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::container_type
ClosedOrbitFinder<Value_type, Size_type, Stepper>::getInverseBendingRadius(const value_type& angle)
{
- if (angle != 0.0)
+ if (angle != 0.0) {
return rotate(angle, h_m);
- else
+ } else {
return h_m;
+ }
}
template<typename Value_type, typename Size_type, class Stepper>
inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::container_type
ClosedOrbitFinder<Value_type, Size_type, Stepper>::getPathLength(const value_type& angle)
{
- if (angle != 0.0)
+ if (angle != 0.0) {
return rotate(angle, ds_m);
- else
+ } else {
return ds_m;
+ }
}
template<typename Value_type, typename Size_type, class Stepper>
inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::container_type
ClosedOrbitFinder<Value_type, Size_type, Stepper>::getFieldIndex(const value_type& angle)
{
- if (angle != 0.0)
+ if (angle != 0.0) {
return rotate(angle, fidx_m);
+ }
return fidx_m;
}
@@ -353,10 +354,11 @@ template<typename Value_type, typename Size_type, class Stepper>
inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::container_type
ClosedOrbitFinder<Value_type, Size_type, Stepper>::getOrbit(value_type angle)
{
- if (angle != 0.0)
+ if (angle != 0.0) {
return rotate(angle, r_m);
- else
+ } else {
return r_m;
+ }
}
template<typename Value_type, typename Size_type, class Stepper>
@@ -365,9 +367,9 @@ inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::container_typ
{
container_type pr = pr_m;
- if (angle != 0.0)
+ if (angle != 0.0) {
pr = rotate(angle, pr);
-
+ }
// change units from meters to \beta * \gamma
/* in Gordon paper:
*
@@ -409,8 +411,6 @@ typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::value_type
// PRIVATE MEMBER FUNCTIONS
// -----------------------------------------------------------------------------------------------------------------------
-
-
template<typename Value_type, typename Size_type, class Stepper>
bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type accuracy,
size_type maxit,
@@ -482,8 +482,8 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
// iterate until suggested energy (start with minimum energy)
// increase energy by dE
- *gmsg << level3 << "Start iteration to find closed orbit of energy " << E_fin << " MeV "
- << "in steps of " << dE << " MeV." << endl;
+ *gmsg << level3 << "Start iteration to find closed orbit of energy "
+ << E_fin << " MeV " << "in steps of " << dE << " MeV." << endl;
for (; E <= E_fin + eps; E += dE) {
@@ -503,7 +503,7 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
// r pr z pz
init = {beta * acon, 0.0, 0.0, 1.0};
if (rguess >= 0.0) {
- init[0] = rguess * 0.001;
+ init[0] = rguess;
}
guess = FIRST;
} else if (guess == FIRST) {
@@ -545,11 +545,10 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
pn1 = pr_m[0] / p;
if ( isTuneMode ) {
-
this->computeOrbitProperties(E);
value_type reo = this->getOrbit( cycl_m->getPHIinit())[0];
value_type peo = this->getMomentum(cycl_m->getPHIinit())[0];
- std::pair<value_type , value_type > tunes = this->getTunes();
+ std::pair<value_type, value_type> tunes = this->getTunes();
*gmsg << std::left
<< "* ----------------------------" << endl
@@ -591,8 +590,9 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
/* domain_m = true --> only integrated over a single sector
* --> multiply by cycl_m->getSymmetry() to get correct phase_m
*/
- if (domain_m)
+ if (domain_m) {
phase_m *= cycl_m->getSymmetry();
+ }
// returns true if converged, otherwise false
return error < accuracy;
@@ -618,8 +618,7 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbitOfEnergy_m(
// index for reaching next element of the arrays r and pr (no nicer way found yet)
size_type idx = 0;
// observer for storing the current value after each ODE step (e.g. Runge-Kutta step) into the containers of r and pr
- auto store = [&](state_type& y, const value_type /*t*/)
- {
+ auto store = [&](state_type& y, const value_type /*t*/) {
r_m[idx] = y[0];
pr_m[idx] = y[1];
vz_m[idx] = y[6];
@@ -649,16 +648,16 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbitOfEnergy_m(
size_type niterations = 0;
// energy dependent values
- value_type en = E / E0_m; // en = E/E0 = E/(mc^2) (E0 is potential energy)
+ value_type en = E / E0_m; // en = E/E0 = E/(mc^2) (E0 is potential energy)
value_type gamma = en + 1.0;
- value_type p = acon_m(wo_m) * std::sqrt(en * (2.0 + en)); // momentum [p] = m; Gordon, formula (3)
+ value_type p = acon_m(wo_m) * std::sqrt(en * (2.0 + en)); // momentum [p] = m; Gordon, formula (3)
value_type gamma2 = gamma * gamma; // = gamma^2
value_type invgamma4 = 1.0 / (gamma2 * gamma2); // = 1/gamma^4
- value_type p2 = p * p; // p^2 = p*p
+ value_type p2 = p * p; // p^2 = p*p
// helper constants
- value_type pr2; // squared radial momentum (pr^2 = pr*pr)
- value_type ptheta, invptheta; // azimuthal momentum
+ value_type pr2; // squared radial momentum (pr^2 = pr*pr)
+ value_type ptheta, invptheta; // azimuthal momentum
// define the six ODEs (using lambda function)
function_t orbit_integration = [&](const state_type &y,
@@ -779,9 +778,11 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbitOfEnergy_m(
} while ((error > accuracy) && (niterations++ < maxit));
- if (error > accuracy)
- *gmsg << "findOrbit not converged after " << niterations << " iterations with error: " << error << ". Needed accuracy " << accuracy << endl;
-
+ if (error > accuracy) {
+ *gmsg << "findOrbit not converged after " << niterations
+ << " iterations with error: " << error
+ << ". Needed accuracy " << accuracy << endl;
+ }
return (error < accuracy);
}
@@ -805,14 +806,14 @@ Value_type ClosedOrbitFinder<Value_type, Size_type, Stepper>::computeTune(const
value_type muPrime;
if (abscos > 1.0) {
// store the number of crossings
- if (uneven)
+ if (uneven) {
ncross = ncross + 1;
-
+ }
// Gordon, formula (36b)
- muPrime = -std::acosh(abscos); // mu'
+ muPrime = -std::acosh(abscos); // mu'
} else {
- muPrime = (uneven) ? std::acos(-cos) : std::acos(cos); // mu'
+ muPrime = (uneven) ? std::acos(-cos) : std::acos(cos); // mu'
/* It has to be fulfilled: 0<= mu' <= pi
* But since |cos(mu)| <= 1, we have
* -1 <= cos(mu) <= 1 --> 0 <= mu <= pi (using above programmed line), such
@@ -837,9 +838,9 @@ Value_type ClosedOrbitFinder<Value_type, Size_type, Stepper>::computeTune(const
/* domain_m = true --> only integrated over a single sector --> multiply by cycl_m->getSymmetry() to
* get correct tune.
*/
- if (domain_m)
+ if (domain_m) {
mu *= cycl_m->getSymmetry();
-
+ }
return mu * Physics::u_two_pi;
}
@@ -903,15 +904,12 @@ ClosedOrbitFinder<Value_type, Size_type, Stepper>::rotate(value_type angle, cons
container_type orbitPropertyCopy = orbitProperty;
- // compute the number of steps per degree
- value_type deg_step = N_m / 360.0;
-
if (angle < 0.0) {
- angle = 360.0 + angle;
+ angle = Physics::two_pi + angle;
}
// compute starting point
- unsigned int start = deg_step * angle;
+ unsigned int start = angle / dtheta_m;
start %= orbitProperty.size();
@@ -922,7 +920,6 @@ ClosedOrbitFinder<Value_type, Size_type, Stepper>::rotate(value_type angle, cons
std::copy_n(orbitProperty.begin(), start, orbitPropertyCopy.end() - start);
return orbitPropertyCopy;
-
}
#endif
diff --git a/src/Distribution/Distribution.cpp b/src/Distribution/Distribution.cpp
index 478f64eb760aaa12d234930f859911df38c31993..b5286c0db6599661345a9c7dd9972c1be54765ae 100644
--- a/src/Distribution/Distribution.cpp
+++ b/src/Distribution/Distribution.cpp
@@ -1226,7 +1226,7 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles,
*gmsg << "* NSTEPS = " << Nint << endl
<< "* HN = " << CyclotronElement->getCyclHarm()
- << " PHIINIT = " << CyclotronElement->getPHIinit() << endl
+ << " PHIINIT = " << CyclotronElement->getPHIinit() * Units::rad2deg << endl
<< "* FIELD MAP = " << CyclotronElement->getFieldMapFN() << endl
<< "* ----------------------------------------------------" << endl;
@@ -1260,8 +1260,8 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles,
typedef boost::numeric::odeint::runge_kutta4<container_t> rk4_t;
typedef ClosedOrbitFinder<double,unsigned int, rk4_t> cof_t;
- cof_t cof(massIneV*Units::eV2MeV, charge, Nint, CyclotronElement, full, Nsectors);
- cof.findOrbit(accuracy, maxitCOF, E_m*Units::eV2MeV, denergy, rguess, true);
+ cof_t cof(massIneV * Units::eV2MeV, charge, Nint, CyclotronElement, full, Nsectors);
+ cof.findOrbit(accuracy, maxitCOF, E_m * Units::eV2MeV, denergy, rguess, true);
throw EarlyLeaveException("Distribution::createMatchedGaussDistribution()",
"Do only tune calculation.");
@@ -1271,11 +1271,11 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles,
std::unique_ptr<SigmaGenerator> siggen = std::unique_ptr<SigmaGenerator>(
new SigmaGenerator(I_m,
- Attributes::getReal(itsAttr[Attrib::Distribution::EX])*Units::m2mm * Units::rad2mrad,
- Attributes::getReal(itsAttr[Attrib::Distribution::EY])*Units::m2mm * Units::rad2mrad,
- Attributes::getReal(itsAttr[Attrib::Distribution::ET])*Units::m2mm * Units::rad2mrad,
- E_m*Units::eV2MeV,
- massIneV*Units::eV2MeV,
+ Attributes::getReal(itsAttr[Attrib::Distribution::EX]) * Units::m2mm * Units::rad2mrad,
+ Attributes::getReal(itsAttr[Attrib::Distribution::EY]) * Units::m2mm * Units::rad2mrad,
+ Attributes::getReal(itsAttr[Attrib::Distribution::ET]) * Units::m2mm * Units::rad2mrad,
+ E_m * Units::eV2MeV,
+ massIneV * Units::eV2MeV,
charge,
CyclotronElement,
Nint,
@@ -1301,15 +1301,14 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles,
*gmsg << "* Sigma-Matrix " << endl;
for (unsigned int i = 0; i < siggen->getSigma().size1(); ++ i) {
- *gmsg << std::setprecision(4) << std::setw(11) << siggen->getSigma()(i,0);
+ *gmsg << std::setprecision(4) << std::setw(11) << siggen->getSigma()(i,0);
for (unsigned int j = 1; j < siggen->getSigma().size2(); ++ j) {
if (std::abs(siggen->getSigma()(i,j)) < 1.0e-15) {
- *gmsg << " & " << std::setprecision(4) << std::setw(11) << 0.0;
+ *gmsg << " & " << std::setprecision(4) << std::setw(11) << 0.0;
}
else{
- *gmsg << " & " << std::setprecision(4) << std::setw(11) << siggen->getSigma()(i,j);
+ *gmsg << " & " << std::setprecision(4) << std::setw(11) << siggen->getSigma()(i,j);
}
-
}
*gmsg << " \\\\" << endl;
}
@@ -1317,11 +1316,11 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles,
generateMatchedGauss(siggen->getSigma(), numberOfParticles, massIneV);
// update injection radius and radial momentum
- CyclotronElement->setRinit(siggen->getInjectionRadius() * Units::m2mm);
+ CyclotronElement->setRinit(siggen->getInjectionRadius());
CyclotronElement->setPRinit(siggen->getInjectionMomentum());
}
else {
- *gmsg << "* Not converged for " << E_m*Units::eV2MeV << " MeV" << endl;
+ *gmsg << "* Not converged for " << E_m * Units::eV2MeV << " MeV" << endl;
throw OpalException("Distribution::CreateMatchedGaussDistribution",
"didn't find any matched distribution.");
diff --git a/src/Elements/OpalCyclotron.cpp b/src/Elements/OpalCyclotron.cpp
index 2e5165585d42a79daf2f0b5997d348f7e120fae0..d73b4008105ada8521402e7c4a66760c4ba841c0 100644
--- a/src/Elements/OpalCyclotron.cpp
+++ b/src/Elements/OpalCyclotron.cpp
@@ -142,26 +142,28 @@ void OpalCyclotron::update() {
double harmnum = Attributes::getReal(itsAttr[CYHARMON]);
double symmetry = Attributes::getReal(itsAttr[SYMMETRY]);
- double rinit = Attributes::getReal(itsAttr[RINIT]);
+ double bscale = Attributes::getReal(itsAttr[BSCALE]);
+
+ double rinit = Units::mm2m * Attributes::getReal(itsAttr[RINIT]);
+ double phiinit = Units::deg2rad * Attributes::getReal(itsAttr[PHIINIT]);
+ double zinit = Units::mm2m * Attributes::getReal(itsAttr[ZINIT]);
double prinit = Attributes::getReal(itsAttr[PRINIT]);
- double phiinit = Attributes::getReal(itsAttr[PHIINIT]);
- double zinit = Attributes::getReal(itsAttr[ZINIT]);
double pzinit = Attributes::getReal(itsAttr[PZINIT]);
- double bscale = Attributes::getReal(itsAttr[BSCALE]);
- double minz = Attributes::getReal(itsAttr[MINZ]);
- double maxz = Attributes::getReal(itsAttr[MAXZ]);
- double minr = Attributes::getReal(itsAttr[MINR]);
- double maxr = Attributes::getReal(itsAttr[MAXR]);
+ double minz = Units::mm2m * Attributes::getReal(itsAttr[MINZ]);
+ double maxz = Units::mm2m * Attributes::getReal(itsAttr[MAXZ]);
+ double minr = Units::mm2m * Attributes::getReal(itsAttr[MINR]);
+ double maxr = Units::mm2m * Attributes::getReal(itsAttr[MAXR]);
- double fmLowE = Attributes::getReal(itsAttr[FMLOWE]);
- double fmHighE = Attributes::getReal(itsAttr[FMHIGHE]);
+ double fmLowE = Units::GeV2MeV * Attributes::getReal(itsAttr[FMLOWE]);
+ double fmHighE = Units::GeV2MeV * Attributes::getReal(itsAttr[FMHIGHE]);
bool spiral_flag = Attributes::getBool(itsAttr[SPIRAL]);
- double trimCoilThreshold = Attributes::getReal(itsAttr[TRIMCOILTHRESHOLD]);
+ double trimCoilThreshold = Units::T2kG * Attributes::getReal(itsAttr[TRIMCOILTHRESHOLD]);
cycl->setFieldMapFN(fmapfm);
cycl->setSymmetry(symmetry);
+ cycl->setBScale(bscale);
cycl->setRinit(rinit);
cycl->setPRinit(prinit);
@@ -169,8 +171,6 @@ void OpalCyclotron::update() {
cycl->setZinit(zinit);
cycl->setPZinit(pzinit);
- cycl->setBScale(bscale);
-
cycl->setCyclotronType(type);
cycl->setCyclHarm(harmnum);
@@ -179,8 +179,8 @@ void OpalCyclotron::update() {
cycl->setMinZ(minz);
cycl->setMaxZ(maxz);
- cycl->setFMLowE(fmLowE * Units::GeV2MeV);
- cycl->setFMHighE(fmHighE * Units::GeV2MeV);
+ cycl->setFMLowE(fmLowE);
+ cycl->setFMHighE(fmHighE);
cycl->setSpiralFlag(spiral_flag);
cycl->setTrimCoilThreshold(trimCoilThreshold);
diff --git a/src/Elements/OpalRingDefinition.cpp b/src/Elements/OpalRingDefinition.cpp
index f555d879edc3675d8c2293b5293614eb660a3521..7e9af9ccfbd5a872be13fe4d7c226bb7828fea0b 100644
--- a/src/Elements/OpalRingDefinition.cpp
+++ b/src/Elements/OpalRingDefinition.cpp
@@ -1,77 +1,84 @@
-/*
- * Copyright (c) 2012, Chris Rogers
- * All rights reserved.
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions are met:
- * 1. Redistributions of source code must retain the above copyright notice,
- * this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright notice,
- * this list of conditions and the following disclaimer in the documentation
- * and/or other materials provided with the distribution.
- * 3. Neither the name of STFC nor the names of its contributors may be used to
- * endorse or promote products derived from this software without specific
- * prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
- * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
- * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
- * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
- * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
- * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
- * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
- * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
- * POSSIBILITY OF SUCH DAMAGE.
- */
-
+//
+// Class OpalRingDefinition
+// The Opal Ring element.
+//
+// Copyright (c) 2012 - 2023, Chris Rogers, RAL-UKRI, England, UK
+// All rights reserved
+//
+// This file is part of OPAL.
+//
+// OPAL is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// You should have received a copy of the GNU General Public License
+// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
+//
#include "Elements/OpalRingDefinition.h"
-#include <limits>
-
#include "AbsBeamline/Ring.h"
#include "Attributes/Attributes.h"
#include "Physics/Units.h"
+#include "Utilities/OpalException.h"
+
+#include <limits>
OpalRingDefinition::OpalRingDefinition() :
OpalElement(SIZE, "RINGDEFINITION",
"The \"RINGDEFINITION\" element defines basic ring parameters.") {
- itsAttr[HARMONIC_NUMBER] = Attributes::makeReal("HARMONIC_NUMBER",
- "The assumed harmonic number of the ring (i.e. number of bunches in the ring on a given turn) (default = 1).", 1.0);
- itsAttr[LAT_RINIT] = Attributes::makeReal("LAT_RINIT",
- "The initial radius of the first element to be placed in the ring [m].");
- itsAttr[LAT_PHIINIT] = Attributes::makeReal("LAT_PHIINIT", "The initial angle around the ring of the first element to be placed. [deg]");
- itsAttr[LAT_THETAINIT] = Attributes::makeReal("LAT_THETAINIT", "The angle relative to the tangent of the ring for the first element to be placed [deg].");
- itsAttr[BEAM_PHIINIT] = Attributes::makeReal("BEAM_PHIINIT",
- "The initial angle around the ring of the beam [deg].");
- itsAttr[BEAM_THETAINIT] = Attributes::makeReal("BEAM_THETAINIT",
- "Defines an angular offset of the beam relative to the tangent vector, in the x-y plane [deg] (default = 0).", 0.0);
- itsAttr[BEAM_PRINIT] = Attributes::makeReal("BEAM_PRINIT",
- "An initial pr momentum offset of the beam.");
- itsAttr[BEAM_RINIT] = Attributes::makeReal("BEAM_RINIT",
- "The initial radius of the beam [m].");
- itsAttr[SYMMETRY] = Attributes::makeReal("SYMMETRY",
- "The rotational symmetry of the lattice.", 1.0);
- itsAttr[SCALE] = Attributes::makeReal("SCALE", "Scale the fields by a multiplicative factor", 1.0);
+ itsAttr[HARMONIC_NUMBER] = Attributes::makeReal
+ ("HARMONIC_NUMBER", "The assumed harmonic number of the ring (i.e. number of bunches in the ring on a given turn) (default = 1).", 1.0);
+
+ itsAttr[LAT_RINIT] = Attributes::makeReal
+ ("LAT_RINIT", "The initial radius of the first element to be placed in the ring [m].");
+
+ itsAttr[LAT_PHIINIT] = Attributes::makeReal
+ ("LAT_PHIINIT", "The initial angle around the ring of the first element to be placed. [deg]");
+
+ itsAttr[LAT_THETAINIT] = Attributes::makeReal
+ ("LAT_THETAINIT", "The angle relative to the tangent of the ring for the first element to be placed [deg].");
+
+ itsAttr[BEAM_PHIINIT] = Attributes::makeReal
+ ("BEAM_PHIINIT", "The initial angle around the ring of the beam [deg].");
+
+ itsAttr[BEAM_THETAINIT] = Attributes::makeReal
+ ("BEAM_THETAINIT", "Defines an angular offset of the beam relative to the tangent vector, in the x-y plane [deg] (default = 0).", 0.0);
+
+ itsAttr[BEAM_PRINIT] = Attributes::makeReal
+ ("BEAM_PRINIT", "An initial pr momentum offset of the beam.");
+
+ itsAttr[BEAM_RINIT] = Attributes::makeReal
+ ("BEAM_RINIT", "The initial radius of the beam [m].");
+
+ itsAttr[SYMMETRY] = Attributes::makeReal
+ ("SYMMETRY", "The rotational symmetry of the lattice.", 1.0);
+
+ itsAttr[SCALE] = Attributes::makeReal
+ ("SCALE", "Scale the fields by a multiplicative factor", 1.0);
+
// should be in RF cavity definition; this comes from cyclotron definition,
// but not right
- itsAttr[RFFREQ] = Attributes::makeReal("RFFREQ",
- "The nominal RF frequency of the ring [MHz].");
+ itsAttr[RFFREQ] = Attributes::makeReal
+ ("RFFREQ", "The nominal RF frequency of the ring [MHz].");
+
// I see also makeBool, but dont know how it works; no registerBoolAttribute
- itsAttr[IS_CLOSED] = Attributes::makeBool("IS_CLOSED",
- "Set to 'false' to disable checking for closure of the ring");
- itsAttr[MIN_R] = Attributes::makeReal("MIN_R",
- "Minimum allowed radius during tracking [m]. If not defined, any radius is allowed. If MIN_R is defined, MAX_R must also be defined.");
- itsAttr[MAX_R] = Attributes::makeReal("MAX_R",
- "Maximum allowed radius during tracking [m]. If not defined, any radius is allowed. If MAX_R is defined, MIN_R must also be defined.");
+ itsAttr[IS_CLOSED] = Attributes::makeBool
+ ("IS_CLOSED", "Set to 'false' to disable checking for closure of the ring");
+
+ itsAttr[MIN_R] = Attributes::makeReal
+ ("MIN_R", "Minimum allowed radius during tracking [m]. If not defined, any radius is allowed. If MIN_R is defined, MAX_R must also be defined.");
+
+ itsAttr[MAX_R] = Attributes::makeReal
+ ("MAX_R", "Maximum allowed radius during tracking [m]. If not defined, any radius is allowed. If MAX_R is defined, MIN_R must also be defined.");
registerOwnership();
setElement(new Ring("RING"));
}
-OpalRingDefinition* OpalRingDefinition::clone(const std::string &name) {
+OpalRingDefinition* OpalRingDefinition::clone(const std::string& name) {
return new OpalRingDefinition(name, this);
}
@@ -79,7 +86,7 @@ void OpalRingDefinition::print(std::ostream& out) const {
OpalElement::print(out);
}
-OpalRingDefinition::OpalRingDefinition(const std::string &name, OpalRingDefinition *parent):
+OpalRingDefinition::OpalRingDefinition(const std::string& name, OpalRingDefinition* parent):
OpalElement(name, parent) {
setElement(new Ring(name));
}
@@ -87,12 +94,13 @@ OpalRingDefinition::OpalRingDefinition(const std::string &name, OpalRingDefiniti
OpalRingDefinition::~OpalRingDefinition() {}
void OpalRingDefinition::update() {
- Ring *ring = dynamic_cast<Ring*>(getElement());
- ring->setBeamPhiInit(Attributes::getReal(itsAttr[BEAM_PHIINIT]));
- ring->setBeamThetaInit(Attributes::getReal(itsAttr[BEAM_THETAINIT]));
+ Ring* ring = dynamic_cast<Ring*>(getElement());
+
+ ring->setBeamPhiInit(Attributes::getReal(itsAttr[BEAM_PHIINIT]) * Units::deg2rad);
+ ring->setBeamThetaInit(Attributes::getReal(itsAttr[BEAM_THETAINIT]) * Units::deg2rad);
ring->setBeamPRInit(Attributes::getReal(itsAttr[BEAM_PRINIT]));
- ring->setBeamRInit(Attributes::getReal(itsAttr[BEAM_RINIT])*Units::m2mm);
- ring->setLatticeRInit(Attributes::getReal(itsAttr[LAT_RINIT])*Units::m2mm);
+ ring->setBeamRInit(Attributes::getReal(itsAttr[BEAM_RINIT]));
+ ring->setLatticeRInit(Attributes::getReal(itsAttr[LAT_RINIT]) * Units::m2mm);
ring->setLatticePhiInit(Attributes::getReal(itsAttr[LAT_PHIINIT])*Units::deg2rad);
ring->setLatticeThetaInit(Attributes::getReal(itsAttr[LAT_THETAINIT])*Units::deg2rad);
@@ -102,22 +110,18 @@ void OpalRingDefinition::update() {
ring->setHarmonicNumber(Attributes::getReal(itsAttr[HARMONIC_NUMBER]));
ring->setRFFreq(Attributes::getReal(itsAttr[RFFREQ]));
ring->setIsClosed(Attributes::getBool(itsAttr[IS_CLOSED]));
- double minR = -1;
- double maxR = -1;
-
- if (itsAttr[MIN_R]) {
- minR = Attributes::getReal(itsAttr[MIN_R]);
- if (!itsAttr[MAX_R]) {
- throw (""); // EXCEPTION
- }
- }
- if (itsAttr[MAX_R]) {
- maxR = Attributes::getReal(itsAttr[MAX_R]);
- if (!itsAttr[MIN_R]) {
- throw (""); // EXCEPTION
- }
+
+ if (itsAttr[MIN_R] && itsAttr[MAX_R]) {
+ double minR = Attributes::getReal(itsAttr[MIN_R]);
+ double maxR = Attributes::getReal(itsAttr[MAX_R]);
ring->setRingAperture(minR, maxR);
+ } else if (itsAttr[MIN_R] && !itsAttr[MAX_R]) {
+ throw OpalException("OpalRingDefinition::update",
+ "If MIN_R is defined, MAX_R must also be defined.");
+ } else if (!itsAttr[MIN_R] && itsAttr[MAX_R]) {
+ throw OpalException("OpalRingDefinition::update",
+ "If MAX_R is defined, MIN_R must also be defined.");
}
setElement(ring);
-}
\ No newline at end of file
+}
diff --git a/src/Elements/OpalRingDefinition.h b/src/Elements/OpalRingDefinition.h
index 6dbb47491782f9c30278e15889061850c11bcd2d..0e219987b29281e77aebc5013851f7ad1ef2882b 100644
--- a/src/Elements/OpalRingDefinition.h
+++ b/src/Elements/OpalRingDefinition.h
@@ -1,30 +1,20 @@
-/*
- * Copyright (c) 2012, Chris Rogers
- * All rights reserved.
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions are met:
- * 1. Redistributions of source code must retain the above copyright notice,
- * this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright notice,
- * this list of conditions and the following disclaimer in the documentation
- * and/or other materials provided with the distribution.
- * 3. Neither the name of STFC nor the names of its contributors may be used to
- * endorse or promote products derived from this software without specific
- * prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
- * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
- * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
- * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
- * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
- * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
- * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
- * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
- * POSSIBILITY OF SUCH DAMAGE.
- */
-
+//
+// Class OpalRingDefinition
+// The Opal Ring element.
+//
+// Copyright (c) 2012 - 2023, Chris Rogers, RAL-UKRI, England, UK
+// All rights reserved
+//
+// This file is part of OPAL.
+//
+// OPAL is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// You should have received a copy of the GNU General Public License
+// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
+//
#ifndef OPAL_OpalRingDefinition_HH
#define OPAL_OpalRingDefinition_HH
@@ -39,7 +29,7 @@ class Ring;
*/
class OpalRingDefinition: public OpalElement {
- public:
+public:
/** Enumeration maps to UI parameters */
enum {
LAT_RINIT = COMMON,
@@ -66,20 +56,21 @@ class OpalRingDefinition: public OpalElement {
virtual ~OpalRingDefinition();
/** Inherited copy constructor */
- virtual OpalRingDefinition *clone(const std::string &name);
+ virtual OpalRingDefinition* clone(const std::string& name);
/** Receive parameters from the parser and hand them off to the Ring */
void update();
/** Calls print on the OpalElement */
- virtual void print(std::ostream &) const;
- private:
+ virtual void print(std::ostream&) const;
+
+private:
// Not implemented.
- OpalRingDefinition(const OpalRingDefinition &);
- void operator=(const OpalRingDefinition &);
+ OpalRingDefinition(const OpalRingDefinition&);
+ void operator=(const OpalRingDefinition&);
// Clone constructor.
- OpalRingDefinition(const std::string &name, OpalRingDefinition *parent);
+ OpalRingDefinition(const std::string& name, OpalRingDefinition* parent);
};
#endif // OPAL_OpalRingDefinition_HH
diff --git a/src/Steppers/BorisPusher.h b/src/Steppers/BorisPusher.h
index eb7d9b44702770563e6fd3e84d434a8356d837cb..1c08aecb54472d798be79c5fa734fcd90c30a7fb 100644
--- a/src/Steppers/BorisPusher.h
+++ b/src/Steppers/BorisPusher.h
@@ -113,7 +113,7 @@ inline void BorisPusher::kick(const Vector_t &/*R*/, Vector_t &P,
P += cross(P, t);
*/
- double const gamma = sqrt(1.0 + dot(P, P));
+ double const gamma = std::sqrt(1.0 + dot(P, P));
Vector_t const t = 0.5 * dt * charge * Physics::c * Physics::c / (gamma * mass) * Bf;
Vector_t const w = P + cross(P, t);
Vector_t const s = 2.0 / (1.0 + dot(t, t)) * t;
@@ -135,7 +135,7 @@ inline void BorisPusher::push(Vector_t &R, const Vector_t &P, const double &/* d
* R[i] += 0.5 * P[i] * recpgamma;
* \endcode
*/
- R += 0.5 * P / sqrt(1.0 + dot(P, P));
+ R += 0.5 * P / std::sqrt(1.0 + dot(P, P));
}
#endif
diff --git a/src/Steppers/LF2.h b/src/Steppers/LF2.h
index 12b0bb32f2a173e1b7a70df6b7f5c9673277bd6d..d5645c2579c9d7e7296799339546209d51926bfb 100644
--- a/src/Steppers/LF2.h
+++ b/src/Steppers/LF2.h
@@ -18,27 +18,25 @@
#ifndef LF2_H
#define LF2_H
-#include "Stepper.h"
+#include "Steppers/Stepper.h"
#include "Physics/Physics.h"
/// Leap-Frog 2nd order
template <typename FieldFunction, typename ... Arguments>
class LF2 : public Stepper<FieldFunction, Arguments...> {
-
+
public:
-
LF2(const FieldFunction& fieldfunc) : Stepper<FieldFunction, Arguments ...>(fieldfunc) { }
-
+
private:
bool doAdvance_m(PartBunchBase<double, 3>* bunch,
const size_t& i,
const double& t,
const double dt,
Arguments& ... args) const;
-
-
+
void push_m(Vector_t& R, const Vector_t& P, const double& h) const;
-
+
bool kick_m(PartBunchBase<double, 3>* bunch, const size_t& i,
const double& t, const double& h,
Arguments& ... args) const;
diff --git a/src/Steppers/LF2.hpp b/src/Steppers/LF2.hpp
index 81ecc6706e1d7b7680274462325b486f7f9e7d5b..2031fa6e30a9c558b34c5f4ca7ae7c6f70c70d0c 100644
--- a/src/Steppers/LF2.hpp
+++ b/src/Steppers/LF2.hpp
@@ -15,8 +15,8 @@
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
-#include "BorisPusher.h"
#include "Physics/Units.h"
+#include "Steppers/BorisPusher.h"
template <typename FieldFunction, typename ... Arguments>
bool LF2<FieldFunction, Arguments ...>::doAdvance_m(PartBunchBase<double, 3>* bunch,
@@ -28,12 +28,12 @@ bool LF2<FieldFunction, Arguments ...>::doAdvance_m(PartBunchBase<double, 3>* bu
bool flagNoDeletion = true;
// push for first LF2 half step
- push_m(bunch->R[i], bunch->P[i], 0.5 * dt * Units::ns2s);
+ push_m(bunch->R[i], bunch->P[i], 0.5 * dt);
- flagNoDeletion = kick_m(bunch, i, t, dt * Units::ns2s, args ...);
+ flagNoDeletion = kick_m(bunch, i, t, dt, args ...);
// push for second LF2 half step
- push_m(bunch->R[i], bunch->P[i], 0.5 * dt * Units::ns2s);
+ push_m(bunch->R[i], bunch->P[i], 0.5 * dt);
return flagNoDeletion;
}
@@ -43,7 +43,7 @@ template <typename FieldFunction, typename ... Arguments>
void LF2<FieldFunction, Arguments ...>::push_m(Vector_t& R, const Vector_t& P,
const double& h) const
{
- double const gamma = sqrt(1.0 + dot(P, P));
+ double const gamma = std::sqrt(1.0 + dot(P, P));
double const c_gamma = Physics::c / gamma;
Vector_t const v = P * c_gamma;
R += h * v;
@@ -74,4 +74,4 @@ bool LF2<FieldFunction, Arguments ...>::kick_m(PartBunchBase<double, 3>* bunch,
h, M, q);
return true;
-}
\ No newline at end of file
+}
diff --git a/src/Steppers/RK4.h b/src/Steppers/RK4.h
index 84aa01c861e58735b71c957f9f01b6a487364b49..320684ad8e2993753d3a395b1a4d6e60b9b30253 100644
--- a/src/Steppers/RK4.h
+++ b/src/Steppers/RK4.h
@@ -18,16 +18,15 @@
#ifndef RK4_H
#define RK4_H
-#include "Stepper.h"
#include "Physics/Physics.h"
#include "Physics/Units.h"
+#include "Steppers/Stepper.h"
/// 4-th order Runnge-Kutta stepper
template <typename FieldFunction, typename ... Arguments>
class RK4 : public Stepper<FieldFunction, Arguments...> {
-
+
public:
-
RK4(const FieldFunction& fieldfunc) : Stepper<FieldFunction, Arguments ...>(fieldfunc) { }
private:
@@ -36,7 +35,6 @@ private:
const double& t,
const double dt,
Arguments& ... args) const;
-
/**
*
*
@@ -53,15 +51,12 @@ private:
double* yp,
const size_t& i,
Arguments& ... args) const;
-
-
+
void copyTo(const Vector_t& R, const Vector_t& P, double* x) const;
-
+
void copyFrom(Vector_t& R, Vector_t& P, double* x) const;
-
+
const double mass_coeff = 1.0e9 * Units::GeV2kg; // from GeV/c^2 to basic unit: GV*C*s^2/m^2, (1.0e9 converts V*C*s^2/m^2 to GV*C*s^2/m^2)
- const double c_mmtns = Physics::c * Units::m2mm / Units::s2ns;
- const double c_mtns = Physics::c / Units::s2ns;
};
#include "RK4.hpp"
diff --git a/src/Steppers/RK4.hpp b/src/Steppers/RK4.hpp
index f9b689e34b8db709e9863ce8a834e770aa540a9e..47a97f159a0efdb9439a99f82aa16678e03a5385 100644
--- a/src/Steppers/RK4.hpp
+++ b/src/Steppers/RK4.hpp
@@ -29,11 +29,11 @@ bool RK4<FieldFunction, Arguments ...>::doAdvance_m(PartBunchBase<double, 3>* bu
// t Independent variable (usually time)
// dt Step size (usually time step)
// i index of particle
-
+
double x[6];
-
+
this->copyTo(bunch->R[i], bunch->P[i], &x[0]);
-
+
double deriv1[6];
double deriv2[6];
double deriv3[6];
@@ -41,10 +41,10 @@ bool RK4<FieldFunction, Arguments ...>::doAdvance_m(PartBunchBase<double, 3>* bu
double xtemp[6];
// Evaluate f1 = f(x,t).
-
bool outOfBound = derivate_m(bunch, x, t, deriv1 , i, args ...);
- if (outOfBound)
+ if (outOfBound) {
return false;
+ }
// Evaluate f2 = f( x+dt*f1/2, t+dt/2 ).
const double half_dt = 0.5 * dt;
@@ -54,30 +54,35 @@ bool RK4<FieldFunction, Arguments ...>::doAdvance_m(PartBunchBase<double, 3>* bu
xtemp[j] = x[j] + half_dt * deriv1[j];
outOfBound = derivate_m(bunch, xtemp, t_half, deriv2 , i, args ...);
- if (outOfBound)
+ if (outOfBound) {
return false;
+ }
// Evaluate f3 = f( x+dt*f2/2, t+dt/2 ).
- for(int j = 0; j < 6; ++j)
+ for (int j = 0; j < 6; ++j) {
xtemp[j] = x[j] + half_dt * deriv2[j];
+ }
outOfBound = derivate_m(bunch, xtemp, t_half, deriv3 , i, args ...);
- if (outOfBound)
+ if (outOfBound) {
return false;
+ }
// Evaluate f4 = f( x+dt*f3, t+dt ).
double t_full = t + dt;
- for(int j = 0; j < 6; ++j)
+ for (int j = 0; j < 6; ++j) {
xtemp[j] = x[j] + dt * deriv3[j];
+ }
outOfBound = derivate_m(bunch, xtemp, t_full, deriv4 , i, args ...);
- if (outOfBound)
+ if (outOfBound) {
return false;
+ }
// Return x(t+dt) computed from fourth-order R-K.
- for(int j = 0; j < 6; ++j)
+ for (int j = 0; j < 6; ++j) {
x[j] += dt / 6.*(deriv1[j] + deriv4[j] + 2.*(deriv2[j] + deriv3[j]));
-
+ }
this->copyFrom(bunch->R[i], bunch->P[i], &x[0]);
return true;
@@ -100,31 +105,29 @@ bool RK4<FieldFunction, Arguments ...>::derivate_m(PartBunchBase<double, 3>* bun
externalB = Vector_t(0.0, 0.0, 0.0);
externalE = Vector_t(0.0, 0.0, 0.0);
- for(int j = 0; j < 3; ++j)
+ for (int j = 0; j < 3; ++j) {
tempR(j) = y[j];
-
+ }
bunch->R[i] = tempR;
bool outOfBound = this->fieldfunc_m(t, i, externalE, externalB, args ...);
double qtom = bunch->Q[i] / (bunch->M[i] * mass_coeff); // m^2/s^2/GV
- double tempgamma = sqrt(1 + (y[3] * y[3] + y[4] * y[4] + y[5] * y[5]));
-
- yp[0] = c_mtns / tempgamma * y[3]; // [m/ns]
- yp[1] = c_mtns / tempgamma * y[4]; // [m/ns]
- yp[2] = c_mtns / tempgamma * y[5]; // [m/ns]
+ double tempgamma = std::sqrt(1 + (y[3] * y[3] + y[4] * y[4] + y[5] * y[5]));
-/*
- yp[0] = c_mmtns / tempgamma * y[3]; // [mm/ns]
- yp[1] = c_mmtns / tempgamma * y[4]; // [mm/ns]
- yp[2] = c_mmtns / tempgamma * y[5]; // [mm/ns]
-*/
+ yp[0] = Physics::c / tempgamma * y[3];
+ yp[1] = Physics::c / tempgamma * y[4];
+ yp[2] = Physics::c / tempgamma * y[5];
yp[3] = (externalE(0) / Physics::c + (externalB(2) * y[4] - externalB(1) * y[5]) / tempgamma) * qtom; // [1/ns]
yp[4] = (externalE(1) / Physics::c - (externalB(2) * y[3] - externalB(0) * y[5]) / tempgamma) * qtom; // [1/ns];
yp[5] = (externalE(2) / Physics::c + (externalB(1) * y[3] - externalB(0) * y[4]) / tempgamma) * qtom; // [1/ns];
+ yp[3] /= Units::ns2s;
+ yp[4] /= Units::ns2s;
+ yp[5] /= Units::ns2s;
+
return outOfBound;
}
@@ -146,7 +149,7 @@ void RK4<FieldFunction, Arguments ...>::copyFrom(Vector_t& R,
Vector_t& P,
double* x) const
{
- for(int j = 0; j < 3; j++) {
+ for (int j = 0; j < 3; j++) {
R(j) = x[j]; // [x,y,z] (mm)
P(j) = x[j+3]; // [px,py,pz] (beta*gamma)
}