Commit 3408e6fb authored by snuverink_j's avatar snuverink_j
Browse files

Merge branch 'use-std__round' into 'master'

Use std::round(x) instead of std::floor(x+0.5)

See merge request !297
parents 80ed400a a965ca09
...@@ -15,6 +15,7 @@ ...@@ -15,6 +15,7 @@
#include <boost/filesystem.hpp> #include <boost/filesystem.hpp>
#include <cmath>
#include <iostream> #include <iostream>
#include <limits> #include <limits>
...@@ -72,7 +73,7 @@ OrbitThreader::OrbitThreader(const PartData &ref, ...@@ -72,7 +73,7 @@ OrbitThreader::OrbitThreader(const PartData &ref,
logger_m.open(fileName, std::ios_base::app); logger_m.open(fileName, std::ios_base::app);
} }
loggingFrequency_m = std::max(1.0, floor(1e-11 / std::abs(dt_m) + 0.5)); loggingFrequency_m = std::max(1.0, std::round(1e-11 / std::abs(dt_m)));
} else { } else {
loggingFrequency_m = std::numeric_limits<size_t>::max(); loggingFrequency_m = std::numeric_limits<size_t>::max();
} }
...@@ -455,4 +456,4 @@ double OrbitThreader::computeMaximalImplicitDrift() { ...@@ -455,4 +456,4 @@ double OrbitThreader::computeMaximalImplicitDrift() {
maxDrift = std::min(maxIntegSteps_m * std::abs(dt_m) * Physics::c, maxDrift); maxDrift = std::min(maxIntegSteps_m * std::abs(dt_m) * Physics::c, maxDrift);
return maxDrift; return maxDrift;
} }
\ No newline at end of file
...@@ -28,6 +28,8 @@ ...@@ -28,6 +28,8 @@
#include "gsl/gsl_interp.h" #include "gsl/gsl_interp.h"
#include "gsl/gsl_spline.h" #include "gsl/gsl_spline.h"
#include <cmath>
#include <iostream> #include <iostream>
#include <fstream> #include <fstream>
...@@ -616,7 +618,7 @@ double RFCavity::getAutoPhaseEstimateFallback(double E0, double t0, double q, do ...@@ -616,7 +618,7 @@ double RFCavity::getAutoPhaseEstimateFallback(double E0, double t0, double q, do
dphi = dphi / 17.5; dphi = dphi / 17.5;
} }
phimax = phimax - floor(phimax / Physics::two_pi + 0.5) * Physics::two_pi; phimax = phimax - std::round(phimax / Physics::two_pi) * Physics::two_pi;
phimax = fmod(phimax, Physics::two_pi); phimax = fmod(phimax, Physics::two_pi);
const int prevPrecision = Ippl::Info->precision(8); const int prevPrecision = Ippl::Info->precision(8);
...@@ -711,7 +713,7 @@ double RFCavity::getAutoPhaseEstimate(const double &E0, const double &t0, const ...@@ -711,7 +713,7 @@ double RFCavity::getAutoPhaseEstimate(const double &E0, const double &t0, const
return tmp_phi; return tmp_phi;
} }
phi = tmp_phi - floor(tmp_phi / Physics::two_pi + 0.5) * Physics::two_pi; phi = tmp_phi - std::round(tmp_phi / Physics::two_pi) * Physics::two_pi;
for(unsigned int i = 1; i < N; ++ i) { for(unsigned int i = 1; i < N; ++ i) {
E[i] = E[i - 1]; E[i] = E[i - 1];
......
...@@ -25,6 +25,8 @@ ...@@ -25,6 +25,8 @@
#include "gsl/gsl_interp.h" #include "gsl/gsl_interp.h"
#include "gsl/gsl_spline.h" #include "gsl/gsl_spline.h"
#include <cmath>
#include <iostream> #include <iostream>
#include <fstream> #include <fstream>
...@@ -116,14 +118,14 @@ void TravelingWave::addKR(int i, double t, Vector_t &K) { ...@@ -116,14 +118,14 @@ void TravelingWave::addKR(int i, double t, Vector_t &K) {
wtf = frequency_m * t + phase_m; wtf = frequency_m * t + phase_m;
CoreFieldmap_m->getFieldstrength(tmpA0, tmpE, tmpB); CoreFieldmap_m->getFieldstrength(tmpA0, tmpE, tmpB);
CoreFieldmap_m->getFieldDerivative(tmpA0, tmpE_diff, tmpB_diff, zdir); CoreFieldmap_m->getFieldDerivative(tmpA0, tmpE_diff, tmpB_diff, zdir);
k = scale_m * (tmpE_diff(2) * cos(wtf) - b * frequency_m * tmpE(2) * sin(wtf) / Physics::c); k = scale_m * (tmpE_diff(2) * std::cos(wtf) - b * frequency_m * tmpE(2) * std::sin(wtf) / Physics::c);
} else if(tmpA0(2) < startExitField_m) { } else if(tmpA0(2) < startExitField_m) {
wtf = frequency_m * t + phaseCore1_m; wtf = frequency_m * t + phaseCore1_m;
tmpA0(2) -= startCoreField_m; tmpA0(2) -= startCoreField_m;
const double z = tmpA0(2); const double z = tmpA0(2);
tmpA0(2) = tmpA0(2) - PeriodLength_m * floor(tmpA0(2) / PeriodLength_m); tmpA0(2) = tmpA0(2) - PeriodLength_m * std::floor(tmpA0(2) / PeriodLength_m);
tmpA0(2) += startCoreField_m; tmpA0(2) += startCoreField_m;
tmpE = 0.0; tmpE = 0.0;
tmpB = 0.0; tmpB = 0.0;
...@@ -132,11 +134,11 @@ void TravelingWave::addKR(int i, double t, Vector_t &K) { ...@@ -132,11 +134,11 @@ void TravelingWave::addKR(int i, double t, Vector_t &K) {
CoreFieldmap_m->getFieldstrength(tmpA0, tmpE, tmpB); CoreFieldmap_m->getFieldstrength(tmpA0, tmpE, tmpB);
CoreFieldmap_m->getFieldDerivative(tmpA0, tmpE_diff, tmpB_diff, zdir); CoreFieldmap_m->getFieldDerivative(tmpA0, tmpE_diff, tmpB_diff, zdir);
k = scaleCore_m * (tmpE_diff(2) * cos(wtf) - b * frequency_m * tmpE(2) * sin(wtf) / Physics::c); k = scaleCore_m * (tmpE_diff(2) * std::cos(wtf) - b * frequency_m * tmpE(2) * std::sin(wtf) / Physics::c);
wtf = frequency_m * t + phaseCore2_m; wtf = frequency_m * t + phaseCore2_m;
tmpA0(2) = z + CellLength_m; tmpA0(2) = z + CellLength_m;
tmpA0(2) = tmpA0(2) - PeriodLength_m * floor(tmpA0(2) / PeriodLength_m); tmpA0(2) = tmpA0(2) - PeriodLength_m * std::floor(tmpA0(2) / PeriodLength_m);
tmpA0(2) += startCoreField_m; tmpA0(2) += startCoreField_m;
tmpE = 0.0; tmpE = 0.0;
tmpB = 0.0; tmpB = 0.0;
...@@ -145,7 +147,7 @@ void TravelingWave::addKR(int i, double t, Vector_t &K) { ...@@ -145,7 +147,7 @@ void TravelingWave::addKR(int i, double t, Vector_t &K) {
CoreFieldmap_m->getFieldstrength(tmpA0, tmpE, tmpB); CoreFieldmap_m->getFieldstrength(tmpA0, tmpE, tmpB);
CoreFieldmap_m->getFieldDerivative(tmpA0, tmpE_diff, tmpB_diff, zdir); CoreFieldmap_m->getFieldDerivative(tmpA0, tmpE_diff, tmpB_diff, zdir);
k += scaleCore_m * (tmpE_diff(2) * cos(wtf) - b * frequency_m * tmpE(2) * sin(wtf) / Physics::c); k += scaleCore_m * (tmpE_diff(2) * std::cos(wtf) - b * frequency_m * tmpE(2) * std::sin(wtf) / Physics::c);
} else { } else {
...@@ -158,7 +160,7 @@ void TravelingWave::addKR(int i, double t, Vector_t &K) { ...@@ -158,7 +160,7 @@ void TravelingWave::addKR(int i, double t, Vector_t &K) {
CoreFieldmap_m->getFieldstrength(tmpA0, tmpE, tmpB); CoreFieldmap_m->getFieldstrength(tmpA0, tmpE, tmpB);
CoreFieldmap_m->getFieldDerivative(tmpA0, tmpE_diff, tmpB_diff, zdir); CoreFieldmap_m->getFieldDerivative(tmpA0, tmpE_diff, tmpB_diff, zdir);
k = scale_m * (tmpE_diff(2) * cos(wtf) - b * frequency_m * tmpE(2) * sin(wtf) / Physics::c); k = scale_m * (tmpE_diff(2) * std::cos(wtf) - b * frequency_m * tmpE(2) * std::sin(wtf) / Physics::c);
} }
...@@ -229,19 +231,19 @@ bool TravelingWave::apply(const Vector_t &R, const Vector_t &/*P*/, const double ...@@ -229,19 +231,19 @@ bool TravelingWave::apply(const Vector_t &R, const Vector_t &/*P*/, const double
if (tmpR(2) < startCoreField_m) { if (tmpR(2) < startCoreField_m) {
if (!CoreFieldmap_m->isInside(tmpR)) return true; if (!CoreFieldmap_m->isInside(tmpR)) return true;
tmpcos = (scale_m + scaleError_m) * cos(frequency_m * t + phase_m + phaseError_m); tmpcos = (scale_m + scaleError_m) * std::cos(frequency_m * t + phase_m + phaseError_m);
tmpsin = -(scale_m + scaleError_m) * sin(frequency_m * t + phase_m + phaseError_m); tmpsin = -(scale_m + scaleError_m) * std::sin(frequency_m * t + phase_m + phaseError_m);
} else if (tmpR(2) < startExitField_m) { } else if (tmpR(2) < startExitField_m) {
Vector_t tmpE2(0.0, 0.0, 0.0), tmpB2(0.0, 0.0, 0.0); Vector_t tmpE2(0.0, 0.0, 0.0), tmpB2(0.0, 0.0, 0.0);
tmpR(2) -= startCoreField_m; tmpR(2) -= startCoreField_m;
const double z = tmpR(2); const double z = tmpR(2);
tmpR(2) = tmpR(2) - PeriodLength_m * floor(tmpR(2) / PeriodLength_m); tmpR(2) = tmpR(2) - PeriodLength_m * std::floor(tmpR(2) / PeriodLength_m);
tmpR(2) += startCoreField_m; tmpR(2) += startCoreField_m;
if (!CoreFieldmap_m->isInside(tmpR)) return true; if (!CoreFieldmap_m->isInside(tmpR)) return true;
tmpcos = (scaleCore_m + scaleCoreError_m) * cos(frequency_m * t + phaseCore1_m + phaseError_m); tmpcos = (scaleCore_m + scaleCoreError_m) * std::cos(frequency_m * t + phaseCore1_m + phaseError_m);
tmpsin = -(scaleCore_m + scaleCoreError_m) * sin(frequency_m * t + phaseCore1_m + phaseError_m); tmpsin = -(scaleCore_m + scaleCoreError_m) * std::sin(frequency_m * t + phaseCore1_m + phaseError_m);
CoreFieldmap_m->getFieldstrength(tmpR, tmpE, tmpB); CoreFieldmap_m->getFieldstrength(tmpR, tmpE, tmpB);
E += tmpcos * tmpE; E += tmpcos * tmpE;
B += tmpsin * tmpB; B += tmpsin * tmpB;
...@@ -250,17 +252,17 @@ bool TravelingWave::apply(const Vector_t &R, const Vector_t &/*P*/, const double ...@@ -250,17 +252,17 @@ bool TravelingWave::apply(const Vector_t &R, const Vector_t &/*P*/, const double
tmpB = 0.0; tmpB = 0.0;
tmpR(2) = z + CellLength_m; tmpR(2) = z + CellLength_m;
tmpR(2) = tmpR(2) - PeriodLength_m * floor(tmpR(2) / PeriodLength_m); tmpR(2) = tmpR(2) - PeriodLength_m * std::floor(tmpR(2) / PeriodLength_m);
tmpR(2) += startCoreField_m; tmpR(2) += startCoreField_m;
tmpcos = (scaleCore_m + scaleCoreError_m) * cos(frequency_m * t + phaseCore2_m + phaseError_m); tmpcos = (scaleCore_m + scaleCoreError_m) * std::cos(frequency_m * t + phaseCore2_m + phaseError_m);
tmpsin = -(scaleCore_m + scaleCoreError_m) * sin(frequency_m * t + phaseCore2_m + phaseError_m); tmpsin = -(scaleCore_m + scaleCoreError_m) * std::sin(frequency_m * t + phaseCore2_m + phaseError_m);
} else { } else {
tmpR(2) -= mappedStartExitField_m; tmpR(2) -= mappedStartExitField_m;
if (!CoreFieldmap_m->isInside(tmpR)) return true; if (!CoreFieldmap_m->isInside(tmpR)) return true;
tmpcos = (scale_m + scaleError_m) * cos(frequency_m * t + phaseExit_m + phaseError_m); tmpcos = (scale_m + scaleError_m) * std::cos(frequency_m * t + phaseExit_m + phaseError_m);
tmpsin = -(scale_m + scaleError_m) * sin(frequency_m * t + phaseExit_m + phaseError_m); tmpsin = -(scale_m + scaleError_m) * std::sin(frequency_m * t + phaseExit_m + phaseError_m);
} }
CoreFieldmap_m->getFieldstrength(tmpR, tmpE, tmpB); CoreFieldmap_m->getFieldstrength(tmpR, tmpE, tmpB);
...@@ -281,19 +283,19 @@ bool TravelingWave::applyToReferenceParticle(const Vector_t &R, const Vector_t & ...@@ -281,19 +283,19 @@ bool TravelingWave::applyToReferenceParticle(const Vector_t &R, const Vector_t &
if (tmpR(2) < startCoreField_m) { if (tmpR(2) < startCoreField_m) {
if (!CoreFieldmap_m->isInside(tmpR)) return true; if (!CoreFieldmap_m->isInside(tmpR)) return true;
tmpcos = scale_m * cos(frequency_m * t + phase_m); tmpcos = scale_m * std::cos(frequency_m * t + phase_m);
tmpsin = -scale_m * sin(frequency_m * t + phase_m); tmpsin = -scale_m * std::sin(frequency_m * t + phase_m);
} else if (tmpR(2) < startExitField_m) { } else if (tmpR(2) < startExitField_m) {
Vector_t tmpE2(0.0, 0.0, 0.0), tmpB2(0.0, 0.0, 0.0); Vector_t tmpE2(0.0, 0.0, 0.0), tmpB2(0.0, 0.0, 0.0);
tmpR(2) -= startCoreField_m; tmpR(2) -= startCoreField_m;
const double z = tmpR(2); const double z = tmpR(2);
tmpR(2) = tmpR(2) - PeriodLength_m * floor(tmpR(2) / PeriodLength_m); tmpR(2) = tmpR(2) - PeriodLength_m * std::floor(tmpR(2) / PeriodLength_m);
tmpR(2) += startCoreField_m; tmpR(2) += startCoreField_m;
if (!CoreFieldmap_m->isInside(tmpR)) return true; if (!CoreFieldmap_m->isInside(tmpR)) return true;
tmpcos = scaleCore_m * cos(frequency_m * t + phaseCore1_m); tmpcos = scaleCore_m * std::cos(frequency_m * t + phaseCore1_m);
tmpsin = -scaleCore_m * sin(frequency_m * t + phaseCore1_m); tmpsin = -scaleCore_m * std::sin(frequency_m * t + phaseCore1_m);
CoreFieldmap_m->getFieldstrength(tmpR, tmpE, tmpB); CoreFieldmap_m->getFieldstrength(tmpR, tmpE, tmpB);
E += tmpcos * tmpE; E += tmpcos * tmpE;
B += tmpsin * tmpB; B += tmpsin * tmpB;
...@@ -302,18 +304,18 @@ bool TravelingWave::applyToReferenceParticle(const Vector_t &R, const Vector_t & ...@@ -302,18 +304,18 @@ bool TravelingWave::applyToReferenceParticle(const Vector_t &R, const Vector_t &
tmpB = 0.0; tmpB = 0.0;
tmpR(2) = z + CellLength_m; tmpR(2) = z + CellLength_m;
tmpR(2) = tmpR(2) - PeriodLength_m * floor(tmpR(2) / PeriodLength_m); tmpR(2) = tmpR(2) - PeriodLength_m * std::floor(tmpR(2) / PeriodLength_m);
tmpR(2) += startCoreField_m; tmpR(2) += startCoreField_m;
tmpcos = scaleCore_m * cos(frequency_m * t + phaseCore2_m); tmpcos = scaleCore_m * std::cos(frequency_m * t + phaseCore2_m);
tmpsin = -scaleCore_m * sin(frequency_m * t + phaseCore2_m); tmpsin = -scaleCore_m * std::sin(frequency_m * t + phaseCore2_m);
} else { } else {
tmpR(2) -= mappedStartExitField_m; tmpR(2) -= mappedStartExitField_m;
if (!CoreFieldmap_m->isInside(tmpR)) return true; if (!CoreFieldmap_m->isInside(tmpR)) return true;
tmpcos = scale_m * cos(frequency_m * t + phaseExit_m); tmpcos = scale_m * std::cos(frequency_m * t + phaseExit_m);
tmpsin = -scale_m * sin(frequency_m * t + phaseExit_m); tmpsin = -scale_m * std::sin(frequency_m * t + phaseExit_m);
} }
CoreFieldmap_m->getFieldstrength(tmpR, tmpE, tmpB); CoreFieldmap_m->getFieldstrength(tmpR, tmpE, tmpB);
...@@ -373,11 +375,11 @@ void TravelingWave::initialise(PartBunchBase<double, 3> *bunch, double &startFie ...@@ -373,11 +375,11 @@ void TravelingWave::initialise(PartBunchBase<double, 3> *bunch, double &startFie
endField = startField + startExitField_m + PeriodLength_m / 2.0; endField = startField + startExitField_m + PeriodLength_m / 2.0;
length_m = endField - startField; length_m = endField - startField;
scaleCore_m = scale_m / sin(2.0 * pi * Mode_m); scaleCore_m = scale_m / std::sin(2.0 * pi * Mode_m);
scaleCoreError_m = scaleError_m / sin(2.0 * pi * Mode_m); scaleCoreError_m = scaleError_m / std::sin(2.0 * pi * Mode_m);
phaseCore1_m = phase_m + pi * Mode_m / 2.0; phaseCore1_m = phase_m + pi * Mode_m / 2.0;
phaseCore2_m = phase_m + pi * Mode_m * 1.5; phaseCore2_m = phase_m + pi * Mode_m * 1.5;
phaseExit_m = phase_m - 2.0 * pi * ((NumCells_m - 1) * Mode_m - floor((NumCells_m - 1) * Mode_m)); phaseExit_m = phase_m - 2.0 * pi * ((NumCells_m - 1) * Mode_m - std::floor((NumCells_m - 1) * Mode_m));
} else { } else {
endField = startField - 1e-3; endField = startField - 1e-3;
...@@ -438,10 +440,10 @@ double TravelingWave::getAutoPhaseEstimate(const double &E0, const double &t0, c ...@@ -438,10 +440,10 @@ double TravelingWave::getAutoPhaseEstimate(const double &E0, const double &t0, c
CoreFieldmap_m->getOnaxisEz(F); CoreFieldmap_m->getOnaxisEz(F);
if(F.size() == 0) return 0.0; if(F.size() == 0) return 0.0;
N1 = static_cast<int>(floor(F.size() / 4.)) + 1; N1 = static_cast<int>(std::floor(F.size() / 4.)) + 1;
N2 = F.size() - 2 * N1 + 1; N2 = F.size() - 2 * N1 + 1;
N3 = 2 * N1 + static_cast<int>(floor((NumCells_m - 1) * N2 * Mode_m)) - 1; N3 = 2 * N1 + static_cast<int>(std::floor((NumCells_m - 1) * N2 * Mode_m)) - 1;
N4 = static_cast<int>(floor(0.5 + N2 * Mode_m)); N4 = static_cast<int>(std::round(N2 * Mode_m));
Dz = F[N1 + N2].first - F[N1].first; Dz = F[N1 + N2].first - F[N1].first;
t.resize(N3, t0); t.resize(N3, t0);
...@@ -454,7 +456,7 @@ double TravelingWave::getAutoPhaseEstimate(const double &E0, const double &t0, c ...@@ -454,7 +456,7 @@ double TravelingWave::getAutoPhaseEstimate(const double &E0, const double &t0, c
} }
for(int i = N1; i < N3 - N1 + 1; ++ i) { for(int i = N1; i < N3 - N1 + 1; ++ i) {
int I = (i - N1) % N2 + N1; int I = (i - N1) % N2 + N1;
double Z = (F[I].first + F[I - 1].first) / 2. + floor((i - N1) / N2) * Dz; double Z = (F[I].first + F[I - 1].first) / 2. + std::floor((i - N1) / N2) * Dz;
E[i] = E0 + Z * scaleCore_m / mass; E[i] = E0 + Z * scaleCore_m / mass;
E2[i] = E[i]; E2[i] = E[i];
} }
...@@ -490,11 +492,11 @@ double TravelingWave::getAutoPhaseEstimate(const double &E0, const double &t0, c ...@@ -490,11 +492,11 @@ double TravelingWave::getAutoPhaseEstimate(const double &E0, const double &t0, c
} }
if(std::abs(B) > 0.0000001) { if(std::abs(B) > 0.0000001) {
tmp_phi = atan(A / B); tmp_phi = std::atan(A / B);
} else { } else {
tmp_phi = Physics::pi / 2; tmp_phi = Physics::pi / 2;
} }
if(q * (A * sin(tmp_phi) + B * cos(tmp_phi)) < 0) { if(q * (A * std::sin(tmp_phi) + B * std::cos(tmp_phi)) < 0) {
tmp_phi += Physics::pi; tmp_phi += Physics::pi;
} }
...@@ -518,7 +520,7 @@ double TravelingWave::getAutoPhaseEstimate(const double &E0, const double &t0, c ...@@ -518,7 +520,7 @@ double TravelingWave::getAutoPhaseEstimate(const double &E0, const double &t0, c
<< "Ekin= " << E[N3 - 1] << " MeV" << std::setprecision(prevPrecision) << "\n" << endl); << "Ekin= " << E[N3 - 1] << " MeV" << std::setprecision(prevPrecision) << "\n" << endl);
return tmp_phi; return tmp_phi;
} }
phi = tmp_phi - floor(tmp_phi / Physics::two_pi + 0.5) * Physics::two_pi; phi = tmp_phi - std::round(tmp_phi / Physics::two_pi) * Physics::two_pi;
for(int i = 1; i < N1; ++ i) { for(int i = 1; i < N1; ++ i) {
...@@ -596,13 +598,13 @@ std::pair<double, double> TravelingWave::trackOnAxisParticle(const double &p0, ...@@ -596,13 +598,13 @@ std::pair<double, double> TravelingWave::trackOnAxisParticle(const double &p0,
if (out) *out << std::setw(18) << z if (out) *out << std::setw(18) << z
<< std::setw(18) << Util::getEnergy(p, mass) << std::setw(18) << Util::getEnergy(p, mass)
<< std::endl; << std::endl;
double dz = 0.5 * p / sqrt(1.0 + p * p) * cdt; double dz = 0.5 * p / std::sqrt(1.0 + p * p) * cdt;
while(z + dz < startCoreField_m + zbegin) { while(z + dz < startCoreField_m + zbegin) {
z += dz; z += dz;
double ez = scale_m / Ezmax * cos(phase) * gsl_spline_eval(onAxisInterpolants, z, onAxisAccel); double ez = scale_m / Ezmax * std::cos(phase) * gsl_spline_eval(onAxisInterpolants, z, onAxisAccel);
p += ez * q * cdt / mass; p += ez * q * cdt / mass;
dz = 0.5 * p / sqrt(1.0 + p * p) * cdt; dz = 0.5 * p / std::sqrt(1.0 + p * p) * cdt;
z += 0.5 * p / sqrt(1.0 + p * p) * cdt; z += 0.5 * p / std::sqrt(1.0 + p * p) * cdt;
phase += dphi; phase += dphi;
t += dt; t += dt;
} }
...@@ -612,19 +614,19 @@ std::pair<double, double> TravelingWave::trackOnAxisParticle(const double &p0, ...@@ -612,19 +614,19 @@ std::pair<double, double> TravelingWave::trackOnAxisParticle(const double &p0,
z += dz; z += dz;
double tmpz = z - (startCoreField_m + zbegin); double tmpz = z - (startCoreField_m + zbegin);
tmpz -= PeriodLength_m * floor(tmpz / PeriodLength_m); tmpz -= PeriodLength_m * std::floor(tmpz / PeriodLength_m);
tmpz += startCoreField_m + zbegin; tmpz += startCoreField_m + zbegin;
double ez = scaleCore_m / Ezmax * cos(phase) * gsl_spline_eval(onAxisInterpolants, tmpz, onAxisAccel); double ez = scaleCore_m / Ezmax * std::cos(phase) * gsl_spline_eval(onAxisInterpolants, tmpz, onAxisAccel);
tmpz = z - (startCoreField_m + zbegin) + CellLength_m; tmpz = z - (startCoreField_m + zbegin) + CellLength_m;
tmpz -= PeriodLength_m * floor(tmpz / PeriodLength_m); tmpz -= PeriodLength_m * std::floor(tmpz / PeriodLength_m);
tmpz += startCoreField_m + zbegin; tmpz += startCoreField_m + zbegin;
ez += scaleCore_m / Ezmax * cos(phase2) * gsl_spline_eval(onAxisInterpolants, tmpz, onAxisAccel); ez += scaleCore_m / Ezmax * std::cos(phase2) * gsl_spline_eval(onAxisInterpolants, tmpz, onAxisAccel);
p += ez * q * cdt / mass; p += ez * q * cdt / mass;
dz = 0.5 * p / sqrt(1.0 + p * p) * cdt; dz = 0.5 * p / std::sqrt(1.0 + p * p) * cdt;
z += dz; z += dz;
phase += dphi; phase += dphi;
phase2 += dphi; phase2 += dphi;
...@@ -638,11 +640,11 @@ std::pair<double, double> TravelingWave::trackOnAxisParticle(const double &p0, ...@@ -638,11 +640,11 @@ std::pair<double, double> TravelingWave::trackOnAxisParticle(const double &p0,
while(z + dz < startExitField_m + 0.5 * PeriodLength_m + zbegin) { while(z + dz < startExitField_m + 0.5 * PeriodLength_m + zbegin) {
z += dz; z += dz;
double tmpz = z - (startExitField_m + zbegin); double tmpz = z - (startExitField_m + zbegin);
tmpz -= PeriodLength_m * floor(tmpz / PeriodLength_m); tmpz -= PeriodLength_m * std::floor(tmpz / PeriodLength_m);
tmpz += 3.0 * PeriodLength_m / 2.0; tmpz += 3.0 * PeriodLength_m / 2.0;
double ez = scale_m / Ezmax * cos(phase) * gsl_spline_eval(onAxisInterpolants, tmpz, onAxisAccel); double ez = scale_m / Ezmax * std::cos(phase) * gsl_spline_eval(onAxisInterpolants, tmpz, onAxisAccel);
p += ez * q * cdt / mass; p += ez * q * cdt / mass;
dz = 0.5 * p / sqrt(1.0 + p * p) * cdt; dz = 0.5 * p / std::sqrt(1.0 + p * p) * cdt;
z += dz; z += dz;
phase += dphi; phase += dphi;
t += dt; t += dt;
...@@ -651,7 +653,7 @@ std::pair<double, double> TravelingWave::trackOnAxisParticle(const double &p0, ...@@ -651,7 +653,7 @@ std::pair<double, double> TravelingWave::trackOnAxisParticle(const double &p0,
gsl_spline_free(onAxisInterpolants); gsl_spline_free(onAxisInterpolants);
gsl_interp_accel_free(onAxisAccel); gsl_interp_accel_free(onAxisAccel);
const double beta = sqrt(1. - 1 / (p * p + 1.)); const double beta = std::sqrt(1. - 1 / (p * p + 1.));
const double tErr = (z - (startExitField_m + 0.5 * PeriodLength_m + zbegin)) / (Physics::c * beta); const double tErr = (z - (startExitField_m + 0.5 * PeriodLength_m + zbegin)) / (Physics::c * beta);
return std::pair<double, double>(p, t - tErr); return std::pair<double, double>(p, t - tErr);
......
...@@ -3,6 +3,7 @@ ...@@ -3,6 +3,7 @@
#include "H5hut.h" #include "H5hut.h"
#include "Physics/Physics.h" #include "Physics/Physics.h"
#include <cmath>
#include <fstream> #include <fstream>
#include <ios> #include <ios>
...@@ -97,7 +98,7 @@ void FM3DH5Block_nonscale::readMap() { ...@@ -97,7 +98,7 @@ void FM3DH5Block_nonscale::readMap() {
long field_size = 0; long field_size = 0;
int Nnodes = Ippl::getNodes();//min(20, Ippl::getNodes()); int Nnodes = Ippl::getNodes();//min(20, Ippl::getNodes());
int Nz_avrg = static_cast<int>(floor(0.5 + num_gridpz_m / Nnodes)); int Nz_avrg = static_cast<int>(std::round((double)num_gridpz_m / Nnodes));
int Nz_diff = Nz_avrg * Nnodes - num_gridpz_m; int Nz_diff = Nz_avrg * Nnodes - num_gridpz_m;
int signNz = Nz_diff > 0 ? 1 : -1; int signNz = Nz_diff > 0 ? 1 : -1;
int *Nz_read_start = new int[Ippl::getNodes() + 1]; int *Nz_read_start = new int[Ippl::getNodes() + 1];
......
...@@ -127,8 +127,8 @@ void FM3DMagnetoStatic::readMap() { ...@@ -127,8 +127,8 @@ void FM3DMagnetoStatic::readMap() {
if (normalize_m) { if (normalize_m) {
double Bymax = 0.0; double Bymax = 0.0;
// find maximum field // find maximum field
unsigned int centerX = static_cast<unsigned int>(std::floor(-xbegin_m / hx_m + 0.5)); unsigned int centerX = static_cast<unsigned int>(std::round(-xbegin_m / hx_m));
unsigned int centerY = static_cast<unsigned int>(std::floor(-ybegin_m / hy_m + 0.5)); unsigned int centerY = static_cast<unsigned int>(std::round(-ybegin_m / hy_m));
for(unsigned int k = 0; k < num_gridpz_m; ++ k) { for(unsigned int k = 0; k < num_gridpz_m; ++ k) {
double By = FieldstrengthBy_m[getIndex(centerX, centerY, k)]; double By = FieldstrengthBy_m[getIndex(centerX, centerY, k)];
if(std::abs(By) > Bymax) { if(std::abs(By) > Bymax) {
......
...@@ -3,6 +3,7 @@ ...@@ -3,6 +3,7 @@
#include "Utilities/GeneralClassicException.h" #include "Utilities/GeneralClassicException.h"
#include "Utilities/Util.h" #include "Utilities/Util.h"
#include <cmath>
#include <fstream> #include <fstream>
#include <ios> #include <ios>
#include <algorithm> #include <algorithm>
...@@ -121,7 +122,7 @@ void FM3DMagnetoStaticExtended::readMap() { ...@@ -121,7 +122,7 @@ void FM3DMagnetoStaticExtended::readMap() {
double Bymax = 0.0; double Bymax = 0.0;
// find maximum field // find maximum field
unsigned int centerX = static_cast<unsigned int>(std::floor(-xbegin_m / hx_m + 0.5)); unsigned int centerX = static_cast<unsigned int>(std::round(-xbegin_m / hx_m));
for(unsigned int k = 0; k < num_gridpz_m; ++ k) { for(unsigned int k = 0; k < num_gridpz_m; ++ k) {
double By = FieldstrengthBy_m[getIndex(centerX, 0, k)]; double By = FieldstrengthBy_m[getIndex(centerX, 0, k)];
if(std::abs(By) > std::abs(Bymax)) { if(std::abs(By) > std::abs(Bymax)) {
......
...@@ -3,6 +3,7 @@ ...@@ -3,6 +3,7 @@
#include "H5hut.h" #include "H5hut.h"
#include "Physics/Physics.h" #include "Physics/Physics.h"
#include <cmath>
#include <fstream> #include <fstream>
#include <ios> #include <ios>
...@@ -100,7 +101,7 @@ void FM3DMagnetoStaticH5Block::readMap() { ...@@ -100,7 +101,7 @@ void FM3DMagnetoStaticH5Block::readMap() {
long field_size = 0; long field_size = 0;
int Nnodes = Ippl::getNodes();//min(20, Ippl::getNodes()); int Nnodes = Ippl::getNodes();//min(20, Ippl::getNodes());
int Nz_avrg = static_cast<int>(floor(0.5 + num_gridpz_m / Nnodes)); int Nz_avrg = static_cast<int>(std::round((double)num_gridpz_m / Nnodes));